[R] multiple comparison function

[KAREN KOTSCHY]

Hi

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.

Cheers
Karen

# 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)  
   
   return(noquote(output))
}

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

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