# Function to perform an automated boxcox-transformation of the dependent variable 
# in various types and classes of models, including generalized and mixed effects models 

# Author: Sebastian G. Vetter-Lang (ORCID: 0000-0001-5374-5872)
# Version 1.1: Bug fix in line 73 and implementation of glmmTMB models
# Date: 2025-12-03


################################################################################
# arguments:

# object            A model of class lm, glm, gls, lme, lmerMod, glmerMod, glmmTMB, or negbin
# data              The data set used in the above model
# lambda            Power transformation parameter for boxcox transformation (see also ?bcPower),
#                   i.e. a range of lambda values to be tested to identify the optimal transformation
# suppress.plot     If TRUE the likelihood plot is suppressed 
# suppress.messages If TRUE the no message on the transformation is printed after the function finished 
# ...               Additional arguments passed to the model

################################################################################
# function: 
boxcox.ac <- function(object, data, lambda=seq(-2, 2, 1/10), suppress.plot=F, suppress.messages=F, ...){

  require(car)
  require(nortest)
  
  dat_ori <- data[,names(data)==as.character(formula(object))[2]]
  form <- as.formula(paste("trans ~", as.character(formula(object))[3]))
  likelihood=NULL
  x <- 0
  
  if (class(object)=="glm" || class(object)=="lm") {
    
    if(family(object)[1]=="binomial") { x <- 1 }
    for (l in lambda){
      if(all(dat_ori>x)) data$trans=bcPower(dat_ori,l)
      if(!all(dat_ori>x)) data$trans=yjPower(dat_ori,l)
      if(family(object)[1]!="gaussian") m=glm(form, family=family(object), data=data, ...)
      if(family(object)[1]=="gaussian") m=lm(form, data=data, ...)
      errtotal=m$res
      normal=ad.test(errtotal)
      likelihood=c(likelihood,normal$statistic)
    }
  } 
  
  if (class(object)[1]=="gls") {
    for (l in lambda){
      if(all(dat_ori>x)) data$trans=bcPower(dat_ori,l)
      if(!all(dat_ori>x)) data$trans=yjPower(dat_ori,l)
      m=gls(form, data=data, ...)
      errtotal=m$res
      normal=ad.test(errtotal)
      likelihood=c(likelihood,normal$statistic)
    }
  } 
  
  if (class(object)[1]=="lme") {
    require(nlme)
    for (l in lambda){
      if(all(dat_ori>x)) data$trans=bcPower(dat_ori,l)
      if(!all(dat_ori>x)) data$trans=yjPower(dat_ori,l)
      testError = tryCatch(lme(form, random=as.formula(as.character(object$call)[4]), data=data, ...), error=function(e){print(paste(e))})
      if(substr(testError[1],1,5)!="Error") {
        m=lme(form, random=as.formula(as.character(object$call)[4]), data=data, ...)
        errtotal=m$res[,1]
        normal=ad.test(errtotal)
        likelihood=c(likelihood,normal$statistic) } else {
          likelihood=c(likelihood,NA)
        }
    }
  }
  
  if (class(object)[1]=="glmerMod" || class(object)[1]=="lmerMod" || class(object)[1]=="lmerModLmerTest") {
    require(lme4)
    require(HLMdiag)
    if(length(grep("glmer", object@call)) >=1) {
      if(family(object)[1]=="binomial")  {
        x <- 1
      }
    }
    for (l in lambda){
      if(all(dat_ori>x)) data$trans=bcPower(dat_ori,l)
      if(!all(dat_ori>x)) data$trans=yjPower(dat_ori,l)
      if(length(grep("glmer", object@call)) >=1) {
        m=glmer(form, family=family(object), data=data)
      } else if (length(grep("lmer", object@call)) >= 1){
        m=lmer(form, data=data)
      }
      errtotal=hlm_resid(m,level=1)$.mar.resid
      normal=ad.test(errtotal)
      likelihood=c(likelihood,normal$statistic)
    }
  }
  
  if (class(object)[1]=="glmmTMB") {
    require(glmmTMB)
    if(family(object)[1]=="binomial")  {
        x <- 1
    }
    for (l in lambda){
      if(all(dat_ori>x)) data$trans=bcPower(dat_ori,l)
      if(!all(dat_ori>x)) data$trans=yjPower(dat_ori,l)
      m=glmmTMB(form, dispformula=object$modelInfo$allForm$dispformula, ziformula=object$modelInfo$allForm$ziformula, family=family(object), data=data)
      errtotal=residuals(m,re.form=~0)
      normal=ad.test(errtotal)
      likelihood=c(likelihood,normal$statistic)
    }
  }

  if (class(object)[1]=="negbin") {
    require(MASS)
    for (l in lambda){
      if(all(dat_ori>x)) data$trans=bcPower(dat_ori,l)
      if(!all(dat_ori>x)) data$trans=yjPower(dat_ori,l)
      m=glm.nb(form, data=data, ...)
      errtotal=m$res
      normal=ad.test(errtotal)
      likelihood=c(likelihood,normal$statistic)
    }
  }
  

  # plot likelihoods
  if(suppress.plot==F) {
  plot(lambda,likelihood, cex=0.75, pch=19)
  }
  lopt=lambda[which(likelihood==min(likelihood, na.rm=T))]
  
  if(lopt == min(lambda)) {
    warning("Optimal lambda NOT identified! Increase lambda range downwards! : boxcox.ac(..., lambda=seq(min, max, 0.1))")
  }
  
  if(lopt == max(lambda)) {
    warning("Optimal lambda NOT identified! Increase lambda range upwards! : boxcox.ac(..., lambda=seq(min, max, 0.1))")
  }
  
  if(lopt > min(lambda) & lopt < max(lambda)) {
    if(lopt>=-0.1 & lopt<=0.1) {
      if(suppress.messages==F) print(paste("log-transformation: opt. lambda = ", lopt, sep=""))
      lopt=0
    }
    
    # transformation of the response variable with lambda.max
    if(all(dat_ori>x)) {
      y_trans=bcPower(dat_ori,lopt)
      if(lopt==0) {
        if(suppress.messages==F) print("log transformation; backtransform with exp()")
        trans <- "log"
      }
      if(lopt!=0) {
        if(suppress.messages==F) print(paste("bcPower transformation; back-transform with ((y*", lopt, "+1)^(1/", lopt, "))", sep=""))
        trans <- "bcPower"
      }
    }  
    if(!all(dat_ori>x)) {
      y_trans=yjPower(dat_ori,lopt)
      if(suppress.messages==F) print(paste("yjPower transformation; back-transform with (((y*", lopt, "+1)^(1/", lopt, "))-1) OR -1*(((abs(y)*(2-", lopt, ")+1)^(1/(2-", lopt, ")))-1) for negative values", sep=""))
      trans <- "yiPower"
    }
    
    if(lopt>=0.9 & lopt<=1.1)  {
      if(suppress.messages==F) print(paste("no transformation nescessary: opt. lambda = ", lopt, sep=""))
      y_trans=dat_ori
      trans="no transfromation"
    }
    
    return(list("boxcox transformed data"=y_trans, "lambda"=lopt, "transformation"=trans))
  }
  
  
}