#  File src/library/stats/R/aov.R
#  Part of the R package, https://www.R-project.org
#
#  Copyright (C) 1995-2020 The R Core Team
#  Copyright (C) 1998 B. D. Ripley
#
#  This program is free software; you can redistribute it and/or modify
#  it under the terms of the GNU General Public License as published by
#  the Free Software Foundation; either version 2 of the License, or
#  (at your option) any later version.
#
#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public License for more details.
#
#  A copy of the GNU General Public License is available at
#  https://www.R-project.org/Licenses/

aov <- function(formula, data = NULL, projections = FALSE, qr = TRUE,
                contrasts = NULL, ...)
{
    Terms <- if(missing(data))
		  terms(formula, "Error")
	     else terms(formula, "Error", data = data)
    indError <- attr(Terms, "specials")$Error
    ## NB: this is only used for n > 1, so singular form makes no sense
    ## in English.  But some languages have multiple plurals.
    if(length(indError) > 1L)
        stop(sprintf(ngettext(length(indError),
                              "there are %d Error terms: only 1 is allowed",
                              "there are %d Error terms: only 1 is allowed"),
                     length(indError)), domain = NA)
    lmcall <- Call <- match.call()
    ## need stats:: for non-standard evaluation
    lmcall[[1L]] <- quote(stats::lm)
    lmcall$singular.ok <- TRUE
    if(projections) qr <- lmcall$qr <- TRUE
    lmcall$projections <- NULL
    if(is.null(indError)) {
	## no Error term
	fit <- eval(lmcall, parent.frame())
	fit$call <- Call # fixup 'lmcall'
	structure(fit,
		  class = c(if(inherits(fit, "mlm")) "maov",
			    "aov", oldClass(fit)),
		  projections = if(projections) proj(fit))
    } else {
	if(pmatch("weights", names(Call), 0L))
            stop("weights are not supported in a multistratum aov() fit")
        deparseb <- function(expr) deparse1(expr, backtick = TRUE)
        ##  Helmert contrasts can be helpful: do we want to force them?
        ##  this version does for the Error model.
        opcons <- options("contrasts")
        options(contrasts = c("contr.helmert", "contr.poly"))
        on.exit(options(opcons))
        allTerms <- Terms
        errorterm <- attr(Terms, "variables")[[1L + indError]]
        intercept <- attr(Terms, "intercept")
        ecall <- lmcall
        ecall$formula <-
            as.formula(paste(deparseb(formula  [[2L]]), "~",
                             deparseb(errorterm[[2L]]), if(!intercept) "- 1"),
                       env = environment(formula))
        ecall$method <- "qr"
        ecall$qr <- TRUE
        ecall$contrasts <- NULL
        er.fit <- eval(ecall, parent.frame())
        options(opcons)
        nmstrata <- attr(terms(er.fit), "term.labels")
        ## remove backticks from simple labels for strata (only)
        nmstrata <- sub("^`(.*)`$", "\\1", nmstrata)
        nmstrata <- c("(Intercept)", nmstrata)
        qr.e <- er.fit$qr
        rank.e <- er.fit$rank
        if(rank.e < NROW(er.fit$coefficients))
            warning("Error() model is singular")
        qty <- er.fit$residuals
        maov <- is.matrix(qty)
        asgn.e <- er.fit$assign[qr.e$pivot[1L:rank.e]]
        ## we want this to label the rows of qtx, not cols of x.
        maxasgn <- length(nmstrata) - 1L
        nobs <- NROW(qty)
	len <- if(nobs > rank.e) {
	    asgn.e[(rank.e+1L):nobs] <- maxasgn + 1L
	    nmstrata <- c(nmstrata, "Within")
	    maxasgn + 2L
	} else maxasgn + 1L
        result <- setNames(vector("list", len), nmstrata)
        lmcall$formula <- form <-
            update(formula, paste(". ~ .-", deparseb(errorterm)))
        Terms <- terms(form)
        lmcall$method <- "model.frame"
        mf <- eval(lmcall, parent.frame())
        xlev <- .getXlevels(Terms, mf)
        resp <- model.response(mf)
        qtx <- model.matrix(Terms, mf, contrasts)
        cons <- attr(qtx, "contrasts")
        dnx <- colnames(qtx)
        asgn.t <- attr(qtx, "assign")
        if(length(wts <- model.weights(mf))) {
            wts <- sqrt(wts)
            resp <- resp * wts
            qtx <- qtx * wts
        }
        qty <- as.matrix(qr.qty(qr.e, resp))
        if((nc <- ncol(qty)) > 1L) {
            dny <- colnames(resp) %||% paste0("Y", 1L:nc)
            dimnames(qty) <- list(seq(nrow(qty)), dny)
        } else dimnames(qty) <- list(seq(nrow(qty)), NULL)
        qtx <- qr.qty(qr.e, qtx)
        dimnames(qtx) <- list(seq(nrow(qtx)) , dnx)
        for(i in seq_along(nmstrata)) {
            select <- asgn.e == (i-1L)
            ni <- sum(select)
            if(!ni) next
            ## helpful to drop constant columns.
            xi <- qtx[select, , drop = FALSE]
            cols <- colSums(xi^2) > 1e-5
            if(any(cols)) {
                xi <- xi[, cols, drop = FALSE]
                attr(xi, "assign") <- asgn.t[cols]
                fiti <- lm.fit(xi, qty[select,,drop=FALSE])
                fiti$terms <- Terms
            } else {
                y <- qty[select,,drop=FALSE]
                fiti <- list(coefficients = numeric(), residuals = y,
                             fitted.values = 0 * y, weights = wts, rank = 0L,
                             df.residual = NROW(y))
            }
            if(projections) fiti$projections <- proj(fiti)
            class(fiti) <- c(if(maov) "maov", "aov", oldClass(er.fit))
            result[[i]] <- fiti
        }
        structure(class = c("aovlist", "listof"),
                  ## drop empty strata:
                  result[!vapply(result, is.null, NA)],
                  error.qr = if(qr) qr.e,
                  call = Call,
                  weights = if(length(wts)) wts,
                  terms = allTerms,
                  contrasts = cons,
                  xlevels = xlev)
    }
}

