R/reaster.R

In aster: Aster Models

reaster <- function(fixed, random, pred, fam, varvar, idvar, root,
    famlist = fam.default(), origin, data, effects, sigma, response)
    UseMethod("reaster")

reaster.default <- function(fixed, random, pred, fam, varvar, idvar, root,
    famlist = fam.default(), origin, data, effects, sigma, response)
{
    stopifnot(is.matrix(fixed))
    stopifnot(is.numeric(fixed))
    stopifnot(is.finite(fixed))

    if (! is.list(random))
        random <- list(random)
    for (irandom in seq(along = random)) {
        r <- random[[irandom]]
        if (! is.matrix(r))
            stop("random must be matrix or list of matrices")
        if (! is.numeric(r))
            stop("random must be numeric matrix or list of such")
        if (! all(is.finite(r)))
            stop("random must be all finite numeric matrix or list of such")
        if (nrow(r) != nrow(fixed))
            stop("fixed and random effect model matrices with different row dimensions")
    }

    stopifnot(is.vector(response))
    stopifnot(is.numeric(response))
    stopifnot(is.finite(response))
    if(length(response) != nrow(fixed))
        stop("length of response vector != row dimension of model matrices")

    stopifnot(is.vector(idvar))
    stopifnot(is.numeric(idvar))
    nind <- length(unique(idvar))

    stopifnot(is.factor(varvar))
    stopifnot(is.character(levels(varvar)))
    nnode <- nlevels(varvar)

    if (nind * nnode != length(response))
        stop("idvar or varvar wrong number of levels")

    varvarmat <- matrix(as.character(varvar), nind, nnode)
    idvarmat <- matrix(as.vector(idvar), nind, nnode)
    foo <- apply(varvarmat, 2, function(x) length(unique(x)))
    bar <- apply(idvarmat, 1, function(x) length(unique(x)))
    if (! (all(foo == 1) & all(bar == 1)))
        stop("idvar or varvar wrong structure")

    stopifnot(is.vector(pred))
    stopifnot(is.numeric(pred))
    stopifnot(pred %in% seq(0, nlevels(varvar)))

    stopifnot(is.vector(fam))
    stopifnot(is.numeric(fam))
    stopifnot(fam %in% seq(1, length(famlist)))
    famlist.blurfle <- sapply(famlist, as.character)
    famlist.problem <- grepl("negative.binomial", famlist.blurfle)
    if (any(famlist.problem[fam]))
        warning("negative binomial incompatible with random effects,",
            " see help(reaster)")

    stopifnot(is.vector(root))
    stopifnot(length(root) == length(response))
    stopifnot(is.numeric(root))
    stopifnot(is.finite(root))
    stopifnot(root > 0)

    if (! missing(origin)) {
        stopifnot(is.vector(origin))
        stopifnot(length(origin) == length(response))
        stopifnot(is.numeric(origin))
        stopifnot(is.finite(origin))
    }

    if (! missing(sigma)) {
        stopifnot(is.vector(sigma))
        stopifnot(is.numeric(sigma))
        stopifnot(is.finite(sigma))
        stopifnot(length(sigma) == length(as.list(random)))
    }

    ##### fit fixed effect aster model (even if effects supplied)
    ncoef <- ncol(fixed)
    modmat <- array(fixed, dim = c(nind, nnode, ncoef))
    dimnames(modmat) <- list(NULL, NULL, colnames(fixed))
    y <- matrix(response, nind, nnode)
    asterargs <- list(x = y, root = matrix(root, nind, nnode),
        pred = pred, fam = fam, modmat = modmat, famlist = famlist)
    if (! missing(origin))
        asterargs$origin <- matrix(origin, nind, nnode)
    aout <- try(do.call(aster, asterargs))
    if (inherits(aout, "try-error") || (! aout$converged)) {
        asterargs$method <- "nlm"
        aout <- try(do.call(aster, asterargs))
        if (inherits(aout, "try-error") || (! aout$converged))
            stop("cannot get convergence in fixed effects model")
    }

    ##### figure out what columns to drop from fixed effects matrix
    ikeep <- match(names(aout$coefficients), colnames(fixed))
    dropped <- colnames(fixed)[- ikeep]
    nfix.orig <- ncol(fixed)
    fixed <- fixed[ , ikeep]
    nfix <- ncol(fixed)
    idx.fixed <- seq(1, nfix)

    nrand <- sapply(random, ncol)
    blurfle <- names(random)
    if (is.null(blurfle)) {
        blurfle <- paste("R", seq(along = random), sep = "")
    } else {
        if (length(unique(blurfle)) != length(blurfle))
            stop("duplicate names for argument random")
    }
    splitter <- factor(rep(blurfle, times = nrand), levels = blurfle)
    idx.random <- nfix + seq(1, sum(nrand))

    if (! missing(effects)) {
        stopifnot(is.vector(effects))
        stopifnot(is.numeric(effects))
        stopifnot(is.finite(effects))
        if (length(dropped) != 0)
            if (length(effects) == nfix.orig + sum(nrand)) {
                warning("shortening effects to agree with columns dropped from fixed")
                effects.fixed <- effects[1:nfix.orig]
                effects.random <- effects[- (1:nfix.orig)]
                effects.fixed <- effects.fixed[ikeep]
                effects <- c(effects.fixed, effects.random)
            }
        if (length(effects) != ncol(fixed) + sum(sapply(random, ncol)))
            stop("length(effects) not sum of column dimensions of fixed and random effects model matrices")
    }

    if (missing(sigma) | missing(effects)) {
        sigma.start <- rep(1, length(nrand))
        eff.start <- c(aout$coefficients, rep(0, sum(nrand)))
        tout <- try(trust(objfun = penmlogl, eff.start,
            rinit = 1, rmax = 10,
            sigma = sigma.start, fixed = fixed, random = random,
            obj = aout, iterlim = 1000), silent = TRUE)
        if (inherits(tout, "try-error") || (! tout$converged))
            stop("cannot get convergence in penalized likelihood step")
        ### FIX ME ###
        ### need to try another optimization method when trust fails

        eff <- tout$argument * tout$scale
        eff.split <- split(eff[idx.random], splitter)
        sigma.start <- sapply(eff.split, function(x) sqrt(mean(x^2)))
    } else {
        sigma.start <- sigma
        eff.start <- effects
    }

    cache <- new.env(parent = emptyenv())
    oout <- suppressWarnings(try(optim(sigma.start, pickle, parm = eff.start,
        fixed = fixed, random = random, obj = aout, y = y, cache = cache),
        silent = TRUE))
    # this used to be
    #
    #     if (inherits(oout, "try-error") || oout$convergence != 0)
    #         stop("step 1 part 1 (optim Nelder-Mead with pickle) failed")
    #
    # but Marcus Warwell had an example where this optim with method
    # Nelder-Mead failed with out$convergence = 10 (indicates degeneracy
    # of the Nelder-Mead simplex, so says the documentation), and it is
    # very unclear what we are supposed to do with that.  Anyhoo, we
    # are only using this result as a starting point for more polishing,
    # so I hope this does not matter.

    if (inherits(oout, "try-error"))
        stop("step 1 part 1 (optim Nelder-Mead with pickle) failed")

    sigma.mle <- oout$par
    trout <- trust(objfun = penmlogl, cache$parm, rinit = 1, rmax = 10,
        sigma = sigma.mle, fixed = fixed, random = random, obj = aout, y = y)
    if (! trout$converged)
        stop("step 1 part 2 (trust with penmlogl) failed")
    parm.mle <- trout$argument
    zwz.mle <- makezwz(sigma.mle, parm = parm.mle,
        fixed = fixed, random = random, obj = aout, y = y)
    save.iter <- NULL
    value.mle <- NaN

    # now iterate pickle3 and trust
    foo <- function(alphaceesigma)
        pickle3(alphaceesigma, fixed = fixed, random = random,
            obj = aout, y = y, zwz = zwz.mle, deriv = 2)
    repeat {
        tout <- trust(foo, c(parm.mle, sigma.mle), rinit = 1, rmax = 10,
            iterlim = 1000)
        if (! tout$converged)
            stop("step 2 (trust with pickle3) failed")
        save.iter <- c(save.iter, tout$iterations)
        sigma.old <- as.vector(sigma.mle)
        sigma.mle <- tout$argument[nfix + sum(nrand) + seq(along = nrand)]
        parm.mle <- tout$argument[seq(1, nfix + sum(nrand))]
        zwz.mle <- makezwz(sigma.mle, parm.mle,
            fixed = fixed, random = random, obj = aout, y = y)
        value.mle <- tout$value
        if (isTRUE(all.equal(sigma.old, as.vector(sigma.mle)))) break
    }

    alpha.mle <- parm.mle[1:nfix]
    c.mle <- parm.mle[- (1:nfix)]

    # fix up negative sigma

    idx <- rep(seq(along = sigma.mle), times = nrand)
    for (k in seq(along = sigma.mle))
        if (sigma.mle[k] < 0) {
            sigma.mle[k] <- (- sigma.mle[k])
            eek <- (idx == k)
            c.mle[eek] <- (- c.mle[eek])
    }

    a.mle <- rep(sigma.mle, times = nrand)
    b.mle <- a.mle * c.mle

    if (! is.null(colnames(fixed)))
        names(alpha.mle) <- colnames(fixed)
    if (! is.null(colnames(random)))
        names(sigma.mle) <- names(random)
    if (all(! is.null(sapply(random, colnames)))) {
        names(b.mle) <- Reduce(c, Map(colnames, random))
        names(c.mle) <- names(b.mle)
    }

    # fix up estimated zeros
    # following code used to be in summary.reaster
    # but does not belong there
    if (missing(origin)) origin <- NULL
    nu.mle <- sigma.mle^2
    if (is.null(origin)) {
        iz <- is.zero(c(alpha.mle, b.mle, nu.mle), fixed, random,
            obj = aout, y, zwz = zwz.mle)
    } else {
        iz <- is.zero(c(alpha.mle, b.mle, nu.mle), fixed, random,
            obj = aout, y, origin = origin, zwz = zwz.mle)
    }
    nu.mle[iz] <- 0
    sigma.mle[iz] <- 0
    izbee <- rep(iz, times = nrand)
    b.mle[izbee] <- 0
    c.mle[izbee] <- 0

    result <- list(obj = aout, fixed = fixed, random = random,
        dropped = dropped, sigma = sigma.mle, nu = nu.mle,
        c = c.mle, b = b.mle, alpha = alpha.mle, zwz = zwz.mle,
        response = response, origin = origin,
        iterations = save.iter, counts = oout$counts,
        deviance = 2.0 * value.mle)
    class(result) <- c("reaster", "asterOrReaster")
    return(result)
}

