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):

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.


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.

  • 1
    $\begingroup$ Have you already decided against using t-tests? What I mean is 1) throw out the data from the control group, 2) ignore the different levels of treatment, 3) for each person average over levels of prePostFup, 4) for each person average over hours 1,2 (= data group 1) as well as over hours 3,4 (= data group 2), 5) run t-test for 2 dependent groups. Since Maxwell & Delaney (2004) as well as Kirk (1995) discourage doing contrasts with a pooled error term in within-designs, this could be a simple alternative. $\endgroup$ – caracal Aug 29 '11 at 19:48
  • $\begingroup$ I would like to do contrast analyses and not pooled t tests. The reason is that contrasts (despite their problems) seem to be the standard procedure in Psychology and are what the readers/reviewer/supervisors want. Furthermore, they are relatively straight forward to do in SPSS. However, despite my 2 years as an active R user so far I have been unable to achieve it with R. Now I have to do some contrasts and I do not want to go back to SPSS just for this. When R is the future (which I think it is), contrasts must be possible. $\endgroup$ – Henrik Aug 30 '11 at 8:46

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.

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), 
m <- lm(value ~ treatment*hour*phase + treatment*hour*phase*id, data=dd)

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)
| cite | improve this answer | |
  • $\begingroup$ That also works! And thanks for the link to the heplots vignette, that really is a nice summary of what's going on in terms of the general linear model. $\endgroup$ – caracal Aug 30 '11 at 17:50
  • $\begingroup$ Thanks a lot. I will accept this answer (instead of the other great answer), as it includes some thought on sphericity correction. $\endgroup$ – Henrik Aug 31 '11 at 20:14
  • $\begingroup$ Note to future readers: The sphericity correction is equally applicable to the other solution as well. $\endgroup$ – Aaron left Stack Overflow Aug 31 '11 at 20:27

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.

| cite | improve this answer | |
  • $\begingroup$ Nice answer. This is more or less what I would have done if I had the patience to work it all out. $\endgroup$ – Aaron left Stack Overflow Aug 30 '11 at 15:38
  • $\begingroup$ Thanks for your detailed answer. Although it seems a little bit to unhandy in practice. $\endgroup$ – Henrik Aug 31 '11 at 20:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.