Repeated Measures ANOVA – Using LME/LMER in R for Two Within-Subject Factors

anovalme4-nlmemixed modelrrepeated measures

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.

Best Answer

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/r-help@stat.math.ethz.ch/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.

Related Question