| # File src/library/stats/R/TukeyHSD.R |
| # Part of the R package, https://www.R-project.org |
| # |
| # Copyright (C) 2000-2001 Douglas M. Bates |
| # Copyright (C) 2002-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/ |
| |
| ### |
| ### Tukey multiple comparisons for R |
| ### |
| |
| TukeyHSD <- |
| function(x, which, ordered = FALSE, conf.level = 0.95, ...) |
| UseMethod("TukeyHSD") |
| |
| TukeyHSD.aov <- |
| function(x, which = seq_along(tabs), ordered = FALSE, |
| conf.level = 0.95, ...) |
| { |
| mm <- model.tables(x, "means") |
| if(is.null(mm$n)) |
| stop("no factors in the fitted model") |
| tabs <- mm$tables |
| if(names(tabs)[1L] == "Grand mean") tabs <- tabs[-1L] |
| tabs <- tabs[which] |
| ## mm$n need not be complete -- factors only -- so index by names |
| nn <- mm$n[names(tabs)] |
| nn_na <- is.na(nn) |
| if(all(nn_na)) |
| stop("'which' specified no factors") |
| if(any(nn_na)) { |
| warning("'which' specified some non-factors which will be dropped") |
| tabs <- tabs[!nn_na] |
| nn <- nn[!nn_na] |
| } |
| out <- setNames(vector("list", length(tabs)), names(tabs)) |
| MSE <- sum(x$residuals^2)/x$df.residual |
| for (nm in names(tabs)) { |
| tab <- tabs[[nm]] |
| means <- as.vector(tab) |
| nms <- if(length(dim(tab)) > 1L) { |
| dn <- dimnames(tab) |
| apply(do.call("expand.grid", dn), 1L, paste, collapse = ":") |
| } else names(tab) |
| n <- nn[[nm]] |
| ## expand n to the correct length if necessary |
| if (length(n) < length(means)) n <- rep.int(n, length(means)) |
| if (as.logical(ordered)) { |
| ord <- order(means) |
| means <- means[ord] |
| n <- n[ord] |
| if (!is.null(nms)) nms <- nms[ord] |
| } |
| center <- outer(means, means, "-") |
| keep <- lower.tri(center) |
| center <- center[keep] |
| width <- qtukey(conf.level, length(means), x$df.residual) * |
| sqrt((MSE/2) * outer(1/n, 1/n, "+"))[keep] |
| est <- center/(sqrt((MSE/2) * outer(1/n, 1/n, "+"))[keep]) |
| pvals <- ptukey(abs(est), length(means), x$df.residual, |
| lower.tail = FALSE) |
| dnames <- list(NULL, c("diff", "lwr", "upr","p adj")) |
| if (!is.null(nms)) dnames[[1L]] <- outer(nms, nms, paste, sep = "-")[keep] |
| out[[nm]] <- array(c(center, center - width, center + width,pvals), |
| c(length(width), 4L), dnames) |
| } |
| class(out) <- c("TukeyHSD", "multicomp") # multicomp is historical |
| attr(out, "orig.call") <- x$call |
| attr(out, "conf.level") <- conf.level |
| attr(out, "ordered") <- ordered |
| out |
| } |
| |
| print.TukeyHSD <- function(x, digits = getOption("digits"), ...) |
| { |
| cat(" Tukey multiple comparisons of means\n") |
| cat(" ", format(100*attr(x, "conf.level"), 2), |
| "% family-wise confidence level\n", sep = "") |
| if (attr(x, "ordered")) |
| cat(" factor levels have been ordered\n") |
| cat("\nFit: ", deparse(attr(x, "orig.call"), 500L), "\n\n", sep = "") |
| xx <- unclass(x) |
| attr(xx, "orig.call") <- attr(xx, "conf.level") <- attr(xx, "ordered") <- NULL |
| xx[] <- lapply(xx, function(z, digits) |
| {z[, "p adj"] <- round(z[, "p adj"], digits); z}, |
| digits = digits) |
| print.default(xx, digits, ...) |
| invisible(x) |
| } |
| |
| plot.TukeyHSD <- function (x, ...) |
| { |
| for (i in seq_along(x)) { |
| xi <- x[[i]][, -4L, drop = FALSE] # drop p-values |
| yvals <- nrow(xi):1L |
| dev.hold(); on.exit(dev.flush()) |
| ## xlab, main are set below, so block them from ... |
| plot(c(xi[, "lwr"], xi[, "upr"]), rep.int(yvals, 2L), type = "n", |
| axes = FALSE, xlab = "", ylab = "", main = NULL, ...) |
| axis(1, ...) |
| axis(2, at = nrow(xi):1, labels = dimnames(xi)[[1L]], srt = 0, ...) |
| abline(h = yvals, lty = 1, lwd = 0.5, col = "lightgray") |
| abline(v = 0, lty = 2, lwd = 0.5, ...) |
| segments(xi[, "lwr"], yvals, xi[, "upr"], yvals, ...) |
| segments(as.vector(xi), rep.int(yvals - 0.1, 3L), as.vector(xi), |
| rep.int(yvals + 0.1, 3L), ...) |
| title(main = paste0(format(100 * attr(x, "conf.level"), digits = 2L), |
| "% family-wise confidence level\n"), |
| xlab = paste("Differences in mean levels of", names(x)[i])) |
| box() |
| dev.flush(); on.exit() |
| } |
| } |