print.aov <-
function(x, intercept = FALSE, tol = sqrt(.Machine$double.eps), ...)
{
    if(!is.null(cl <- x$call)) {
        cat("Call:\n   ")
        dput(cl, control = NULL)
    }
    qrx <- if(x$rank) qr(x)
    asgn <- x$assign[qrx$pivot[1L:x$rank]]
    effects <- x$effects
    if(!is.null(effects))
        effects <- as.matrix(effects)[seq_along(asgn),,drop=FALSE]
    rdf <- x$df.residual
    resid <- as.matrix(x$residuals)
    wt <- x$weights
    if(!is.null(wt)) resid <- resid * sqrt(wt)
    RSS <- colSums(resid^2)
    uasgn <- unique(asgn)
    nmeffect <- c("(Intercept)", attr(x$terms, "term.labels"))[1+uasgn]
    nterms <- length(uasgn)
    nresp <- NCOL(effects)
    df <- numeric(nterms)
    ss <- matrix(NA, nterms, nresp)
    if(nterms) {
        for(i in seq(nterms)) {
            ai <- asgn == uasgn[i]
            df[i] <- sum(ai)
            ef <- effects[ai,, drop=FALSE]
            ss[i, ] <- if(sum(ai) > 1) colSums(ef^2) else ef^2
        }
        keep <- df > 0L
        if(!intercept && uasgn[1L] == 0) keep[1L] <- FALSE
        nmeffect <- nmeffect[keep]
        df <- df[keep]
        ss <- ss[keep, , drop = FALSE]
        nterms <- length(df)
    }
    cat("\nTerms:\n")
    if(nterms == 0L) {
        ## empty model
        if(rdf > 0L) {
            ss <- RSS
            ssp <- sapply(ss, format)
            if(!is.matrix(ssp)) ssp <- t(ssp)
            tmp <- as.matrix(c(ssp, format(rdf)))
            if(length(ss) > 1L) {
                rn <- colnames(x$fitted.values)
                if(is.null(rn)) rn <- paste("resp", seq_along(ss))
            } else rn <- "Sum of Squares"
            dimnames(tmp) <- list(c(rn, "Deg. of Freedom"), "Residuals")
            print(tmp, quote = FALSE, right = TRUE)
            cat("\n")
            rs <- sqrt(RSS/rdf)
            cat(if(length(rs) > 1L) "Residual standard errors:"
                 else "Residual standard error:", sapply(rs, format))
            cat("\n")
        } else
        print(matrix(0, 2L, 1L, dimnames =
                     list(c("Sum of Squares", "Deg. of Freedom"), "<empty>")))
    } else {
        if(rdf > 0L) {
            nterms <- nterms + 1L
            df <- c(df, rdf)
            ss <- rbind(ss, RSS)
            nmeffect <- c(nmeffect, "Residuals")
        }
        ssp <- apply(zapsmall(ss), 2L, format)
        tmp <- t(cbind(ssp, format(df)))
        if(ncol(effects) > 1L) {
            rn <- colnames(x$coefficients)
            if(is.null(rn)) rn <- paste("resp", seq(ncol(effects)))
        } else rn <- "Sum of Squares"
        dimnames(tmp) <- list(c(rn, "Deg. of Freedom"), nmeffect)
        print(tmp, quote = FALSE, right = TRUE)
        rank <- x$rank
#        int <- attr(x$terms, "intercept")
#        nobs <- NROW(x$residuals) - !(is.null(int) || int == 0)
        cat("\n")
        if(rdf > 0L) {
            rs <- sqrt(RSS/rdf)
            cat(if(length(rs) > 1L) "Residual standard errors:"
                else "Residual standard error:", sapply(rs, format))
            cat("\n")
        }
        coef <- as.matrix(x$coefficients)[, 1L]
        R <- qrx$qr
        R <- R[1L:min(dim(R)), , drop=FALSE]
        R[lower.tri(R)] <- 0
        if(rank < (nc <- length(coef))) {
            cat(paste(nc - rank, "out of", nc, "effects not estimable\n"))
            R <- R[, 1L:rank, drop = FALSE]
        }
        d2 <- sum(abs(diag(R)))
        diag(R) <- 0
        if(sum(abs(R))/d2 > tol)
            cat("Estimated effects may be unbalanced\n")
        else cat("Estimated effects are balanced\n")
        if(nzchar(mess <- naprint(x$na.action))) cat(mess, "\n", sep = "")
    }
    invisible(x)
}

