########################################################
## ANOVAHC
## Heteroscedastic Type II ANOVA for regression models
##
## Ben Meuleman, 2025
########################################################

## REQUIRED PACKAGES
library(clubSandwich)
library(foreach)
library(insight)

## FUNCTION
AnovaHC <- function(model,random=NULL) {

  ## MODEL INFORMATION
  modinfo <- terms(model)
  data <- get_data(model)  

  #predvars <- attr(modinfo,"term.labels")[attr(modinfo,"order")==1]
  predvars <- find_predictors(model)$conditional
  classes <- sapply(data[,predvars],is.numeric)

  options(contrasts=c("contr.sum","contr.sum"))
  tempdata <- data
  tempdata[,predvars[which(classes)]] <- scale(tempdata[,predvars[which(classes)]])

  tempmod <- lm(formula(modinfo),data=tempdata)
  effects <- attr(model.matrix(tempmod),"assign")
  parnames <- names(coef(tempmod))
  cluster <- if(is.null(random)) { random <- "case-wise" ; cluster <- 1:nrow(tempdata) } else { cluster <- data[,random] }

  ## CYCLE THROUGH EFFECTS
  out <- foreach(i=1:max(effects),.combine="rbind") %do% {
    wtest <- Wald_test(tempmod, constraints = constrain_zero(parnames[which(effects==i)]), vcov = "CR2", cluster=cluster, test = "HTZ")
    data.frame(Effect=attr(modinfo,"term.labels")[i],as.data.frame(wtest))
  }
  rm(tempdata,tempmod)
  out <- as.data.frame(out)
  out
  out$sig <- cut(out$p_val,breaks=c(0,0.001,0.01,0.05,0.1,1),labels=c("***","**","*",".",""))
  options(contrasts=c("contr.treatment","contr.treatment"))

  ## OUTPUT
  cat("Heteroscedastic Type II ANOVA breakdown:","\nCR2 estimator with Satterthwaite's correction for DDF","\n\n","Cluster variable:",random,"\n\n")
  out
}