24

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.

amoeba
  • 104,745
mark999
  • 3,258

2 Answers2

17

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.

Community
  • 1
Aaron left Stack Overflow
  • 7,198
  • Aaron, I don't think that your lmer code is correct. The OPs aov call is simply a standard repeated-measures design, which one would analyze with lmer as lmer(Y~A*B+(1|subject)). (Though see also this answer for more complicated models that permit estimation of across-Ss effect variance and correlations: http://stats.stackexchange.com/questions/13166/rs-lmer-cheat-sheet/13173#13173) – Mike Lawrence Aug 03 '11 at 22:52
  • 4
    The OP's aov call has three random effects; to replicate that with lmer my above code is correct. Your lmer code only has one random effect. Which is correct will depend on the context. – Aaron left Stack Overflow Aug 04 '11 at 02:45
  • Also note that the answer you linked to does not have any examples of crossed random effects. – Aaron left Stack Overflow Aug 04 '11 at 02:54
6

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]

John
  • 22,838
  • 9
  • 53
  • 87
  • 1
    I believe that if everything is balanced the result using a mixed model (that is, getting estimates with ML or REML) and the result using an ANOVA (that is, getting estimates with moments) will be nearly identical. The trouble in this case is the syntax to get the same model fit using both methods. – Aaron left Stack Overflow Aug 03 '11 at 18:48
  • I'm not sure what you're trying to accomplish. It looked like you were just trying to learn how to replicate results to better understand the relationship. What you want to do cannot be done with nlme. I just looked at lmer and it's not possible there either (although at least it reports the MS for your effects identically to the ANOVA). If you want ANOVA results, just do an ANOVA. With real data, done correctly, the two are almost never exactly the same. – John Aug 03 '11 at 18:58
  • 1
    Also note that the first attempt is not correct because it doesn't properly account for the crossed random effects. – Aaron left Stack Overflow Aug 03 '11 at 18:59
  • 1
    John, thank you for your answer. My reason for doing this was that I read somewhere on this site that repeated measures ANOVA is generally not recommended any more, with mixed-effects models being preferred. For some reason I was under the impression that the two methods would give the same results for a balanced design, and I was trying to confirm this. – mark999 Aug 03 '11 at 20:54
  • Aaron, I took it as the correct answer to what would be considered roughly equivalent. It's typically what's recommended as a first step in replicating repeated measures. There is no 'correct' as in perfect match. Adding in more random effects won't solve the problem. I note that one of the answers you refer to recommends the solution you wrote. However, that's absolutely no different an ANOVA result (the model is different but not ANOVA) than what I said was correct. I suspect the author was trying to just match what the OP was asking but it's not a sensible model. – John Aug 03 '11 at 22:21
  • @mark999: I believe you are correct that the two methods give the same results for a balanced design. I was more sure of that before you asked, though. It may be that we haven't yet got the syntax right, but it could also be that only some of the results are the same. That is, maybe in cases like this only the estimates are the same but the inference could be different. I suggest that because I know the inference in mixed models is a somewhat open question, and maybe we're getting bit by that here. – Aaron left Stack Overflow Aug 04 '11 at 02:59
  • @John: Yes, it is roughly equivalent, and I'd do it also as a first step in replicating a repeated measures design. I'm however not convinced that there's no perfect match and also, I don't believe that the ANOVA-type inference is in general the same as in the model you wrote. (Though it may be in this data set as some of the estimated variances are zero.) – Aaron left Stack Overflow Aug 04 '11 at 03:01
  • You're right that the ANOVA type inference isn't correct for the model, especially for interactions. But if the effects come out similarly it's very similar for the ANOVA of the model. – John Aug 04 '11 at 07:23