########################################################
## ANOVARQ
## Type III ANOVA breakdown for quantile regression models
##
## Ben Meuleman, 2025
########################################################

## LOAD REQUIRED PACKAGES
library(quantreg)
library(insight)
library(psych)
library(foreach)

## QUANTILE ANOVA FUNCTION
AnovaRQ <- function(model,test="Wald",type="within") {
  
  ## MODEL INFORMATION AND REFIT
  form <- formula(model)
  formcharacter <- paste(as.character(form)[2],as.character(form)[1],as.character(form)[3],sep="")
  tau <- model$tau
  data <- get_data(model)
  predvars <- find_predictors(model)$conditional

  mm <- data.frame(model.matrix(lm(form,data=data))[,-1],Y=data[,find_response(model)])
  source.model <- rq(Y~.,data=mm,tau=tau)

  ## EFFECT AND PARAMETER IDENTIFICATION
  effects <- c(fparse(form)$x,unlist(lapply(fparse(form)$prod,paste,collapse=":")))
  parms <- find_parameters(model)$conditional[-1]
  effect.order <- attr(model.matrix(lm(form,data=data)),"assign")[-1]

  ## WITHIN-QUANTILE ANOVA
  if(type=="within") {
    anova.out <- foreach(t=1:length(tau),.combine="rbind") %:%
      foreach(i=1:length(effects),.combine="rbind") %do% {
        cureff <- effects[i]
        curorder <- !c(effect.order==i,TRUE)
        newform <- as.formula(paste("Y~",paste(colnames(mm)[curorder],collapse="+"),sep=""))
        source.model <- rq(Y~.,data=mm,tau=tau[t])
        cur.model <- rq(newform,data=mm,tau=tau[t])
        anova(cur.model,source.model,test=test)$table
      }
    out <- data.frame(tau=rep(tau,each=length(effects)),Effect=rep(effects,times=length(tau)),anova.out)
    out$sig <- cut(out$pvalue,breaks=c(0,0.001,0.01,0.05,0.1,1),labels=c("***","**","*",".",""),include.lowest=TRUE)
    colnames(out) <- c("tau","effect","DF1","DF2","F-value","p-value","sig")
    options(contrasts=c("contr.treatment","contr.treatment"))

    ## OUTPUT
    cat("Quantile regression Type III ANOVA breakdown\n\nWithin-quantile",test,"tests:\n\n")
    print(out)
    cat("------------------------\nSignif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1\n------------------------\n\n")
    out
  }

  ##BETWEEN-QAUNTILE ANOVA
  if(type=="between") {
    tau.pairs <- combn(tau,2)
    anova.out <- foreach(i=1:ncol(tau.pairs),.combine="rbind") %do% {
      cur.model <- rq(form,data=data,tau=tau.pairs[,i])
      anova(cur.model,test=test,joint=FALSE)$table
    }
    out <- data.frame(tau.pair=rep(apply(tau.pairs,2,paste,collapse="-"),each=length(parms)),Effect=rep(parms,times=length(tau)),anova.out)
    out$sig <- cut(out$pvalue,breaks=c(0,0.001,0.01,0.05,0.1,1),labels=c("***","**","*",".",""),include.lowest=TRUE)
    rownames(out) <- 1:nrow(out)
    colnames(out) <- c("tau.pair","parameter","DF1","DF2","F-value","p-value","sig")
    options(contrasts=c("contr.treatment","contr.treatment"))

    ## OUTPUT
    cat("Quantile regression Type III ANOVA breakdown\n\nBetween-quantile",test,"tests:\n\n")
    print(out)
    cat("------------------------\nSignif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1\n------------------------\n\n")
    invisible(out)
  }
}