reaster.formula <- function(fixed, random, pred, fam, varvar, idvar, root,
    famlist = fam.default(), origin, data, effects, sigma, response)
{
    stopifnot(inherits(fixed, "formula"))
    if (! is.list(random))
        random <- list(random)
    if (!  all(sapply(random, function(x) inherits(x, "formula"))))
        stop("random must be formula or list of formulas")
    stopifnot(missing(data) || is.data.frame(data))

    save.fixed <- fixed
    save.random <- random

    oldopt <- options(na.action = na.fail)
    on.exit(options(oldopt))

    ##### stuff copied from glm.R and not understood #####
    ##### see also http://developer.r-project.org/model-fitting-functions.txt

    call <- match.call()
    if(missing(data))
        data <- environment(fixed)
    mf <- match.call(expand.dots = FALSE)
    m <- match(c("fixed", "data", "varvar", "idvar", "root"), names(mf), 0)
    mf <- mf[c(1, m)]
    mf$drop.unused.levels <- TRUE
    mf[[1]] <- as.name("model.frame")
    fred <- as.list(mf)
    names(fred)[2] <- "formula"
    mf <- as.call(fred)
    mf <- eval.parent(mf)

    mt <- attr(mf, "terms")
    if (missing(response))
        response <- model.response(mf, "numeric")
    if (is.empty.model(mt)) {
        stop("empty model")
    } else {
        fixed <- model.matrix(mt, mf)
    }
    rownames(fixed) <- NULL

    varvar <- mf[["(varvar)"]]
    idvar <- mf[["(idvar)"]]
    root <- mf[["(root)"]]

    ##### end of stuff copied from glm.R and not understood #####

    for (irandom in seq(along = random)) {
        f <- random[[irandom]]
        foo <- lm(f, data = data, method = "model.frame")
        bar <- model.matrix(f, data = foo)
        rownames(bar) <- NULL
        random[[irandom]] <- bar
    }

    result <- NextMethod("reaster", response = response)

    result$formula <- list(fixed = save.fixed, random = save.random)
    result$call <- call
    class(result) <- c("reaster.formula", "reaster", "asterOrReaster")
    return(result)
}