summary.aov <- function(object, intercept = FALSE, split,
                        expand.split = TRUE, keep.zero.df = TRUE, ...)
{
    splitInteractions <- function(split, factors, names, asgn, df.names)
    {
        ns <- names(split)
        for(i in unique(asgn)) {
            if(i == 0 || names[i+1L] %in% ns) next
            f <- rownames(factors)[factors[, i] > 0]
            sp <- f %in% ns
            if(any(sp)) {              # some marginal terms are split
                if(sum(sp) > 1L) {
                    old <- split[ f[sp] ]
		    nn <- setNames(nm = f[sp])
                    marg <- lapply(nn, function(x)
                                   df.names[asgn == (match(x, names) - 1L)])
                    term.coefs <- strsplit(df.names[asgn == i], ":", fixed=TRUE)
                    ttc <- sapply(term.coefs, function(x) x[sp])
                    rownames(ttc) <- nn
		    splitnames <-
			setNames(nm = apply(expand.grid(lapply(old, names)), 1L,
				 function(x) paste(x, collapse=".")))
                    tmp <- sapply(nn, function(i)
                                  names(old[[i]])[match(ttc[i, ], marg[[i]])] )
                    tmp <- apply(tmp, 1L, function(x) paste(x, collapse="."))
                    new <- lapply(splitnames, function(x) match(x, tmp))
                    split[[ names[i+1L] ]] <-
                        new[sapply(new, function(x) length(x) > 0L)]
                } else {
                    old <- split[[ f[sp] ]]
                    marg.coefs <- df.names[asgn == (match(f[sp], names) - 1L)]
                    term.coefs <- strsplit(df.names[asgn == i], ":", fixed=TRUE)
                    ttc <- sapply(term.coefs, function(x) x[sp])
                    new <- lapply(old, function(x)
                                  seq_along(ttc)[ttc %in% marg.coefs[x]])
                    split[[ names[i+1L] ]] <- new
                }
            }
        }
        split
    }

    asgn <- object$assign[object$qr$pivot[1L:object$rank]]
    uasgn <- unique(asgn)
    nterms <- length(uasgn)
    effects <- object$effects
    if(!is.null(effects))
        effects <- as.matrix(effects)[seq_along(asgn),,drop=FALSE]
    rdf <- object$df.residual
    nmeffect <- c("(Intercept)", attr(object$terms, "term.labels"))
    coef <- as.matrix(object$coefficients)
    resid <- as.matrix(object$residuals)
    wt <- object$weights
    if(!is.null(wt)) resid <- resid * sqrt(wt)
    nresp <- NCOL(resid)
    ans <- vector("list", nresp)
    if(nresp > 1) {
        names(ans) <- character(nresp)
        for (y in 1L:nresp) {
            cn <- colnames(resid)[y]
            if(is.null(cn) || cn == "") cn <- y
            names(ans)[y] <- paste(" Response", cn)
        }
    }

    if(!is.null(effects) && !missing(split)) {
        ns <- names(split)
        if(!is.null(Terms <- object$terms)) {
            if(!is.list(split))
                stop("the 'split' argument must be a list")
            if(!all(ns %in% nmeffect)) {
                na <- sum(!ns %in% nmeffect)
                stop(sprintf(ngettext(na,
                                      "unknown name %s in the 'split' list",
                                      "unknown names %s in the 'split' list"),
                             paste(sQuote(ns[na]), collapse = ", ")),
                     domain = NA)
            }
        }
        if(expand.split) {
            df.names <- names(coef(object))
            split <- splitInteractions(split, attr(Terms, "factors"),
                                       nmeffect, asgn, df.names)
            ns <- names(split)
        }
    }

    for (y in 1L:nresp) {
        if(is.null(effects)) {
            nterms <- 0L
            df <- ss <- ms <- numeric()
            nmrows <- character()
        } else {
            df <- ss <- numeric()
            nmrows <- character()
            for(i in seq(nterms)) {
                ai <- (asgn == uasgn[i])
                df <- c(df, sum(ai))
                ss <- c(ss, sum(effects[ai, y]^2))
                nmi <- nmeffect[1 + uasgn[i]]
                nmrows <- c(nmrows, nmi)
                if(!missing(split) && !is.na(int <- match(nmi, ns))) {
		    df <- c(df, lengths(split[[int]]))
                    if(is.null(nms <- names(split[[int]])))
                        nms <- paste0("C", seq_along(split[[int]]))
                    ss <- c(ss, unlist(lapply(split[[int]],
                                              function(i, e)
                                              sum(e[i]^2), effects[ai, y])))
                    nmrows <- c(nmrows, paste0("  ", nmi, ": ", nms))
                }
            }
        }
        if(rdf > 0L) {
            df <- c(df, rdf)
            ss <- c(ss, sum(resid[, y]^2))
            nmrows <- c(nmrows,  "Residuals")
        }
        nt <- length(df)
        ms <- ifelse(df > 0L, ss/df, NA)
        x <- list(Df = df, "Sum Sq" = ss, "Mean Sq" = ms)
        if(rdf > 0L) {
            TT <- ms/ms[nt]
            TP <- pf(TT, df, rdf, lower.tail = FALSE)
            TT[nt] <- TP[nt] <- NA
            x$"F value" <- TT
            x$"Pr(>F)" <- TP
            ## 'nterms' ~= 'Residuals' have no P-value
        }
        class(x) <- c("anova", "data.frame")
        attr(x, "row.names") <- format(nmrows)
        if(!keep.zero.df) x <- x[df > 0L, ]
        pm <- pmatch("(Intercept)", row.names(x), 0L)
        if(!intercept && pm > 0L) x <- x[-pm ,]
        ans[[y]] <- x
    }
    class(ans) <- c("summary.aov", "listof")
    attr(ans, "na.action") <- object$na.action
    ans
}

