anova.asterOrReaster <- function(object, ..., tolerance = .Machine$double.eps ^ 0.75) { dotargs <- list(...) if (length(dotargs) == 0) stop("need at least two models to compare") allargs <- c(list(object), dotargs) if (! all(sapply(allargs, function(x) inherits(x, "asterOrReaster")))) stop("some arguments not of class \"asterOrReaster\"") return(anovaAsterOrReasterList(allargs, tolerance)) } anovaAsterOrReasterList <- function(objectlist, tolerance = .Machine$double.eps ^ 0.75) { stopifnot(is.list(objectlist)) if (length(objectlist) < 2) stop("must compare two or more models") if (! all(sapply(objectlist, function(x) inherits(x, "asterOrReaster")))) stop("some components not of class \"asterOrReaster\"") nmodels <- length(objectlist) stopifnot(is.numeric(tolerance)) stopifnot(length(tolerance) == 1) stopifnot(tolerance > 0) # attempt to check that models are nested # FIXME: this is bad, but don't check offset vectors for now # assume default offset vector was used get.fixed <- function(x) { if(inherits(x, "reaster")) { return(x$fixed) } if(inherits(x, "aster")) { bar <- x$modmat nparm <- dim(bar)[3] baz <- matrix(as.vector(bar), ncol = nparm) colnames(baz) <- dimnames(bar)[[3]] return(baz) } stop("model object not class \"aster\" or \"reaster\"") } fixed <- lapply(objectlist, get.fixed) get.random <- function(x) { if(inherits(x, "reaster")) return(x$random) if(inherits(x, "aster")) return(list()) stop("model object not class \"aster\" or \"reaster\"") } random <- lapply(objectlist, get.random) nnode.fixed <- sapply(fixed, nrow) nnode.random <- lapply(fixed, function(x) sapply(x, nrow)) nnode <- c(nnode.fixed, unlist(nnode.random)) if (length(unique(nnode)) != 1) stop("model matrices do not all have the same row dimension") nnode <- unique(nnode) df.fixed <- sapply(fixed, ncol) df.random <- sapply(random, length) if (any(diff(df.fixed) < 0)) stop("degrees of freedom for fixed effects not nondecreasing") if (any(diff(df.random) < 0)) stop("degrees of freedom for random effects not nondecreasing") # check model matrix for fixed effects nesting qrList <- lapply(fixed, qr) for (i in 2:nmodels) { modmat1 <- fixed[[i - 1]] qr2 <- qrList[[i]] resid1 <- qr.resid(qr2, modmat1) norm.modmat1 <- apply(modmat1, 2, function(x) sqrt(sum(x^2))) norm.resid1 <- apply(resid1, 2, function(x) sqrt(sum(x^2))) if (any(norm.resid1 / norm.modmat1 > tolerance)) stop("model matrices for fixed effects not nested\nmodel ", i - 1, " and model ", i) } # check model matrix for random effects nesting # here we hope simple check on column labels is o. k. collab.random <- lapply(random, function(x) lapply(x, colnames)) foo.random <- lapply(collab.random, function(x) sapply(x, is.null)) if (any(unlist(foo.random))) { warning("some model matrices for random effects have no colnames\nnot checking for nested models") } else { for (i in 2:nmodels) { c1 <- collab.random[[i - 1]] c2 <- collab.random[[i]] p1 <- sapply(c1, paste, collapse = " ") p2 <- sapply(c2, paste, collapse = " ") if (! all(p1 %in% p2)) stop("model matrices for random effects not nested\nmodel ", i - 1, " and model ", i) } } resdf.fixed <- diff(df.fixed) resdf.random <- diff(df.random) if (any(resdf.random > 1)) stop("don't know how to compare models differing by two or more variance components") resdev <- sapply(objectlist, function(x) x$deviance) get.formulae <- function(x) { if(inherits(x, "aster")) { qux <- x$formula if (is.null(qux)) return("(no formulas)") return(as.character(deparse(qux))) } if(inherits(x, "reaster")) { qux <- x$formula if (is.null(qux)) return("(no formulas)") qux <- unlist(qux) quacks <- sapply(qux, function(x) as.character(deparse(x))) return(paste(quacks, collapse = ", ")) } stop("model object not class \"aster\" or \"reaster\"") } formulae <- sapply(objectlist, get.formulae) if (all(df.random == 0)) { # all aster models table <- data.frame(df.fixed, - resdev, c(NA, resdf.fixed), c(NA, diff(- resdev))) dimnames(table) <- list(1:nmodels, c("Model Df", "Model Dev", "Df", "Deviance")) title <- "Analysis of Deviance Table\n" topnote <- paste("Model ", format(1:nmodels), ": ", formulae, collapse = "\n", sep = "") table <- cbind(table, "P(>|Chi|)" = pchisq(table[ , "Deviance"], table[ , "Df"], lower.tail = FALSE)) table <- structure(table, heading = c(title, topnote), class = c("anova", "data.frame")) } else { # some, perhaps all reaster models table <- data.frame(df.fixed, df.random, - resdev, c(NA, resdf.fixed), c(NA, resdf.random), c(NA, diff(- resdev))) dimnames(table) <- list(1:nmodels, c("Mod Df Fix", "Mod Df Rand", "Mod Dev", "Df Fix", "Df Rand", "Deviance")) title <- "Analysis of Deviance Table\n" topnote <- paste("Model ", format(1:nmodels), ": ", formulae, collapse = "\n", sep = "") d <- table[ , "Deviance"] nf <- table[ , "Df Fix"] nr <- table[ , "Df Rand"] p <- (1 / 2) * pchisq(d, nf, lower.tail = FALSE) + (1 / 2) * pchisq(d, nf + nr, lower.tail = FALSE) table <- cbind(table, "P-value" = p) table <- structure(table, heading = c(title, topnote), class = c("anova", "data.frame")) } return(table) }