| # File src/library/stats/R/model.tables.R |
| # Part of the R package, https://www.R-project.org |
| # |
| # Copyright 1998 B. D. Ripley |
| # Copyright (C) 1998-2015 The R Core Team |
| # |
| # 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/ |
| |
| model.tables <- function(x, ...) UseMethod("model.tables") |
| |
| model.tables.aov <- function(x, type = "effects", se = FALSE, cterms, ...) |
| { |
| if(inherits(x, "maov")) |
| stop("'model.tables' is not implemented for multiple responses") |
| type <- match.arg(type, c("effects", "means", "residuals")) |
| if(type == "residuals") |
| stop(gettextf("type '%s' is not implemented yet", type), domain = NA) |
| prjs <- proj(x, unweighted.scale = TRUE) |
| if(is.null(x$call)) stop("this fit does not inherit from \"lm\"") |
| mf <- model.frame(x) |
| factors <- attr(prjs, "factors") |
| nf <- names(factors) |
| dn.proj <- setNames(as.list(nf), nf) |
| m.factors <- factors |
| t.factor <- attr(prjs, "t.factor") |
| vars <- colnames(t.factor) |
| which <- match(vars, names(dn.proj)) |
| which <- which[!is.na(which)] |
| dn.proj <- dn.proj[which] |
| m.factors <- m.factors[which] |
| ## with cterms, can specify subset of tables by name |
| if(!missing(cterms)) { |
| if(any(is.na(match(cterms, names(factors))))) |
| stop("'cterms' argument must match terms in model object") |
| dn.proj <- dn.proj[cterms] |
| m.factors <- m.factors[cterms] |
| } |
| if(type == "means") { |
| dn.proj <- |
| lapply(dn.proj, |
| function(x, mat, vn) |
| c("(Intercept)", |
| vn[(t(mat) %*% (as.logical(mat[, x]) - 1)) == 0]), |
| t.factor, vars) |
| } |
| tables <- make.tables.aovproj(dn.proj, m.factors, prjs, mf) |
| |
| ## This was reordering some interaction terms, e.g. N + V:N |
| ##n <- replications(paste("~", paste(names(tables), collapse = "+")), |
| ## data = mf) |
| n <- NULL |
| for(xx in names(tables)) n <- c(n, replications(paste("~", xx), data=mf)) |
| if(se) |
| if(is.list(n)) { |
| message("Design is unbalanced - use se.contrast() for se's") |
| se <- FALSE |
| } else se.tables <- se.aov(x, n, type = type) |
| if(type == "means" && "(Intercept)" %in% colnames(prjs)) { |
| gmtable <- mean(prjs[,"(Intercept)"]) |
| class(gmtable) <- "mtable" |
| tables <- c("Grand mean" = gmtable, tables) |
| } |
| result <- list(tables = tables, n = n) |
| if(se) result$se <- se.tables |
| attr(result, "type") <- type |
| class(result) <- c("tables_aov", "list.of") |
| result |
| } |
| |
| se.aov <- function(object, n, type = "means") |
| { |
| ## for balanced designs only |
| rdf <- object$df.residual |
| rse <- sqrt(sum(object$residuals^2)/rdf) |
| if(type == "effects") result <- rse/sqrt(n) |
| if(type == "means") |
| result <- |
| lapply(n, |
| function(x, d) { |
| nn <- unique(x) |
| nn <- nn[!is.na(nn)] |
| mat <- outer(nn, nn, function(x, y) 1/x + 1/y) |
| dimnames(mat) <- list(paste(nn), paste(nn)) |
| d * sqrt(mat) |
| }, d=rse) |
| attr(result, "type") <- type |
| class(result) <- "mtable" |
| result |
| } |
| |
| |
| model.tables.aovlist <- function(x, type = "effects", se = FALSE, ...) |
| { |
| type <- match.arg(type, c("effects", "means", "residuals")) |
| if(type == "residuals") |
| stop(gettextf("type '%s' is not implemented yet", type), domain = NA) |
| prjs <- proj(x, unweighted.scale = TRUE) |
| mf <- model.frame.aovlist(x) |
| factors <- lapply(prjs, attr, "factors") |
| dn.proj <- unlist(lapply(factors, names), recursive = FALSE) |
| m.factors <- unlist(factors, recursive = FALSE) |
| dn.strata <- rep.int(names(factors), lengths(factors)) |
| names(dn.strata) <- names(m.factors) <- names(dn.proj) <- unlist(dn.proj) |
| t.factor <- attr(prjs, "t.factor") |
| efficiency <- FALSE |
| if(type == "effects" || type == "means") { |
| if(anyDuplicated(names(dn.proj)[names(dn.proj) != "Residuals"])) { |
| efficiency <- eff.aovlist(x) |
| ## Elect to use the effects from the lowest stratum: |
| ## usually expect this to be highest efficiency |
| eff.used <- apply(efficiency, 2L, |
| function(x, ind = seq_len(x)) { |
| temp <- (x > 0) |
| if(sum(temp) == 1) temp |
| else max(ind[temp]) == ind |
| }) |
| } |
| } |
| if(any(efficiency)) { |
| if(is.list(eff.used)) |
| stop("design is unbalanced so cannot proceed") |
| which <- match(outer(rownames(efficiency), |
| colnames(efficiency), paste)[eff.used], |
| paste(dn.strata, dn.proj)) |
| efficiency <- efficiency[eff.used] |
| } else which <- match(colnames(t.factor), names(dn.proj)) |
| which <- which[!is.na(which)] |
| dn.proj <- dn.proj[which] |
| dn.strata <- dn.strata[which] |
| m.factors <- m.factors[which] |
| if(type == "means") { |
| t.factor <- t.factor[, names(dn.proj), drop = FALSE] |
| dn.proj <- |
| lapply(dn.proj, |
| function(x, mat, vn) |
| vn[(t(mat) %*% (as.logical(mat[, x]) - 1)) == 0], |
| t.factor, colnames(t.factor)) |
| } |
| tables <- |
| if(any(efficiency)) { |
| names(efficiency) <- names(dn.proj) |
| make.tables.aovprojlist(dn.proj, dn.strata, m.factors, prjs, mf, |
| efficiency) |
| } |
| else make.tables.aovprojlist(dn.proj, dn.strata, m.factors, prjs, mf) |
| if(type == "means") { |
| gmtable <- mean(prjs[["(Intercept)"]]) |
| class(gmtable) <- "mtable" |
| tables <- lapply(tables, "+", gmtable) |
| tables <- c("Grand mean" = gmtable, tables) |
| } |
| # n <- replications(attr(x, "call"), data = mf) |
| n <- replications(terms(x), data = mf) |
| if(se) |
| if(type == "effects" && is.list(n)) { |
| message("Standard error information not returned as design is unbalanced. \nStandard errors can be obtained through 'se.contrast'.") |
| se <- FALSE |
| } else if(type != "effects") { |
| warning(gettextf("SEs for type '%s' are not yet implemented", |
| type), domain = NA) |
| se <- FALSE |
| } else { |
| se.tables <- se.aovlist(x, dn.proj, dn.strata, factors, mf, |
| efficiency, n, type = type) |
| } |
| result <- list(tables = tables, n = n) |
| if(se) result$se <- se.tables |
| attr(result, "type") <- type |
| class(result) <- c("tables_aov", "list.of") |
| result |
| } |
| |
| se.aovlist <- function(object, dn.proj, dn.strata, factors, mf, efficiency, n, |
| type = "diff.means", ...) |
| { |
| if(type != "effects") |
| stop(gettextf("SEs for type '%s' are not yet implemented", type), |
| domain = NA) |
| RSS <- sapply(object, function(x) sum(x$residuals^2)/x$df.residual) |
| res <- vector(length = length(n), mode = "list") |
| names(res) <- names(n) |
| for(i in names(n)) { |
| sse <- RSS[[dn.strata[dn.proj[[i]]]]] |
| if(any(efficiency)) |
| sse <- sse/efficiency[i] |
| res[[i]] <- as.vector(sqrt(sse/n[i])) |
| class(res[[i]]) <- "mtable" |
| } |
| attr(res, "type") <- type |
| res |
| } |
| |
| |
| make.tables.aovproj <- |
| function(proj.cols, mf.cols, prjs, mf, fun = "mean", prt = FALSE, ...) |
| { |
| tables <- setNames(vector("list", length(proj.cols)), names(proj.cols)) |
| for(i in seq_along(tables)) { |
| terms <- proj.cols[[i]] |
| terms <- terms[terms %in% colnames(prjs)] |
| data <- |
| if(length(terms) == 1L) prjs[, terms] |
| else prjs[, terms] %*% as.matrix(rep.int(1, length(terms))) |
| tables[[i]] <- tapply(data, mf[mf.cols[[i]]], |
| get(fun, mode="function")) |
| class(tables[[i]]) <- "mtable" |
| if(prt) print(tables[i], ..., quote = FALSE) |
| } |
| tables |
| } |
| |
| |
| make.tables.aovprojlist <- |
| function(proj.cols, strata.cols, model.cols, projections, model, eff, |
| fun = "mean", prt = FALSE, ...) |
| { |
| tables <- setNames(vector("list", length(proj.cols)), names(proj.cols)) |
| if(!missing(eff)) { |
| for(i in seq_along(tables)) { |
| terms <- proj.cols[[i]] |
| if(all(is.na(eff.i <- match(terms, names(eff))))) |
| eff.i <- rep.int(1, length(terms)) |
| if(length(terms) == 1L) |
| data <- projections[[strata.cols[i]]][, terms]/ eff[eff.i] |
| else { |
| if(length(strata <- unique(strata.cols[terms])) == 1L) |
| data <- projections[[strata]][, terms] %*% |
| as.matrix(1/eff[eff.i]) |
| else { |
| mat <- NULL |
| for(j in strata) { |
| mat <- cbind(mat, projections[[j]][, terms[!is.na(match(terms, |
| names(strata.cols)[strata.cols == j]))]]) |
| } |
| data <- mat %*% as.matrix(1/eff[eff.i]) |
| } |
| } |
| tables[[i]] <- tapply(data, model[model.cols[[i]]], |
| get(fun, mode="function")) |
| attr(tables[[i]], "strata") <- strata.cols[i] |
| class(tables[[i]]) <- "mtable" |
| if(prt) print(tables[i], ..., quote = FALSE) |
| } |
| } else for(i in seq_along(tables)) { |
| terms <- proj.cols[[i]] |
| if(length(terms) == 1L) data <- projections[[strata.cols[i]]][, terms] |
| else { |
| if(length(strata <- unique(strata.cols[terms])) == 1L) |
| data <- projections[[strata]][, terms] %*% |
| as.matrix(rep.int(1, length(terms))) |
| else { |
| mat <- NULL |
| for(j in strata) { |
| mat <- cbind(mat, projections[[j]][, terms[!is.na(match(terms, |
| names(strata.cols)[strata.cols == j]))]]) |
| } |
| data <- mat %*% as.matrix(rep.int(1, length(terms))) |
| } |
| } |
| tables[[i]] <- tapply(data, model[model.cols[[i]]], get(fun)) |
| attr(tables[[i]], "strata") <- strata.cols[i] |
| class(tables[[i]]) <- "mtable" |
| if(prt) print(tables[i], ..., quote = FALSE) |
| } |
| tables |
| } |
| |
| replications <- function(formula, data = NULL, na.action) |
| { |
| if(missing(data) && inherits(formula, "data.frame")) { |
| data <- formula |
| formula <- ~ . |
| } |
| if(!inherits(formula, "terms")) { |
| formula <- as.formula(formula) |
| if(length(formula) < 3L) { |
| f <- y ~ x |
| f[[3L]] <- formula[[2L]] |
| formula <- f |
| } |
| formula <- terms(formula, data = data) |
| } |
| if(missing(na.action)) |
| if(!is.null(tj <- attr(data, "na.action")) && is.function(tj)) |
| na.action <- tj |
| else { |
| naa <- getOption("na.action") |
| if(!is.null(naa)) na.action <- match.fun(naa) |
| else na.action <- na.fail |
| } |
| f <- attr(formula, "factors") |
| o <- attr(formula, "order") |
| labels <- attr(formula, "term.labels") |
| vars <- as.character(attr(formula, "variables"))[-1L] |
| if(is.null(data)) { |
| v <- c(quote(data.frame), attr(formula, "variables")) |
| data <- eval(as.call(v), parent.frame()) |
| } |
| if(!is.function(na.action)) stop("na.action must be a function") |
| data <- na.action(data) |
| class(data) <- NULL |
| n <- length(o) |
| z <- setNames(vector("list", n), labels) |
| dummy <- numeric(.row_names_info(data, 2L)) |
| data <- lapply(data, function(x) if (is.character(x)) as.factor(x) else x) |
| notfactor <- !sapply(data, function(x) inherits(x, "factor")) |
| balance <- TRUE |
| for(i in seq_len(n)) { |
| l <- labels[i] |
| if(o[i] < 1 || substring(l, 1L, 5L) == "Error") { z[[l]] <- NULL; next } |
| select <- vars[f[, i] > 0] |
| if(any(nn <- notfactor[select])) { |
| warning(gettextf("non-factors ignored: %s", |
| paste(names(nn), collapse = ", ")), |
| domain = NA) |
| next |
| } |
| if(length(select)) |
| tble <- tapply(dummy, unclass(data[select]), length) |
| nrep <- unique(as.vector(tble)) |
| if(length(nrep) > 1L) { |
| balance <- FALSE |
| tble[is.na(tble)] <- 0 |
| z[[l]] <- tble |
| } else z[[l]] <- as.vector(nrep) |
| } |
| if(balance) unlist(z) else z |
| } |
| |
| print.tables_aov <- function(x, digits = 4L, ...) |
| { |
| tables.aov <- x$tables |
| n.aov <- x$n |
| se.aov <- if(se <- !is.na(match("se", names(x)))) x$se |
| type <- attr(x, "type") |
| switch(type, |
| effects = cat("Tables of effects\n"), |
| means = cat("Tables of means\n"), |
| residuals = if(length(tables.aov) > 1L) cat( |
| "Table of residuals from each stratum\n")) |
| if(!is.na(ii <- match("Grand mean", names(tables.aov)))) { |
| cat("Grand mean\n") |
| gmtable <- tables.aov[[ii]] |
| print.mtable(gmtable, digits = digits, ...) |
| } |
| for(i in names(tables.aov)) { |
| if(i == "Grand mean") next |
| table <- tables.aov[[i]] |
| cat("\n", i, "\n") |
| if(!is.list(n.aov)) |
| print.mtable(table, digits = digits, ...) |
| else { |
| n <- n.aov[[i]] |
| if(length(dim(table)) < 2L) { |
| table <- rbind(table, n) |
| rownames(table) <- c("", "rep") |
| print(table, digits = digits, ...) |
| } else { |
| ctable <- array(c(table, n), dim = c(dim(table), 2L)) |
| dim.t <- dim(ctable) |
| d <- length(dim.t) |
| ctable <- aperm(ctable, c(1, d, 2:(d - 1))) |
| dim(ctable) <- c(dim.t[1L] * dim.t[d], dim.t[-c(1, d)]) |
| dimnames(ctable) <- |
| c(list(format(c(rownames(table), rep.int("rep", dim.t[1L])))), |
| dimnames(table)[-1L]) |
| ctable <- eval(parse(text = paste( |
| "ctable[as.numeric(t(matrix(seq(nrow(ctable)),ncol=2)))", paste(rep.int(", ", d - 2), collapse = " "), "]"), |
| keep.source = FALSE)) |
| names(dimnames(ctable)) <- names(dimnames(table)) |
| class(ctable) <- "mtable" |
| print.mtable(ctable, digits = digits, ...) |
| } |
| } |
| } |
| if(se) { |
| if(type == "residuals") rn <- "df" else rn <- "replic." |
| switch(attr(se.aov, "type"), |
| effects = cat("\nStandard errors of effects\n"), |
| means = cat("\nStandard errors for differences of means\n"), |
| residuals = cat("\nStandard errors of residuals\n")) |
| if(length(unlist(se.aov)) == length(se.aov)) { |
| ## the simplest case: single replication, unique se |
| # kludge for NA's |
| n.aov <- n.aov[!is.na(n.aov)] |
| se.aov <- unlist(se.aov) |
| cn <- names(se.aov) |
| se.aov <- rbind(format(se.aov, digits = digits), format(n.aov)) |
| dimnames(se.aov) <- list(c(" ", rn), cn) |
| print(se.aov, quote=FALSE, right=TRUE, ...) |
| } else for(i in names(se.aov)) { |
| se <- se.aov[[i]] |
| if(length(se) == 1L) { ## single se |
| se <- rbind(se, n.aov[i]) |
| dimnames(se) <- list(c(i, rn), "") |
| print(se, digits = digits, ...) |
| } else { ## different se |
| dimnames(se)[[1L]] <- "" |
| cat("\n", i, "\n") |
| cat("When comparing means with same levels of:\n") |
| print(se, digits, ...) |
| cat("replic.", n.aov[i], "\n") |
| } |
| } |
| } |
| invisible(x) |
| } |
| |
| eff.aovlist <- function(aovlist) |
| { |
| Terms <- terms(aovlist) |
| if(names(aovlist)[[1L]] == "(Intercept)") aovlist <- aovlist[-1L] |
| pure.error.strata <- sapply(aovlist, function(x) is.null(x$qr)) |
| aovlist <- aovlist[!pure.error.strata] |
| s.labs <- names(aovlist) |
| ## find which terms are in which strata |
| s.terms <- |
| lapply(aovlist, function(x) { |
| asgn <- x$assign[x$qr$pivot[1L:x$rank]] |
| attr(terms(x), "term.labels")[asgn] |
| }) |
| t.labs <- attr(Terms, "term.labels") |
| t.labs <- t.labs[t.labs %in% unlist(s.terms)] |
| eff <- matrix(0, ncol = length(t.labs), nrow = length(s.labs), |
| dimnames = list(s.labs, t.labs)) |
| for(i in names(s.terms)) eff[i, s.terms[[i]] ] <- 1 |
| cs <- colSums(eff) |
| ## if all terms are in just one stratum we are done |
| if(all(cs <= 1)) return(eff[, cs > 0, drop = FALSE]) |
| |
| nm <- t.labs[ cs > 1] |
| pl <- |
| lapply(aovlist, function(x) |
| { |
| asgn <- x$assign[x$qr$pivot[1L:x$rank]] |
| sp <- split(seq_along(asgn), attr(terms(x), "term.labels")[asgn]) |
| sp <- sp[names(sp) %in% nm] |
| sapply(sp, function(x, y) { |
| y <- y[x, x, drop = FALSE] |
| res <- sum(diag(y)^2) |
| if(nrow(y) > 1 && sum(y^2) > 1.01 * res) |
| stop("eff.aovlist: non-orthogonal contrasts would give an incorrect answer") |
| res |
| }, y=x$qr$qr) |
| }) |
| for(i in names(pl)) eff[i, names(pl[[i]]) ] <- pl[[i]] |
| cs <- colSums(eff) |
| eff <- eff/rep(cs, each = nrow(eff)) |
| eff[, cs != 0, drop = FALSE] |
| } |
| |
| |
| model.frame.aovlist <- function(formula, data = NULL, ...) |
| { |
| ## formula is an aovlist object |
| call <- match.call() |
| oc <- attr(formula, "call") |
| Terms <- attr(formula, "terms") |
| rm(formula) |
| indError <- attr(Terms, "specials")$Error |
| errorterm <- attr(Terms, "variables")[[1 + indError]] |
| form <- update.formula(Terms, |
| paste(". ~ .-", deparse(errorterm, width.cutoff=500L, |
| backtick = TRUE), |
| "+", deparse(errorterm[[2L]], width.cutoff=500L, |
| backtick = TRUE))) |
| nargs <- as.list(call) |
| oargs <- as.list(oc) |
| nargs <- nargs[match(c("data", "na.action", "subset"), names(nargs), 0)] |
| args <- oargs[match(c("data", "na.action", "subset"), names(oargs), 0)] |
| args[names(nargs)] <- nargs |
| args$formula <- form |
| env <- environment(Terms) |
| if (is.null(env)) env <- parent.frame() |
| ## need stats:: for non-standard evaluation |
| fcall <- c(list(quote(stats::model.frame)), args) |
| eval(as.call(fcall), env) |
| } |
| |
| print.mtable <- |
| function(x, ..., digits = getOption("digits"), quote = FALSE, right = FALSE) |
| { |
| xxx <- x |
| xx <- attr(x, "Notes") |
| # nn <- names(dimnames(x)) |
| a.ind <- match(names(a <- attributes(x)), c("dim", "dimnames", "names")) |
| a <- a[!is.na(a.ind)] |
| class(x) <- attributes(x) <- NULL |
| attributes(x) <- a |
| # if(length(nn) > 1L) |
| # cat(paste("Dim ",paste(seq(length(nn)), "=", nn, collapse= ", "),"\n")) |
| if(length(x) == 1 && is.null(names(x)) && is.null(dimnames(x))) |
| names(x) <- rep("", length(x)) |
| if(length(dim(x)) && is.numeric(x)) { |
| xna <- is.na(x) |
| x <- format(zapsmall(x, digits)) |
| x[xna] <- " " |
| } |
| print(x, quote = quote, right = right, ...) |
| if(length(xx)) { |
| cat("\nNotes:\n") |
| print(xx) |
| } |
| invisible(xxx) |
| } |
| |
| |