print.summary.aov <-
    function(x, digits = max(3L, getOption("digits") - 3L),
             symbolic.cor = FALSE,
             signif.stars = getOption("show.signif.stars"), ...)
{
    if (length(x) == 1L)
        print(x[[1L]], digits = digits, symbolic.cor = symbolic.cor,
              signif.stars = signif.stars)
    else NextMethod()
    if(nzchar(mess <- naprint(attr(x, "na.action")))) cat(mess, "\n", sep = "")
    invisible(x)
}

## coef.aov() := coef.default _but_ with different default for 'complete'
## --> in ./models.R

# For maov objects, the coefficients are a matrix and
# NAs can't sensibly be removed (PR#16380)
coef.maov <- function(object, ...)
    object$coefficients

alias <- function(object, ...) UseMethod("alias")

alias.formula <- function(object, data, ...)
{
    lm.obj <- if(missing(data)) aov(object) else aov(object, data)
    alias(lm.obj, ...)
}

alias.lm <- function(object, complete = TRUE, partial = FALSE,
                     partial.pattern = FALSE, ...)
{
    CompPatt <- function(x, ...) {
        x[abs(x) < 1e-6] <- 0
        MASS::fractions(x)
    }
    PartPatt <- function(x) {
        z <- zapsmall(x) != 0
        if(any(z)) {
            xx <- abs(signif(x[z], 2))
            ll <- length(unique(xx))
            if(ll > 10L) xx <- cut(xx, 9L) else if(ll == 1L) x[] <- 1
            x[z] <- paste0(ifelse(x[z] > 0, " ", "-"), xx)
        }
        x[!z] <- ""
        collabs <- colnames(x)
        collabs <- if(length(collabs))
                       abbreviate(sub(".", "", collabs, fixed=TRUE), 3L)
                   else 1L:ncol(x)
        colnames(x) <- collabs
        class(x) <- "mtable"
        x
    }
    Model <- object$terms
    attributes(Model) <- NULL
    value <- list(Model = Model)
    R <- qr(object)$qr
    R <- R[1L:min(dim(R)), , drop=FALSE]
    R[lower.tri(R)] <- 0
    d <- dim(R)
    rank <- object$rank
    p <- d[2L]
    if(complete) {                      # full rank, no aliasing
        value$Complete <-
            if(is.null(p) || rank == p) NULL else {
                p1 <- 1L:rank
                X <- R[p1, p1]
                Y <-  R[p1, -p1, drop = FALSE]
                beta12 <- as.matrix(qr.coef(qr(X), Y))
                # dimnames(beta12) <- list(dn[p1], dn[ -p1])
                CompPatt(t(beta12))
            }
    }
    if(partial) {
        ## We only want one aspect of the summary, which we know to be reliable
        tmp <- suppressWarnings(summary.lm(object)$cov.unscaled)
        ses <- sqrt(diag(tmp))
        beta11 <- tmp /outer(ses, ses)
        beta11[row(beta11) >= col(beta11)] <- 0
        beta11[abs(beta11) < 1e-6] <- 0
        if(all(beta11 == 0)) beta11 <- NULL
        else if(partial.pattern) beta11 <- PartPatt(beta11)
        value$Partial <- beta11
    }
    class(value) <- "listof"
    value
}