summary.reaster <- function(object, standard.deviation = TRUE, ...)
{
    stopifnot(inherits(object, "reaster"))
    stopifnot(inherits(object$obj, "aster"))
    stopifnot(is.logical(standard.deviation))
    stopifnot(length(standard.deviation) == 1)

    alpha <- object$alpha
    sigma <- object$sigma
    nu <- object$nu
    cee <- object$c
    bee <- object$b
    zwz <- object$zwz
    fixed <- object$fixed
    random <- object$random
    if (is.matrix(random))
        random <- list(random)
    nfix <- ncol(fixed)
    nrand <- sapply(random, ncol)
    ########## REVIEWED DOWN TO HERE ##########

    obj <- object$obj
    y <- matrix(object$response, nrow = nrow(obj$x), ncol = ncol(obj$x))
    origin <- object$origin

    if (is.null(origin)) {
        iz <- is.zero(c(alpha, bee, nu), fixed, random, obj, y, zwz = zwz)
    } else {
        iz <- is.zero(c(alpha, bee, nu), fixed, random, obj, y,
            origin = origin, zwz = zwz)
    }
    nu[iz] <- 0
    sigma[iz] <- 0
    izbee <- rep(iz, times = nrand)
    bee[izbee] <- 0

    has.se <- c(rep(TRUE, length(alpha)), ! iz)
    if (all(iz)) {
        se.alpha <- sqrt(diag(solve(obj$fisher)))
        se.bee <- rep(NA_real_, length(bee))
        se.sigma <- rep(NA_real_, length(sigma))
        se.nu <- rep(NA_real_, length(nu))
        # want to return "subfish" created in other part
        # not defining subfish here is bug that causes crash
        # when this case (all variance components zero) occurs
        subfish <- obj$fisher
    } else {
        subrandom <- random[! iz]
        subnu <- nu[! iz]
        subbee <- bee[! izbee]
        subzwz <- zwz[! izbee, , drop = FALSE]
        subzwz <- subzwz[ , ! izbee, drop = FALSE]
        if (is.null(origin)) {
            qout <- quickle(c(alpha, subnu), subbee, fixed, subrandom,
                obj, y, zwz = subzwz, deriv = 2)
        } else {
            qout <- quickle(c(alpha, subnu), subbee, fixed, subrandom,
                obj, y, origin = origin, zwz = subzwz, deriv = 2)
        }
        subfish <- qout$hessian
    }
    eout <- eigen(subfish, symmetric = TRUE)
    goodfish <- min(eout$values) > 0
    if (goodfish) {
        # was buggy.  Have to use eigen to invert if use eigen to
        # check if invertible !!!!!!!!
        subfish.inv <- eout$vectors %*% diag(1 / eout$values) %*%
            t(eout$vectors)
    } else {
        warning(paste("estimated Fisher information matrix not positive",
            "definite, making all standard errors infinite"))
        subfish.inv <- matrix(NaN, nrow = nrow(subfish), ncol = ncol(subfish))
        diag(subfish.inv) <- Inf
    }
    se.subparm <- sqrt(diag(subfish.inv))
    se.parm <- rep(NA_real_, length(has.se))
    se.parm[has.se] <- se.subparm
    se.alpha <- se.parm[seq(along = alpha)]
    se.nu <- se.parm[- seq(along = alpha)]
    se.sigma <- se.nu / (2 * sigma)

    foo <- alpha
    foo <- cbind(foo, se.alpha)
    foo <- cbind(foo, foo[ , 1] / foo[ , 2])
    foo <- cbind(foo, 2 * pnorm(- abs(foo[ , 3])))
    rownames(foo) <- colnames(object$fixed)
    colnames(foo) <- c("Estimate", "Std. Error", "z value", "Pr(>|z|)")

    bar <- sigma
    bar <- cbind(bar, se.sigma)
    bar <- cbind(bar, bar[ , 1] / bar[ , 2])
    bar <- cbind(bar, pnorm(- abs(bar[ , 3])))
    blurfle <- names(object$random)
    if (is.null(blurfle))
        blurfle <- paste("sigma", seq(along = object$random), sep = "")
    rownames(bar) <- blurfle
    colnames(bar) <- c("Estimate", "Std. Error", "z value", "Pr(>|z|)/2")

    baz <- nu
    baz <- cbind(baz, se.nu)
    baz <- cbind(baz, baz[ , 1] / baz[ , 2])
    baz <- cbind(baz, pnorm(- abs(baz[ , 3])))
    blurfle <- names(object$random)
    if (is.null(blurfle))
        blurfle <- paste("nu", seq(along = object$random), sep = "")
    rownames(baz) <- blurfle
    colnames(baz) <- c("Estimate", "Std. Error", "z value", "Pr(>|z|)/2")

    return(structure(list(alpha = foo, sigma = bar, nu = baz,
        object = object, standard.deviation = standard.deviation,
        fisher = subfish, vcov = subfish.inv),
        class = "summary.reaster"))
}

