How to specify specific contrasts for repeated measures ANOVA using car?

Question

I am trying to run a repeated measures Anova in R followed by some specific contrasts on that dataset. I think the correct approach would be to use Anova() from the car package.

Lets illustrate my question with the example taken from ?Anova using the OBrienKaiser data (Note: I ommited the gender factor from the example):
We have a design with one between subjects factor, treatment (3 levels: control, A, B), and 2 repeated-measures (within subjects) factors, phase (3 levels: pretest, posttest, followup) and hour (5 levels: 1 to 5).

The standard ANOVA table is given by (in difference to example(Anova) I switched to Type 3 Sums of Squares, that is what my field wants):

require(car)
phase <- factor(rep(c("pretest", "posttest", "followup"), c(5, 5, 5)),
levels=c("pretest", "posttest", "followup"))
hour <- ordered(rep(1:5, 3))
idata <- data.frame(phase, hour)
mod.ok <- lm(cbind(pre.1, pre.2, pre.3, pre.4, pre.5, post.1, post.2, post.3, post.4, post.5, fup.1, fup.2, fup.3, fup.4, fup.5) ~ treatment, data=OBrienKaiser)
av.ok <- Anova(mod.ok, idata=idata, idesign=~phase*hour, type = 3)
summary(av.ok, multivariate=FALSE)

Now, imagine that the highest order interaction would have been significant (which is not the case) and we would like to explore it further with the following contrasts:
Is there a difference between hours 1&2 versus hours 3 (contrast 1) and between hours 1&2 versus hours 4&5 (contrast 2) in the treatment conditions (A&B together)?
In other words, how do I specify these contrasts:

1. ((treatment %in% c("A", "B")) & (hour %in% 1:2)) versus ((treatment %in% c("A", "B")) & (hour %in% 3))
2. ((treatment %in% c("A", "B")) & (hour %in% 1:2)) versus ((treatment %in% c("A", "B")) & (hour %in% 4:5))

My idea would be to run another ANOVA ommitting the non-needed treatment condition (control):

mod2 <- lm(cbind(pre.1, pre.2, pre.3, pre.4, pre.5, post.1, post.2, post.3, post.4, post.5, fup.1, fup.2, fup.3, fup.4, fup.5) ~ treatment, data=OBrienKaiser, subset = treatment != "control")
av2 <- Anova(mod2, idata=idata, idesign=~phase*hour, type = 3)
summary(av2, multivariate=FALSE)

However, I still have no idea how to set up the appropriate within-subject contrast matrix comparing hours 1&2 with 3 and 1&2 with 4&5. And I am not sure if omitting the non-needed treatment group is indeed a good idea as it changes the overall error term.

Before going for Anova() I was also thinking going for lme. However, there are small differences in F and p values between textbook ANOVA and what is returned from anove(lme) due to possible negative variances in standard ANOVA (which are not allowed in lme). Relatedly, somebody pointed me to gls which allows for fitting repeated measures ANOVA, however, it has no contrast argument.

To clarify: I want an F or t test (using type III sums of squares) that answers whether or not the desired contrasts are significant or not.

---

Update:

I already asked a very similar question on R-help, there was no answer.

A similar questions was posed on R-help some time ago. However, the answers did also not solve the problem.

---

Update (2015):

As this question still generates some activity, specifying theses and basically all other contrasts can now be done relatively easy with the afex package in combination with the lsmeans package as described in the afex vignette.

Answer 1