print.aovlist <- function(x, ...)
{
    cl <- attr(x, "call")
    if(!is.null(cl)) {
        cat("\nCall:\n")
        dput(cl)
    }
    if(!is.null(attr(x, "weights")))
        cat("Note: The results below are on the weighted scale\n")
    nx <- names(x)
    if(nx[1L] == "(Intercept)") {
        mn <- x[[1L]]$coefficients
        if(is.matrix(mn)) {
            cat("\nGrand Means:\n")
            print(format(mn[1,]), quote = FALSE)
        } else cat("\nGrand Mean: ", format(mn[1L]), "\n", sep = "")
        nx <- nx[-1L]
    }
    for(ii in seq_along(nx)) {
        i <- nx[ii]
        cat("\nStratum ", ii, ": ", i, "\n", sep = "")
        xi <- x[[i]]
        print(xi, ...)
    }
    invisible(x)
}

summary.aovlist <- function(object, ...)
{
    if(!is.null(attr(object, "weights")))
        cat("Note: The results below are on the weighted scale\n")
    dots <- list(...)
    strata <- names(object)
    if(strata[1L] == "(Intercept)") {
        strata <- strata[-1L]
        object <- object[-1L]
    }
    x <- setNames(vector(length = length(strata), mode = "list"),
		  paste("Error:", strata))
    for(i in seq_along(strata))
        x[[i]] <- do.call("summary", c(list(object = object[[i]]), dots))
    class(x) <- "summary.aovlist"
    x
}