print.summary.reaster <-
    function (x, digits = max(3, getOption("digits") - 3),
        signif.stars = getOption("show.signif.stars"), ...)
{
    if (! is.null(x$object$call)) {
        cat("\nCall:\n")
        cat(paste(deparse(x$object$call), sep="\n", collapse="\n"),
            "\n\n", sep="")
    }

    cat("\nFixed Effects:\n")
    printCoefmat(x$alpha, digits = digits,
        signif.stars = signif.stars, na.print = "NA", ...)

    if (x$standard.deviation) {
        cat("\nSquare Roots of Variance Components (P-values are one-tailed):\n")
        printCoefmat(x$sigma, digits = digits,
            signif.stars = signif.stars, na.print = "NA", ...)
    } else {
        cat("\nVariance Components (P-values are one-tailed):\n")
        printCoefmat(x$nu, digits = digits,
            signif.stars = signif.stars, na.print = "NA", ...)
    }

    if (! is.null(x$bee)) {
        cat("\nRandom Effects:\n")
        printCoefmat(x$bee, digits = digits,
            signif.stars = signif.stars, na.print = "NA", ...)
    }

    return(invisible(x))
}

vcov.reaster <- function(object, standard.deviation = TRUE,
    re.too = FALSE, complete = TRUE, ...)
{
    stopifnot(inherits(object, "reaster"))
    stopifnot(is.logical(standard.deviation))
    stopifnot(length(standard.deviation) == 1)
    stopifnot(is.logical(re.too))
    stopifnot(length(re.too) == 1)
    stopifnot(is.logical(complete))
    stopifnot(length(complete) == 1)

    # should not have to remake and re-evaluate objfun
    # but current code is not using objfun earlier
    objfun <- objfun.factory(object$fixed, object$random,
        object$response, object$obj$pred, object$obj$fam,
        as.vector(object$obj$root), object$zwz,
        standard.deviation = standard.deviation, deriv = 2)
    theta.hat <- if (standard.deviation) {
        c(object$alpha, object$c, object$sigma)
    } else {
        c(object$alpha, object$b, object$nu)
    }
    oout <- objfun(theta.hat)

    is.alpha <- seq_along(theta.hat) <= length(object$alpha)
    is.nu_or_sigma <-
        seq_along(theta.hat) > length(object$alpha) + length(object$b)
    is.bee_or_cee <- (! (is.alpha | is.nu_or_sigma))
    is.zero.nu_or_sigma <- object$nu == 0.0
    if (is.matrix(object$random))
       object$random <- list(object$random)
    nrand <- sapply(object$random, ncol)
    is.zero.bee_or_cee <- rep(is.zero.nu_or_sigma, times = nrand)
    is.zero <- c(rep(FALSE, sum(is.alpha)),
        is.zero.bee_or_cee, is.zero.nu_or_sigma)

    pW <- oout$hessian
    is.nonzero.bee <- is.bee_or_cee & (! is.zero)
    is.nonzero.psi <- is.alpha | ((! is.zero) & is.nu_or_sigma)
    pW_bee_bee <- pW[is.nonzero.bee, is.nonzero.bee]
    pW_psi_psi <- pW[is.nonzero.psi, is.nonzero.psi]
    pW_psi_bee <- pW[is.nonzero.psi, is.nonzero.bee]
    pW_bee_bee_inv <- .Call(C_pos_def_mat_inv, pW_bee_bee)
    if (inherits(pW_bee_bee_inv, "try-error"))
        stop("Hessian of pW not invertible; no standard errors")
    vcov_psi_psi_inv <- pW_psi_psi -
        pW_psi_bee %*% pW_bee_bee_inv %*% t(pW_psi_bee)
    vcov_psi_psi_inv <- (vcov_psi_psi_inv + t(vcov_psi_psi_inv)) / 2.0
    vcov_psi_psi <- .Call(C_pos_def_mat_inv, vcov_psi_psi_inv)
    if (inherits(vcov_psi_psi, "try-error"))
        stop("Hessian of qW not invertible; no standard errors")

    # Now vcov_psi_psi is the asymptotic variance-covariance matrix
    # of the elements of (alpha, nu) or (alpha, sigma) as the case
    # may be (depending on the value of argument standard.deviation)
    # that are not estimated to be zero (because those estimated to
    # be zero are not asymptotic normal)

    # names
    names.alpha <- names(object$alpha)
    if (is.null(names.alpha))
        names.alpha <- paste0("alpha", seq_along(object$alpha))
    names.bee <- names(object$b)
    if (is.null(names.bee))
        names.bee <- paste0(ifelse(standard.deviation, "c", "b"),
            seq_along(object$b))
    names.nu <- names(object$nu)
    if (is.null(names.nu))
        names.nu <- paste0(ifelse(standard.deviation, "sigma", "nu"),
            seq_along(object$nu))

    # still have to deal with arguments re.too and complete

    if ((! re.too) && ((! complete) || (! any(is.zero)))) {
        my.names.nu <- names.nu[! is.zero.nu_or_sigma]
        my.names <- c(names.alpha, my.names.nu)
        stopifnot(length(my.names) == nrow(vcov_psi_psi))
        stopifnot(length(my.names) == ncol(vcov_psi_psi))
        dimnames(vcov_psi_psi) <- list(my.names, my.names)
        attr(vcov_psi_psi, "is.alpha") <-
            is.alpha[(! is.bee_or_cee) & (! is.zero)]
        my.is <- is.nu_or_sigma[! is.bee_or_cee]
        my.zero <- is.zero[! is.bee_or_cee]
        if (standard.deviation) {
            attr(vcov_psi_psi, "is.sigma") <- my.is[! my.zero]
        } else {
            attr(vcov_psi_psi, "is.nu") <- my.is[! my.zero]
        }
        return(vcov_psi_psi)
    }

    if ((! re.too)) {
        # have some components to fill in with zeros
        my.names <- c(names.alpha, names.nu)
        result <- matrix(0.0, length(my.names), length(my.names))
        my.is.nonzero <- is.nonzero.psi[!is.bee_or_cee]
        dimnames(result) <- list(my.names, my.names)
        result[my.is.nonzero, my.is.nonzero] <- vcov_psi_psi
        attr(result, "is.alpha") <- is.alpha[! is.bee_or_cee]
        if (standard.deviation) {
            attr(result, "is.sigma") <- is.nu_or_sigma[! is.bee_or_cee]
        } else {
            attr(result, "is.nu") <- is.nu_or_sigma[! is.bee_or_cee]
        }
        return(result)
    }

    # now have only re.too == TRUE case (with either value of complete)

    # have to assemble matrix V that is equation (4.21) in the Aster Theory book
    # row dimension is total number of variables
    # column dimension is total number of components of alpha and nu
    #     that are not estimated to be equal to zero
    # all rows of vee for variables (alpha, b, and nu) estimated to be
    #     zero are all zeros
    b_star_psi <- (- pW_bee_bee_inv %*% t(pW_psi_bee))
    vee <- matrix(0.0, nrow = length(is.nonzero.psi),
        ncol = sum(is.nonzero.psi))
    vee[is.nonzero.psi, ] <- diag(sum(is.nonzero.psi))
    vee[is.nonzero.bee, ] <- b_star_psi

    if ((! complete) || (! any(is.zero))) {
        my.names <- c(names.alpha, names.bee, names.nu)
        my.names <- my.names[! is.zero]
        my.vee <- vee[! is.zero, ]
        result <- my.vee %*% vcov_psi_psi %*% t(my.vee)
        dimnames(result) <- list(my.names, my.names)
        attr(result, "is.alpha") <- is.alpha[! is.zero]
        if (standard.deviation) {
            attr(result, "is.cee") <- is.bee_or_cee[! is.zero]
            attr(result, "is.sigma") <- is.nu_or_sigma[! is.zero]
        } else {
            attr(result, "is.bee") <- is.bee_or_cee[! is.zero]
            attr(result, "is.nu") <- is.nu_or_sigma[! is.zero]
        }
        return(result)
    }

    # now left on re.too = TRUE and complete = TRUE

    my.vee <- diag(! is.zero) %*% vee
    result <- my.vee %*% vcov_psi_psi %*% t(my.vee)
    my.names <- c(names.alpha, names.bee, names.nu)
    dimnames(result) <- list(my.names, my.names)
    attr(result, "is.alpha") <- is.alpha
    if (standard.deviation) {
        attr(result, "is.cee") <- is.bee_or_cee
        attr(result, "is.sigma") <- is.nu_or_sigma
    } else {
        attr(result, "is.bee") <- is.bee_or_cee
        attr(result, "is.nu") <- is.nu_or_sigma
    }
    return(result)

}

