Skip to content

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

4 messages · Rodolfo Canet-Castello, Peter Dalgaard, Ko-Kang Kevin Wang +1 more

Rodolfo Canet-Castello
# 
Dear friends

I'm trying to switch to R for my current work, but could not find a 
more or less direct procedure to perform multiple range tests such 
as LSD, Duncan, Tukey, etc after an analysis of variance. I've 
checked out the documentation, but could not find any useful 
information for a non-very technical user. Is there any package 
offering prebuilt procedures to perform such tests? I would 
appreciate very much any hint.

I have not subscribed to the list just for this help application. 
Please Cc: your replies, if any, to my address. Thank you very 
much for your help.

*********************************************************
Dr. Rodolfo Canet-Castello
Instituto Valenciano de Investigaciones Agrarias (IVIA)
Dpto. Recursos Naturales
Aptdo. oficial. 46113-Moncada (Valencia). ESPA?A-SPAIN.
Phone: 34-96-1391000      Fax: 34-96-1390240
Web page: http://www.ivia.es
*********************************************************
Web page of Spanish Group of Soil Enzymology:
http://www.ivia.es/soilenzymology/
*********************************************************
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Peter Dalgaard
# 
"Rodolfo Canet-Castello" <rcanet at ivia.es> writes:
For one-way anova, pairwise.t.test will do the LSD (aka. "none"),
Bonferroni, Hochberg, and Holm (default) multiple testing adjustments.

I did have some view to other cases when writing it, but didn't get
around to implementing it (yet).

("Other cases" meaning 

	- Scheffe (does anyone ever use those?)
	- Tukey, Duncan, Newman-Keuls, REGW, ... for *balanced* anova
	  effects
        - pairwise testing of contrasts in general regression models

As far as my preliminary analysis went, the results of all of these
procedures can be represented by a table of adjusted p-values, but the
interface issues are not completely trivial.)
Ko-Kang Kevin Wang
# 
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.

Cheers,

Ko-Kang
Rodolfo Canet-Castello wrote:

            

      
      
      
    

        
--
-----------------------------------------------------------------------------------

 Ko-Kang Wang
 Undergraduate Student
 Computer Science/Statistics Double Major
 University of Auckland
 Auckland 1005
 New Zealand
-----------------------------------------------------------------------------------



-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

  
  


  


      

    

    

      
      

        


  

        

      Douglas Bates
    

    

      
      #
    

  


  

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

            

      
      
      
    

        
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._