print.summary.aovlist <- function(x, ...)
{
    nn <- names(x)
    for (i in nn) {
        cat("\n", i, "\n", sep = "")
        print(x[[i]], ...)
    }
    invisible(x)
}

coef.listof <- function(object, ...)
{
    val <- setNames(vector("list", length(object)), names(object))
    for(i in seq_along(object)) val[[i]] <- coef(object[[i]])
    class(val) <- "listof"
    val
}

se.contrast <- function(object, ...) UseMethod("se.contrast")

se.contrast.aov <-
    function(object, contrast.obj,
             coef = contr.helmert(ncol(contrast))[, 1L],
             data = NULL, ...)
{
    contrast.weight.aov <- function(object, contrast)
    {
        qro <- qr(object)
        asgn <- object$assign[qro$pivot[1L:object$rank]]
        uasgn <- unique(asgn)
        nterms <- length(uasgn)
        nmeffect <- c("(Intercept)",
                      attr(object$terms, "term.labels"))[1L + uasgn]
        effects <- as.matrix(qr.qty(qro, contrast))
        res <- matrix(0, nrow = nterms, ncol = ncol(effects),
                      dimnames = list(nmeffect, colnames(contrast)))
        for(i in seq(nterms)) {
            select <- (asgn == uasgn[i])
            res[i,] <- colSums(effects[seq_along(asgn)[select], , drop = FALSE]^2)
        }
        res
    }
    if(is.null(data)) contrast.obj <- eval(contrast.obj)
    else contrast.obj <- eval(substitute(contrast.obj), data, parent.frame())
    if(!is.matrix(contrast.obj)) { # so a list
        if(!missing(coef)) {
            if(sum(coef) != 0)
                stop("'coef' must define a contrast, i.e., sum to 0")
            if(length(coef) != length(contrast.obj))
                stop("'coef' must have same length as 'contrast.obj'")
        }
        contrast <-
            sapply(contrast.obj, function(x)
               {
                   if(!is.logical(x))
                       stop(gettextf("each element of '%s' must be logical",
                                     substitute(contrasts.list)),
                            domain = NA)
                   x/sum(x)
               })
        if(!length(contrast) || all(is.na(contrast)))
            stop("the contrast defined is empty (has no TRUE elements)")
        contrast <- contrast %*% coef
    } else {
        contrast <- contrast.obj
        if(any(abs(colSums(contrast)) > 1e-8))
            stop("columns of 'contrast.obj' must define a contrast (sum to zero)")
        if(!length(colnames(contrast)))
            colnames(contrast) <- paste("Contrast", seq(ncol(contrast)))
    }
    weights <- contrast.weight.aov(object, contrast)
    rdf <- object$df.residual
    if(rdf == 0L) stop("no degrees of freedom for residuals")
    resid <- as.matrix(object$residuals)
    wt <- object$weights
    if(!is.null(wt)) resid <- resid * sqrt(wt)
    rse <- sum(resid^2)/rdf
    if(!is.matrix(contrast.obj)) sqrt(sum(weights) * rse)
    else sqrt(rse * colSums(weights))
}