This method is generally considered "old-fashioned" so while it may be possible, the syntax is difficult and I suspect fewer people know how to manipulate the anova commands to get what you want. The more common method is using glht with a likelihood-based model from nlme or lme4. (I'm certainly welcome to be proved wrong by other answers though.)

That said, if I needed to do this, I wouldn't bother with the anova commands; I'd just fit the equivalent model using lm, pick out the right error term for this contrast, and compute the F test myself (or equivalently, t test since there's only 1 df). This requires everything to be balanced and have sphericity, but if you don't have that, you should probably be using a likelihood-based model anyway. You might be able to somewhat correct for non-sphericity using the Greenhouse-Geiser or Huynh-Feldt corrections which (I believe) use the same F statistic but modify the df of the error term.

If you really want to use car, you might find the heplot vignettes helpful; they describe how the matrices in the car package are defined.

Using caracal's method (for the contrasts 1&2 - 3 and 1&2 - 4&5), I get

      psiHat      tStat          F         pVal
1 -3.0208333 -7.2204644 52.1351067 2.202677e-09
2 -0.2083333 -0.6098777  0.3719508 5.445988e-01

This is how I'd get those same p-values:

Reshape the data into long format and run lm to get all the SS terms.

library(reshape2)
d <- OBrienKaiser
d$id <- factor(1:nrow(d))
dd <- melt(d, id.vars=c(18,1:2), measure.vars=3:17)
dd$hour <- factor(as.numeric(gsub("[a-z.]*","",dd$variable)))
dd$phase <- factor(gsub("[0-9.]*","", dd$variable), 
                   levels=c("pre","post","fup"))
m <- lm(value ~ treatment*hour*phase + treatment*hour*phase*id, data=dd)
anova(m)

Make an alternate contrast matrix for the hour term.

foo <- matrix(0, nrow=nrow(dd), ncol=4)
foo[dd$hour %in% c(1,2) ,1] <- 0.5
foo[dd$hour %in% c(3) ,1] <- -1
foo[dd$hour %in% c(1,2) ,2] <- 0.5
foo[dd$hour %in% c(4,5) ,2] <- -0.5
foo[dd$hour %in% 1 ,3] <- 1
foo[dd$hour %in% 2 ,3] <- 0
foo[dd$hour %in% 4 ,4] <- 1
foo[dd$hour %in% 5 ,4] <- 0

Check that my contrasts give the same SS as the default contrasts (and the same as from the full model).

anova(lm(value ~ hour, data=dd))
anova(lm(value ~ foo, data=dd))

Get the SS and df for just the two contrasts I want.

anova(lm(value ~ foo[,1], data=dd))
anova(lm(value ~ foo[,2], data=dd))

Get the p-values.

> F <- 73.003/(72.81/52)
> pf(F, 1, 52, lower=FALSE)
[1] 2.201148e-09
> F <- .5208/(72.81/52)
> pf(F, 1, 52, lower=FALSE)
[1] 0.5445999

Optionally adjust for sphericity.

pf(F, 1*.48867, 52*.48867, lower=FALSE)
pf(F, 1*.57413, 52*.57413, lower=FALSE)

Answer 2

If you want/have to use contrasts with the pooled error term from the corresponding ANOVA, you could do the following. Unfortunately, this will be long, and I don't know how to do this more conveniently. Still, I think the results are correct, as they are verified against Maxwell & Delaney (see below).

You want to compare groups of your first within factor hour in an SPF-p.qr design (notation from Kirk (1995): Split-Plot-Factorial design 1 between factor treatment with p groups, first within factor hour with q groups, second within factor prePostFup with r groups). The following assumes identically sized treatment groups and sphericity.

Nj    <- 10                                             # number of subjects per group
P     <- 3                                              # number of treatment groups
Q     <- 5                                              # number of hour groups
R     <- 3                                              # number of PrePostFup groups
id    <- factor(rep(1:(P*Nj), times=Q*R))                                  # subject
treat <- factor(rep(LETTERS[1:P], times=Q*R*Nj), labels=c("CG", "A", "B")) # treatment
hour  <- factor(rep(rep(1:Q, each=P*Nj), times=R))                         # hour
ppf   <- factor(rep(1:R, each=P*Q*Nj), labels=c("pre", "post", "fup"))     # prePostFup
DV    <- round(rnorm(Nj*P*Q*R, 15, 2), 2)               # some data with no effects
dfPQR <- data.frame(id, treat, hour, ppf, DV)           # data frame long format

summary(aov(DV ~ treat*hour*ppf + Error(id/(hour*ppf)), data=dfPQR)) # SPF-p.qr ANOVA

First note that the main effect for hour is the same after averaging over prePostFup, thus switching to the simpler SPF-p.q design which only contains treatment and hour as IVs.

dfPQ <- aggregate(DV ~ id + treat + hour, FUN=mean, data=dfPQR)  # average over ppf
# SPF-p.q ANOVA, note effect for hour is the same as before
summary(aov(DV ~ treat*hour + Error(id/hour), data=dfPQ))

Now note that in the SPF-p.q ANOVA, the effect for hour is tested against the interaction id:hour, i.e., this interaction provides the error term for the test. Now the contrasts for hour groups can be tested just like in a oneway between-subjects ANOVA by simply substituting the error term, and corresponding degrees of freedom. The easy way to get the SS and df of this interaction is to fit the model with lm().

(anRes <- anova(lm(DV ~ treat*hour*id, data=dfPQ)))
SSE    <- anRes["hour:id", "Sum Sq"]     # SS interaction hour:id -> will be error SS
dfSSE  <- anRes["hour:id", "Df"]         # corresponding df

But let's also calculate everything manually here.

# substitute DV with its difference to cell / person / treatment group means
Mjk   <- ave(dfPQ$DV,           dfPQ$treat, dfPQ$hour, FUN=mean)  # cell means
Mi    <- ave(dfPQ$DV, dfPQ$id,                         FUN=mean)  # person means
Mj    <- ave(dfPQ$DV,           dfPQ$treat,            FUN=mean)  # treatment means
dfPQ$IDxIV <- dfPQ$DV - Mi - Mjk + Mj                             # interaction hour:id
(SSE  <- sum(dfPQ$IDxIV^2))               # SS interaction hour:id -> will be error SS
dfSSE <- (Nj*P - P) * (Q-1)               # corresponding df
(MSE  <- SSE / dfSSE)                     # mean square

Now that we have the correct error term, we can build the usual test statistic for planned comparisons: $t = \frac{\hat{\psi} - 0}{||c|| \sqrt{MS_{E}}}$ where $c$ is the contrast vector, $||c||$ is its length, $\hat{\psi} = \sum\limits_{k=1}^{q} c_{k} M_{.k}$ is the contrast estimate, and $MS_{E}$ is the mean square for the hour:id interaction (the suitable error term).

Mj     <- tapply(dfPQ$DV, dfPQ$hour, FUN=mean)  # group means for hour
Nj     <- table(dfPQ$hour)                      # cell sizes for hour (here the same)
cntr   <- rbind(c(1, 1, -2,  0, 0),
                c(1, 1, -1, -1, 0))             # matrix of contrast vectors
psiHat <- cntr   %*% Mj                         # estimates psi-hat
lenSq  <- cntr^2 %*% (1/Nj)                     # squared lengths of contrast vectors
tStat  <- psiHat / sqrt(lenSq*MSE)              # t-statistics
pVal   <- 2*(1-pt(abs(tStat), dfSSE))           # p-values
data.frame(psiHat, tStat, pVal)

For multiple comparisons, you'd have to think about $\alpha$-correction methods, e.g., Bonferroni.

The corresponding calculations for Maxwell & Delaney's (2004) example on p. 599f can be found here. Note that M&D calculate the F-value, to see that the results are identical, you have to square the value for the t-statistic. That code also includes the analysis done with Anova() from car, as well as a manual calculation of the $\hat{\epsilon}$ corrections for the main effect of the within-factor.