[R] post hoc comparison in repeated measure

[Richard M. Heiberger]

This is how I would do it.


arraychip <- read.table("arraychip.dat", sep='\t',
                        header=T, row.names=1)
for (i in 2:4) arraychip[[i]] <- factor(arraychip[[i]])


## run aov:
fit <- aov(x ~ treatment + Time + Error(animal), data=arraychip)
summary(fit)

## single stratum for the same ANOVA
fit2 <- aov(x ~ treatment + animal + Time, data=arraychip)
summary(fit2)

trt.means <- model.tables(fit2, cterms="treatment", type="means")
treat.means <- as.vector(trt.means$tables$treatment)
treat.n <- as.vector(trt.means$n$treatment)
names(treat.means) <- levels(arraychip$treatment)


## this expression works in R and S-Plus
mi.mj <- outer(treat.means, treat.means, "-")
s2 <- summary(fit2)$"Mean Sq"[2]             ## S-Plus
s2 <- summary(fit2)[[1]]["animal","Mean Sq"] ## R
s2.n <- s2 / treat.n
si.sj <- sqrt(outer(s2.n, s2.n, "+"))
q.tukey <- qtukey(.95, 4 ,36) /sqrt(2)

mi.mj
lower <- mi.mj - q.tukey*si.sj
upper <- mi.mj + q.tukey*si.sj
lower
upper


## in S-Plus you can also use
multicomp.default(treat.means,
                  df.residual=summary(fit2)$Df[2],
                  vmat=diag(
                    ## rep(
                    summary(fit2)$"Mean Sq"[2]
                    ##   ,  length(trt.means))
                    / treat.n
                     ),
                  method="tukey")