se.contrast.aovlist <-
    function(object, contrast.obj,
             coef = contr.helmert(ncol(contrast))[, 1L],
             data = NULL, ...)
{
    contrast.weight.aovlist <- function(object, contrast)
    {
        e.qr <- attr(object, "error.qr")
        if(!inherits(e.qr, "qr"))
            stop("'object' does not include an error 'qr' component")
        c.qr <- qr.qty(e.qr, contrast)
        e.assign <- attr(e.qr$qr, "assign")
        n.object <- length(object)
        e.assign <- c(e.assign,
                      rep.int(n.object - 1L, nrow(c.qr) - length(e.assign)))
	res <- setNames(vector("list", n.object), names(object))
        for(j in seq_along(names(object))) {
            strata <- object[[j]]
            if(is.qr(strata$qr)) {
                scontrast <- c.qr[e.assign == (j - 1L), , drop = FALSE]
                effects <- as.matrix(qr.qty(strata$qr, scontrast))
                asgn <- strata$assign[strata$qr$pivot[1L:strata$rank]]
                uasgn <- unique(asgn)
                nm <- c("(Intercept)", attr(strata$terms, "term.labels"))
                res.i <-
                    matrix(0, length(uasgn), ncol(effects),
                           dimnames = list(nm[1L + uasgn], colnames(contrast)))
                for(i in seq_along(uasgn)) {
                    select <- (asgn == uasgn[i])
                    res.i[i, ] <-
                        colSums(effects[seq_along(asgn)[select], , drop = FALSE]^2)
                }
                res[[j]] <- res.i
            }
        }
        res
    }
    SS <- function(aov.object)
    {
        rdf <- aov.object$df.residual
        if(is.null(rdf)) {
            nobs <- length(aov.object$residuals)
            rank <- aov.object$rank
            rdf <- nobs - rank
        }
        if(rdf == 0L) return(NA)
        resid <- as.matrix(aov.object$residuals)
        wt <- aov.object$weights
        if(!is.null(wt)) resid <- resid * sqrt(wt)
        sum(resid^2)/rdf
    }
    if(is.null(attr(object, "error.qr"))) {
        message("Refitting model to allow projection")
        object <- update(object, qr = TRUE)
    }
    contrast.obj <-
        if(is.null(data)) eval(contrast.obj)
        else eval(substitute(contrast.obj), data, parent.frame())
    if(!is.matrix(contrast.obj)) {
        if(!missing(coef)) {
            if(sum(coef) != 0)
                stop("'coef' must define a contrast, i.e., sum to 0")
            if(length(coef) != length(contrast.obj))
                stop("'coef' must have same length as 'contrast.obj'")
        }
        contrast <-
            sapply(contrast.obj,
                   function(x) {
                       if(!is.logical(x))
                           stop(gettextf("each element of '%s' must be logical",
                                         substitute(contrasts.obj)),
                                domain = NA)
                       x/sum(x)
                   })
        if(!length(contrast) || all(is.na(contrast)))
            stop("the contrast defined is empty (has no TRUE elements)")
        contrast <- contrast %*% coef
    }
    else {
        contrast <- contrast.obj
        if(any(abs(colSums(contrast)) > 1e-8))
            stop("columns of 'contrast.obj' must define a contrast(sum to zero)")
        if(!length(colnames(contrast)))
            colnames(contrast) <- paste("Contrast", seq(ncol(contrast)))
    }
    weights <- contrast.weight.aovlist(object, contrast)
    weights <- weights[-match("(Intercept)", names(weights))]
    effic <- eff.aovlist(object)
    ## Need to identify the lowest stratum where each nonzero term appears
    eff.used <- apply(effic, 2L,
                      function(x, ind = seq_along(x)) {
                          temp <- (x > 0)
                          if(sum(temp) == 1) temp
                          else max(ind[temp]) == ind
                      })
    if(is.matrix(eff.used)) {
        strata.nms <- rownames(effic)[row(eff.used)[eff.used]]
        var.nms <- colnames(effic)[col(eff.used)[eff.used]]
    } else {
        strata.nms <- rownames(effic)
        var.nms <- colnames(effic)
    }
    rse.list <- sapply(object[unique(strata.nms)], SS)
    wgt <- matrix(0, nrow = length(var.nms), ncol = ncol(contrast),
                  dimnames = list(var.nms, colnames(contrast)))
    for(i in seq_along(var.nms))
        wgt[i, ] <- weights[[strata.nms[i]]][var.nms[i], , drop = FALSE]
    rse <- rse.list[strata.nms]
    eff <- effic[eff.used]
    drop(sqrt((rse/eff^2) %*% wgt))
}