anova.reaster <- function(object, ...)
{
    dotargs <- list(...)
    if (length(dotargs) == 0)
        stop("need at least two objects of class \"reaster\"")
    allargs <- c(list(object), dotargs)
    if (! all(sapply(allargs, function(x) inherits(x, "reaster"))))
        stop("some arguments not of class \"reaster\"")
    return(anova.reasterlist(allargs))
}

anova.reasterlist <- function(object, ...)
{
    stopifnot(is.list(object))
    stopifnot(length(object) >= 2)
    if (! all(sapply(object, function(x) inherits(x, "reaster"))))
        stop("some components not of class \"reaster\"")

    nmodels <- length(object)

    # attempt to check that models are nested
    # all models have same origin
    ok <- TRUE
    for (i in 2:nmodels) {
        o1 <- object[[i - 1]]$origin
        o2 <- object[[i]]$origin
        ok <- ok & identical(o1, o2)
    }
    if (! ok) warning("not same origin for all models, probably not nested")

    # all fixed effect model matrices have column labels
    ok <- TRUE
    for (i in 1:nmodels) {
        mf <- object[[i]]$fixed
        ok <- ok & (! is.null(colnames(mf)))
    }
    if (! ok) warning("no colnames for some model matrix for fixed effects, cannot check nesting")

    # all fixed effect model matrices column labels are nested
    ok <- TRUE
    for (i in 2:nmodels) {
        mf1 <- object[[i - 1]]$fixed
        mf2 <- object[[i]]$fixed
        ok <- ok & all(colnames(mf1) %in% colnames(mf2))
    }
    if (! ok) warning("colnames for model matrices for fixed effects not nested, models probably not nested")

    # all random effect model matrices have column labels
    ok <- TRUE
    for (i in 1:nmodels) {
        mr <- object[[i]]$random
        cr <- lapply(mr, colnames)
        ok <- ok & (! any(is.null(cr)))
    }
    if (! ok) warning("no colnames for some model matrix for random effects, cannot check nesting")

    # all random effect model matrices column labels are nested
    ok <- TRUE
    for (i in 2:nmodels) {
        mr1 <- object[[i - 1]]$random
        mr2 <- object[[i]]$random
        cr1 <- lapply(mr1, colnames)
        cr2 <- lapply(mr2, colnames)
        for (j1 in seq(along = cr1)) {
            ok.too <- FALSE
            for (j2 in seq(along = cr2)) {
                ok.too <- ok.too | identical(cr1[[j1]], cr2[[j2]])
            }
            ok <- ok & ok.too
        }
    }
    if (! ok) warning("colnames for model matrices for random effects not nested, models probably not nested")

    resdf.fix <- sapply(object, function(x) nrow(x$fixed))
    resdf.rand <- sapply(object, function(x) length(x$random))
    if (any(resdf.rand > 1)) stop("don't yet know how to compare models differing by two or more variance components")
    resdf <- resdf.fix + resdf.rand
    resdev <- sapply(object, function(x) x$deviance)
    table <- data.frame(resdf, resdev, c(NA, diff(resdf.fix)),
        c(NA, diff(resdf.rand)), c(NA, diff(resdev)))
    fred <- function(x) {
        foo <- x$fixed
        if (is.null(foo)) return("(no formulas)")
        result <- as.character(deparse(foo))
        for (i in seq(along = x$random)) {
            foo <- x$random[[i]]
            if (is.null(foo)) return("(no formulas)")
            more.result <- as.character(deparse(foo))
            result <- paste(result, more.result, sep = ", ")
        }
        return(result)
    }
    variables <- sapply(object, fred)
    dimnames(table) <- list(1:nmodels, c("Model Df", "Model Dev", "Df Fix",
        "Df Rand", "Deviance"))
    title <- "Analysis of Deviance Table\n"
    topnote <- paste("Model ", format(1:nmodels),": ",
        variables, collapse = "\n", sep = "")
    mixers <- table[ , "Df Rand"] == 1
    pval <- double(nrow(table))
    pval[! mixers] <- pchisq(table[ , "Deviance"],
        table[ , "Df Fix"], lower.tail = FALSE)
    pval[mixers] <- (0.5) * pchisq(table[ , "Deviance"],
        table[ , "Df Fix"], lower.tail = FALSE) +
        pchisq(table[ , "Deviance"], table[ , "Df Fix"] + 1,
        lower.tail = FALSE)
    table <- cbind(table, "P-value" = pval)
    structure(table, heading = c(title, topnote),
        class = c("anova", "data.frame"))
}