Repeated measures ANOVA with lme/lmer in R for two within-subject factors

Question

I'm trying to use lme from the nlme package to replicate results from aov for repeated measures ANOVAs. I've done this for a single-factor repeated measures experiment and for a two-factor experiment with one between-subjects factor and one within-subjects factor, but I'm having trouble doing it for a two-factor experiment with two within-subjects factors.

An example is shown below. A and B are fixed-effect factors and subject is a random-effect factor.

set.seed(1)
d <- data.frame(
    Y = rnorm(48),
    subject = factor(rep(1:12, 4)),
    A = factor(rep(1:2, each=24)),
    B = factor(rep(rep(1:2, each=12), 2)))

summary(aov(Y ~ A*B + Error(subject/(A*B)), data=d))  # Standard repeated measures ANOVA

library(nlme)
# Attempts:
anova(lme(Y ~ A*B, data=d, random = ~ 1 | subject))  # not same as above
anova(lme(Y ~ A*B, data=d, random = ~ 1 | subject/(A+B)))  # gives error

I could not see an explanation of this in the Pinheiro and Bates book, but I may have overlooked it.

Answer 1

What you're fitting with aov is called a strip plot, and it's tricky to fit with lme because the subject:A and subject:B random effects are crossed.

Your first attempt is equivalent to aov(Y ~ A*B + Error(subject), data=d), which doesn't include all the random effects; your second attempt is the right idea, but the syntax for crossed random effects using lme is very tricky.

Using lme from the nlme package, the code would be

lme(Y ~ A*B, random=list(subject=pdBlocked(list(~1, pdIdent(~A-1), pdIdent(~B-1)))), data=d)

Using lmer from the lme4 package, the code would be something like

lmer(Y ~ A*B + (1|subject) + (1|A:subject) + (1|B:subject), data=d)    

These threads from R-help may be helpful (and to give credit, that's where I got the nlme code from).

http://www.biostat.wustl.edu/archives/html/s-news/2005-01/msg00091.html

http://permalink.gmane.org/gmane.comp.lang.r.lme4.devel/3328

http://www.mail-archive.com/[email protected]/msg10843.html

This last link refers to p.165 of Pinheiro/Bates; that may be helpful too.

EDIT: Also note that in the data set you have, some of variance components are negative, which is not allowed using random effects with lme, so the results differ. A data set with all positive variance components can be created using a seed of 8. The results then agree. See this answer for details.

Also note that lme from nlme does not compute the denominator degrees of freedom correctly, so the F-statistics agree but not the p-values, and lmer from lme4 doesn't try too because it's very tricky in the presence of unbalanced crossed random effects, and may not even be a sensible thing to do. But that's more than I want to get into here.

Answer 2

Your first attempt is the correct answer if that's all you're trying to do. nlme() works out the between and within components, you don't need to specify them.

The problem you're running into isn't because you don't know how to specify the model, it's because repeated measures ANOVA and mixed effects are not the same thing. Sometimes the results from the ANOVA and mixed effects model will match. This is especially the case when you aggregate your data like you would for an ANOVA and calculate both from that. But generally, when done correctly, while the conclusions may be similar the results are almost never the same. Your example data aren't like real repeated measures where you often have replications of each measure within S. When you do an ANOVA typically you aggregate across those replications to get an estimate of the effect for each subject. In mixed effects modelling you do no such thing. You work with the raw data. When you do that you'll find that the results are never the same between ANOVA and lme().

[as an aside, using lmer() (from the lme4 package) instead of lme() give me SS and MS values that exactly match the ANOVA for effects in your example, it's just that the F's are different]