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