| # File src/library/stats/R/aov.R |
| # Part of the R package, https://www.R-project.org |
| # |
| # Copyright (C) 1995-2017 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) deparse(expr, width.cutoff = 500L, 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) |
| if(is.null(dny)) dny <- 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), 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)) |
| } |
