[R] Are multiple range tests (Duncan, Bonferroni, Tukey, ...) available in R?

[Douglas Bates]

Ko-Kang Wang <Ko-Kang at xtra.co.nz> writes:

> Yes, even as an undergraduate student taking a paper on experimental
> design, I found this problem a while ago.
> 
> I am trying to use the current ANOVA table function and
> model.tables, extend them to get LSD, LSD_B, Tukey...etc directly.
> However since I've never written an R package before (done some
> single R functions though), it may take some time to do so.
> 
> By the way, have you tried model.tables()?  You can, I think
> sometimes anywa, get SED from this function.

Thanks for the hint of looking at model.tables().  Since V&R refer to
its only being available in S-PLUS 4.0 and later, I had neglected to
look for it in R.

That function does a lot of the heavy lifting for the multiple
comparisons calculations.  I had been promising my statistics class a
function to do Tukey multiple comparisons and I now have a draft
version enclosed below.  I would appreciate comments on it, especially
if I have the calculations wrong.  I think they are correct but I
would appreciate confirmation.

The function returns a set of confidence intervals on the differences
of the means for the levels of the factors.  If the optional argument
`order' is TRUE the levels of the factor are first put in increasing
order so the `diff' column of the result will be positive and you can
check which pairs of levels are significantly different by checking
which entries in the `lwr' column are positive.

The calculation is designed to accomodate unbalanced designs and
multiple factors.  It is not guaranteed that multiple comparisons will
always make sense in those cases. A sample usage is shown first.

 > library(Devore5)
 > data(xmp10.01)
 > fm1 <- aov(strength ~ type, data = xmp10.01)
 > summary(fm1)
	     Df Sum Sq Mean Sq F value    Pr(>F)    
 type         3 127375   42458  25.094 5.525e-07 ***
 Residuals   20  33839    1692                      
 ---
 Signif. codes:  0  `***'  0.001  `**'  0.01  `*'  0.05  `.'  0.1  ` '  1 
 > model.tables(fm1, "means", se = TRUE)
 Tables of means
 Grand mean

 682.5042 

  type 
     A     B     C     D 
 713.0 756.9 698.1 562.0 

 Standard errors for differences of means
	  type
	 23.75
 replic.     6
 > Tukey(fm1)
 $type
	   diff        lwr         upr
 B-A   43.93333  -22.53671  110.403377
 C-A  -14.93333  -81.40338   51.536711
 D-A -150.98333 -217.45338  -84.513289
 C-B  -58.86667 -125.33671    7.603377
 D-B -194.91667 -261.38671 -128.446623
 D-C -136.05000 -202.52004  -69.579956

 attr(,"class")
 [1] "Tukey"
 > Tukey(fm1, order = TRUE)
 $type
	  diff        lwr       upr
 C-D 136.05000  69.579956 202.52004
 A-D 150.98333  84.513289 217.45338
 B-D 194.91667 128.446623 261.38671
 A-C  14.93333 -51.536711  81.40338
 B-C  58.86667  -7.603377 125.33671
 B-A  43.93333 -22.536711 110.40338

 attr(,"class")
 [1] "Tukey"

### $Id: Tukey.R,v 1.1 2000/05/05 22:31:21 bates Exp $
###
###               Tukey multiple comparisons for R
###
### Copyright 2000-2000  Douglas M. Bates <bates at stat.wisc.edu>
###
### This file is made available under the terms of the GNU General
### Public License, version 2, or at your option, any later version,
### incorporated herein by reference.
###
### 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.
###
### You should have received a copy of the GNU General Public
### License along with this program; if not, write to the Free
### Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
### MA 02111-1307, USA

Tukey <- function(x, ...) UseMethod("Tukey")

Tukey.aov <-
    function(x, order = FALSE, conf.level = 0.95, ...)
{
    mm <- model.tables(x, "means")
    tabs <- mm$tables[-1]
    nn <- mm$n
    out <- vector("list", length(tabs))
    names(out) <- names(tabs)
    MSE <- sum(resid(x)^2)/x$df.residual
    for (nm in names(tabs)) {
        means <- as.vector(tabs[[nm]])
        nms <- names(tabs[[nm]])
        n <- nn[[nm]]
        ## expand n to the correct length if necessary
        if (length(n) < length(means)) n <- rep(n, length(means))
        if (as.logical(order)) {
            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]
        dnames <- list(NULL, c("diff", "lwr", "upr"))
        if (!is.null(nms)) dnames[[1]] <- outer(nms, nms, paste, sep = "-")[keep]
        out[[nm]] <- array(c(center, center - width, center + width),
                           c(length(width), 3), dnames)
    }
    class(out) <- "Tukey"
    out
}

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._