[R] multiple comparison function

KAREN KOTSCHY karen at gecko.biol.wits.ac.za
Tue Jan 4 16:28:54 CET 2000


I asked a while ago about multiple comparison tests for use with 
multi--way ANOVAs. Thanks to all those who replied and gave me 
some ideas and functions.

What I ended up doing was writing my own function for Tukey's 
HSD test (ref: Zar's Biostatistics pg 186-190), as this was the 
procedure I could best understand.

I've included my function here for those who may be interested. I'm 
sure there are many things that can be improved, but it works!
When interpreting the results one does, however, have to check 
which mean is which, as they are numbered according to their 
ranked order, not the order in which they were originally input.


# Tukey Multiple Comparison Test
# Caters for unequal sample sizes. Ref: ZAR pg   
# 186--190
# Arguments: x = vector of means to be compared 
#            ms = within-cells mean square (or   
#                 error MS for one-way ANOVA)
#            na, nb = total no. of data per      
#                     level of factor tested    
#                     (or sample sizes for one- 
#                      way)
#            df = within-cells degrees of        
#                 freedom (or overall df for one-
#                 way)

tukeymc <- function(x, ms, na, nb, df) {
   x <- sort(x)  # rank means
   se <- sqrt(ms/2 * (1/na + 1/nb))

   # Generate index vectors a and b for          
   # sequential subtraction of means
   y <- rev(1:length(x))
   y <- y[-length(x)]
   a <- rep(y, times = y-1)
   mb <- matrix(1, length(x)-1, length(x)-1)
   rb <- row(mb)
   b <- rb[rb + col(mb) <= length(x)]
   rm(y, mb, rb)

   # Compare diffs bt means with value from      
   # tables; alpha=0.05
   q <- (x[a] - x[b]) / se
   qt <- qtukey(0.95, length(x), df)
   results <- ifelse(q >= qt, "diff *", "same")

   # Generate output
   tab <- cbind(round(q, digits=4), results)
   rownames(tab) <- paste("mean", a, "vs",       
         "mean", b)
   rm(a, b)
   colnames(tab) <- c("q", "results")
   output <- list("Tukey's Multiple Comparison   
         Test","Ranked means" = x, "se" = se,    
      "qtables" = qt, tab)  

Karen Kotschy
Centre for Water in the Environment
University of the Witwatersrand

Tel: 011 716-2218
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

More information about the R-help mailing list