#' Summarize an Agnostic Analysis of Variance Model
#'
#' @param object object of class "agnostic.aov".
#' @param intercept logical: should intercept terms be included?
#' @param split an optional named list, with names corresponding to terms in the model. Each component is itself a list with integer components giving contrasts whose contributions are to be summed.
#' @param expand.split logical: should the split apply also to interactions involving the factor?
#' @param keep.zero.df logical: should terms with no degrees of freedom be included?
#' @param ... Arguments to be passed to or from other methods.
#'
#' @return An object of class c("summary.aov", "listof") or "summary.aovlist" respectively.
#'
#' For fits with a single stratum the result will be a list of agnostic ANOVA tables, one for each response (even if there is only one response): the tables are of class "anova" inheriting from class "data.frame". They have columns "Df", "Sum Sq", "Mean Sq", as well as "F value", "P.H0" and "P.H1" if there are non-zero residual degrees of freedom. There is a row for each term in the model, plus one for "Residuals" if there are any.
#'
#' For multistratum fits the return value is a list of such summaries, one for each stratum.
#' @export
#'
#' @examples
#' #Test data
#' obs <- c(54, 60, 59, 45, 40, 55, 47, 33, 51, 66, 44, 34, 36, 61,
#' 49, 48, 50, 54, 62, 42, 48, 61, 60, 44)
#' obs2 <- rnorm(24)
#' groups <- as.factor(rep(c("A", "B", "c", "D"), 6))
#' groups2 <- as.factor(rep(c("A", "D", "B", "C"), 6))
#'
#' mod1 <- agnostic.aov(obs ~ groups)
#' summary(mod1)
#' mod2 <- agnostic.aov(cbind(obs, obs2) ~ groups * groups2)
#' summary(mod2)
summary.agnostic.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)) {
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$"P.H0" <- TP
x$"P.H1" <- pf(TT, df, rdf, lower.tail = TRUE)
}
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
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.