R/anova.aodml.R

In aods3: Analysis of Overdispersed Data using S3 Methods

anova.aodml <- function(object, ...) {

  #### copied from anova method for lmer (packages Matrix / lme4)

  mCall <- match.call(expand.dots = TRUE)
  dots <- list(...)
  modp <- sapply(dots, inherits, "aodml")
  mods <- c(list(object), dots[modp])
  names(mods) <- sapply(as.list(mCall)[c(FALSE, TRUE, modp)], as.character)

	####
	
	n <- length(mods)
  if(n < 2)
    stop("At least 2 valid models are needed.")
  
	dfr <- data.frame(
    df.resid = rep(NA, n), df.model = rep(NA, n), "2 x Log-lik" = rep(NA, n),
  	AIC = rep(NA, n), AICc = rep(NA, n), BIC = rep(NA, n),
  	Test = rep(NA, n), df = rep(NA, n), "LR stat." = rep(NA, n), 
  	"P(>LR stat.)" = rep(NA, n), check.names = FALSE
    )
	
	mod <- vector(mode = "character", length = length(mods))
	
	# utility function to remove white spaces at the beginning and at the end of character strings
  tr <- function(string) gsub("^[[:space:]]+|[[:space:]]+$", "", string)

  for(i in seq(n)){
  	#    fm <- get(nam[i])
    fm <- mods[[i]]
  	
    # model name and formula
  	# paste + tr are necessary when there are interaction terms in the formula
    fb <- paste(tr(deparse(fm$formula)), collapse = " ")
    ft <- deparse(fm$phi.formula)
    mod[i] <- paste(names(mods)[i], ": ", "Mu = ", fb, " ; Phi = ", ft, " ($phi = ", round(fm$phi, 4), ")", sep = "")
  	
    # anova table
    akic <- AIC(fm)
    df.resid <- df.residual(fm)
    df.model <- fm$df.model
    
    dfr$df.model[i] <- df.model
    dfr$df.resid[i] <- df.resid 
    dfr$AIC[i]  <- akic[, 3]
    dfr$AICc[i] <- akic[, 4]
    dfr$BIC[i]  <- -2 * fm$logL + log(df.model + df.resid) * df.model

    dfr[i, 3] <- 2 * fm$logL
  }
  
	dfr$Test <- c(NA, paste(names(mods)[-nrow(dfr)], names(mods)[-1], sep = " vs "))
  dfr$df <- c(NA, abs(diff(dfr$df.model)))
  dfr[, 9]  <- c(NA, abs(diff(dfr[, 3])))
  
	for(i in 2:n)
    dfr[i , 10] <- 1 - pchisq(dfr[i, 9], df = dfr$df[i])
  rownames(dfr) <- names(mods)
	
	structure(list(models = mod, anova.table = dfr), class = "anova.aodml")

}