# What is the lme4::lmer equivalent of a three-way repeated measures ANOVA?

My question is based on this response which showed which lme4::lmer model corresponds to a two-way repeated measures ANOVA:

require(lme4)
set.seed(1234)
d <- data.frame(
y = rnorm(96),
subject = factor(rep(1:12, 4)),
a = factor(rep(1:2, each=24)),
b = factor(rep(rep(1:2, each=12))),
c = factor(rep(rep(1:2, each=48))))

# standard two-way repeated measures ANOVA:
summary(aov(y~a*b+Error(subject/(a*b)), d[d$c == "1",])) # corresponding lmer call: anova(lmer(y ~ a*b+(1|subject) + (1|a:subject) + (1|b:subject), d[d$c == "1",]))


My question now is on how to extend this to the case of a three-way ANOVA:

summary(aov(y~a*b*c+Error(subject/(a*b*c)), d))
## [...]
## Error: subject:a:b:c
##           Df Sum Sq Mean Sq F value Pr(>F)
## a:b:c      1  0.101  0.1014   0.115  0.741
## Residuals 11  9.705  0.8822


The natural extension as well as versions thereof do not match the ANOVA results:

anova(lmer(y ~ a*b*c +(1|subject) + (1|a:subject) + (1|b:subject) + (1|c:subject), d))
## [...]
## a:b:c  1 0.1014  0.1014  0.1500

anova(lmer(y ~ a*b*c +(1|subject) + (1|a:subject) + (1|b:subject) + (1|c:subject) +
(1|a:b:subject) + (1|a:c:subject) + (1|b:c:subject), d))
## [...]
## a:b:c  1 0.1014  0.1014  0.1539


Note that a very similar question has been asked before. However, it was missing example data (which is provided here).

• Are you sure that you don't want your two level model to be y ~ a*b + (1 + a*b|subject), d[d$c == "1",] ? Or perhaps I'm missing something? – Rasmus Bååth Oct 2 '14 at 20:04 • @RasmusBååth Go Ahead and try to fit it, lmer will complain as the random effects are not identified any more. Initially I also thought that this is the model I want, but it is not. If you compare the lmer model I propose for the 2-way case with the standard ANOVA you will see that the F-values exactly match. As said in the response I linked. – Henrik Oct 2 '14 at 20:07 • For the three-way problem, the first lmer model that you wrote (which excludes the random two-way interactions) is not expected to be equivalent to a 3-way RM-ANOVA, but the second one that you wrote (which includes the random two-way interactions) should be. As for why there is a discrepancy even with that model, I have a hunch about what the problem is, going to grab dinner then will look at the toy dataset some more. – Jake Westfall Nov 4 '14 at 1:13 ## 2 Answers The direct answer to your question is that the last model you wrote, anova(lmer(y ~ a*b*c +(1|subject) + (1|a:subject) + (1|b:subject) + (1|c:subject) + (1|a:b:subject) + (1|a:c:subject) + (1|b:c:subject), d))  I believe is "in principle" correct, although it is a strange parameterization that doesn't always seem to work well in actual practice. As for why the output that you get from this model is discrepant with the aov() output, I think there are two reasons. 1. Your simple simulated dataset is pathological in that the best-fitting model is one that implies negative variance components, which mixed models fit by lmer() (and most other mixed model programs) won't allow. 2. Even with a non-pathological dataset, the way you have the model set up, as mentioned above, does not always seem to work well in practice, although I must admit I don't really understand why. It is also just generally strange in my opinion, but that's another story. Let me first demonstrate the parameterization that I prefer on your initial two-way ANOVA example. Assume that your dataset d is loaded. Your model (note that I changed from dummy to contrast codes) was: options(contrasts=c("contr.sum","contr.poly")) mod1 <- lmer(y ~ a*b+(1|subject) + (1|a:subject) + (1|b:subject), data = d[d$c == "1",])
anova(mod1)
# Analysis of Variance Table
#     Df  Sum Sq Mean Sq F value
# a    1 2.20496 2.20496  3.9592
# b    1 0.13979 0.13979  0.2510
# a:b  1 1.23501 1.23501  2.2176


which worked fine here in that it matched the aov() output. The model that I prefer involves two changes: manually contrast-coding the factors so that we are not working with R factor objects (which I recommend doing in 100% of cases), and specifying the random effects differently:

d <- within(d, {
A <- 2*as.numeric(paste(a)) - 3
B <- 2*as.numeric(paste(b)) - 3
C <- 2*as.numeric(paste(c)) - 3
})
mod2 <- lmer(y ~ A*B + (1|subject)+(0+A|subject)+(0+B|subject),
data = d[d$c == "1",]) anova(mod2) # Analysis of Variance Table # Df Sum Sq Mean Sq F value # A 1 2.20496 2.20496 3.9592 # B 1 0.13979 0.13979 0.2510 # A:B 1 1.23501 1.23501 2.2176 logLik(mod1) # 'log Lik.' -63.53034 (df=8) logLik(mod2) # 'log Lik.' -63.53034 (df=8)  The two approaches are totally equivalent in the simple 2-way problem. Now we'll move to a 3-way problem. I mentioned earlier that the example dataset you gave was pathological. So what I want to do before addressing your example dataset is to first generate a dataset from an actual variance components model (i.e., where non-zero variance components are built into the true model). First I will show how my preferred parameterization seems to work better than the one you proposed. Then I will demonstrate another way of estimating the variance components which does not impose that they must be non-negative. Then we will be in a position to see the problem with the original example dataset. The new dataset will be identical in structure except we will have 50 subjects: set.seed(9852903) d2 <- expand.grid(A=c(-1,1), B=c(-1,1), C=c(-1,1), sub=seq(50)) d2 <- merge(d2, data.frame(sub=seq(50), int=rnorm(50), Ab=rnorm(50), Bb=rnorm(50), Cb=rnorm(50), ABb=rnorm(50), ACb=rnorm(50), BCb=rnorm(50))) d2 <- within(d2, { y <- int + (1+Ab)*A + (1+Bb)*B + (1+Cb)*C + (1+ABb)*A*B + (1+ACb)*A*C + (1+BCb)*B*C + A*B*C + rnorm(50*2^3) a <- factor(A) b <- factor(B) c <- factor(C) })  The F-ratios we want to match are: aovMod1 <- aov(y ~ a*b*c + Error(factor(sub)/(a*b*c)), data = d2) tab <- lapply(summary(aovMod1), function(x) x[[1]][1,2:4]) do.call(rbind, tab) # Sum Sq Mean Sq F value # Error: factor(sub) 439.48 8.97 # Error: factor(sub):a 429.64 429.64 32.975 # Error: factor(sub):b 329.48 329.48 27.653 # Error: factor(sub):c 165.44 165.44 17.924 # Error: factor(sub):a:b 491.33 491.33 49.694 # Error: factor(sub):a:c 305.46 305.46 41.703 # Error: factor(sub):b:c 466.09 466.09 40.655 # Error: factor(sub):a:b:c 392.76 392.76 448.101  Here are our two models: mod3 <- lmer(y ~ a*b*c + (1|sub)+(1|a:sub)+(1|b:sub)+(1|c:sub)+ (1|a:b:sub)+(1|a:c:sub)+(1|b:c:sub), data = d2) anova(mod3) # Analysis of Variance Table # Df Sum Sq Mean Sq F value # a 1 32.73 32.73 34.278 # b 1 21.68 21.68 22.704 # c 1 12.53 12.53 13.128 # a:b 1 60.93 60.93 63.814 # a:c 1 50.38 50.38 52.762 # b:c 1 57.30 57.30 60.009 # a:b:c 1 392.76 392.76 411.365 mod4 <- lmer(y ~ A*B*C + (1|sub)+(0+A|sub)+(0+B|sub)+(0+C|sub)+ (0+A:B|sub)+(0+A:C|sub)+(0+B:C|sub), data = d2) anova(mod4) # Analysis of Variance Table # Df Sum Sq Mean Sq F value # A 1 28.90 28.90 32.975 # B 1 24.24 24.24 27.653 # C 1 15.71 15.71 17.924 # A:B 1 43.56 43.56 49.694 # A:C 1 36.55 36.55 41.703 # B:C 1 35.63 35.63 40.655 # A:B:C 1 392.76 392.76 448.101 logLik(mod3) # 'log Lik.' -984.4531 (df=16) logLik(mod4) # 'log Lik.' -973.4428 (df=16)  As we can see, only the second method matches the output from aov(), although the first method is at least in the ballpark. The second method also achieves a higher log-likelihood. I am not sure why these two methods give different results, as again I think they are "in principle" equivalent, but maybe it is for some numerical/computational reasons. Or maybe I am mistaken and they are not equivalent even in principle. Now I will show another way of estimating the variance components based on traditional ANOVA ideas. Basically we will take the expected mean square equations for your design, substitute in the observed values of the mean squares, and solve for the variance components. To get the expected mean squares we will use an R function that I wrote a few years ago, called EMS(), which is documented HERE. Below I assume the function is loaded already. # prepare coefficient matrix r <- 1 # number of replicates s <- 50 # number of subjects a <- 2 # number of levels of A b <- 2 # number of levels of B c <- 2 # number of levels of C CT <- EMS(r ~ a*b*c*s, random="s") expr <- strsplit(CT[CT != ""], split="") expr <- unlist(lapply(expr, paste, collapse="*")) expr <- sapply(expr, function(x) eval(parse(text=x))) CT[CT != ""] <- expr CT[CT == ""] <- 0 mode(CT) <- "numeric" # residual variance and A*B*C*S variance are confounded in # this design, so remove the A*B*C*S variance component CT <- CT[-15,-2] CT # VarianceComponent # Effect e b:c:s a:c:s a:b:s a:b:c c:s b:s a:s b:c a:c a:b s c b a # a 1 0 0 0 0 0 0 4 0 0 0 0 0 0 200 # b 1 0 0 0 0 0 4 0 0 0 0 0 0 200 0 # c 1 0 0 0 0 4 0 0 0 0 0 0 200 0 0 # s 1 0 0 0 0 0 0 0 0 0 0 8 0 0 0 # a:b 1 0 0 2 0 0 0 0 0 0 100 0 0 0 0 # a:c 1 0 2 0 0 0 0 0 0 100 0 0 0 0 0 # b:c 1 2 0 0 0 0 0 0 100 0 0 0 0 0 0 # a:s 1 0 0 0 0 0 0 4 0 0 0 0 0 0 0 # b:s 1 0 0 0 0 0 4 0 0 0 0 0 0 0 0 # c:s 1 0 0 0 0 4 0 0 0 0 0 0 0 0 0 # a:b:c 1 0 0 0 50 0 0 0 0 0 0 0 0 0 0 # a:b:s 1 0 0 2 0 0 0 0 0 0 0 0 0 0 0 # a:c:s 1 0 2 0 0 0 0 0 0 0 0 0 0 0 0 # b:c:s 1 2 0 0 0 0 0 0 0 0 0 0 0 0 0 # e 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 # get mean squares (MSmod <- summary(aov(y ~ a*b*c*factor(sub), data=d2))) # Df Sum Sq Mean Sq # a 1 429.6 429.6 # b 1 329.5 329.5 # c 1 165.4 165.4 # factor(sub) 49 439.5 9.0 # a:b 1 491.3 491.3 # a:c 1 305.5 305.5 # b:c 1 466.1 466.1 # a:factor(sub) 49 638.4 13.0 # b:factor(sub) 49 583.8 11.9 # c:factor(sub) 49 452.2 9.2 # a:b:c 1 392.8 392.8 # a:b:factor(sub) 49 484.5 9.9 # a:c:factor(sub) 49 358.9 7.3 # b:c:factor(sub) 49 561.8 11.5 # a:b:c:factor(sub) 49 42.9 0.9 MS <- MSmod[[1]][,"Mean Sq"] # solve ans <- solve(CT, MS) cbind(rev(ans[c(grep("e",names(ans)),grep("s",names(ans)))])/ c(1,2,2,2,4,4,4,1)) # s 1.0115549 # a:s 1.5191114 # b:s 1.3797937 # c:s 1.0441351 # a:b:s 1.1263331 # a:c:s 0.8060402 # b:c:s 1.3235126 # e 0.8765093 summary(mod4) # Random effects: # Groups Name Variance Std.Dev. # sub (Intercept) 1.0116 1.0058 # sub.1 A 1.5191 1.2325 # sub.2 B 1.3798 1.1746 # sub.3 C 1.0441 1.0218 # sub.4 A:B 1.1263 1.0613 # sub.5 A:C 0.8060 0.8978 # sub.6 B:C 1.3235 1.1504 # Residual 0.8765 0.9362 # Number of obs: 400, groups: sub, 50  Okay, now we will return to the original example. The F-ratios we are trying to match are: aovMod2 <- aov(y~a*b*c+Error(subject/(a*b*c)), data = d) tab <- lapply(summary(aovMod2), function(x) x[[1]][1,2:4]) do.call(rbind, tab) # Sum Sq Mean Sq F value # Error: subject 13.4747 1.2250 # Error: subject:a 1.4085 1.4085 1.2218 # Error: subject:b 3.1180 3.1180 5.5487 # Error: subject:c 6.3809 6.3809 5.2430 # Error: subject:a:b 1.5706 1.5706 2.6638 # Error: subject:a:c 1.0907 1.0907 1.5687 # Error: subject:b:c 1.4128 1.4128 2.3504 # Error: subject:a:b:c 0.1014 0.1014 0.1149  Here are our two models: mod5 <- lmer(y ~ a*b*c + (1|subject)+(1|a:subject)+(1|b:subject)+ (1|c:subject)+(1|a:b:subject)+(1|a:c:subject)+(1|b:c:subject), data = d) anova(mod5) # Analysis of Variance Table # Df Sum Sq Mean Sq F value # a 1 0.8830 0.8830 1.3405 # b 1 3.1180 3.1180 4.7334 # c 1 3.8062 3.8062 5.7781 # a:b 1 1.5706 1.5706 2.3844 # a:c 1 0.9620 0.9620 1.4604 # b:c 1 1.4128 1.4128 2.1447 # a:b:c 1 0.1014 0.1014 0.1539 mod6 <- lmer(y ~ A*B*C + (1|subject)+(0+A|subject)+(0+B|subject)+ (0+C|subject)+(0+A:B|subject)+(0+A:C|subject)+ (0+B:C|subject), data = d) anova(mod6) # Analysis of Variance Table # Df Sum Sq Mean Sq F value # a 1 0.8830 0.8830 1.3405 # b 1 3.1180 3.1180 4.7334 # c 1 3.8062 3.8062 5.7781 # a:b 1 1.5706 1.5706 2.3844 # a:c 1 0.9620 0.9620 1.4604 # b:c 1 1.4128 1.4128 2.1447 # a:b:c 1 0.1014 0.1014 0.1539 logLik(mod5) # 'log Lik.' -135.0351 (df=16) logLik(mod6) # 'log Lik.' -134.9191 (df=16)  In this case the two models yield basically the same results, although the second method has a very slightly higher log-likelihood. Neither method matches aov(). But let's look at what we get when we solve for the variance components as we did above, using the ANOVA procedure that does not constrain variance components to be non-negative (but which can only be used in balanced designs with no continuous predictors and no missing data; the classical ANOVA assumptions). # prepare coefficient matrix r <- 1 # number of replicates s <- 12 # number of subjects a <- 2 # number of levels of A b <- 2 # number of levels of B c <- 2 # number of levels of C CT <- EMS(r ~ a*b*c*s, random="s") expr <- strsplit(CT[CT != ""], split="") expr <- unlist(lapply(expr, paste, collapse="*")) expr <- sapply(expr, function(x) eval(parse(text=x))) CT[CT != ""] <- expr CT[CT == ""] <- 0 mode(CT) <- "numeric" # residual variance and A*B*C*S variance are confounded in # this design, so remove the A*B*C*S variance component CT <- CT[-15,-2] # get mean squares MSmod <- summary(aov(y ~ a*b*c*subject, data=d)) MS <- MSmod[[1]][,"Mean Sq"] # solve ans <- solve(CT, MS) cbind(rev(ans[c(grep("e",names(ans)),grep("s",names(ans)))])/ c(1,2,2,2,4,4,4,1)) # s 0.04284033 # a:s 0.03381648 # b:s -0.04004005 # c:s 0.04184887 # a:b:s -0.03657940 # a:c:s -0.02337501 # b:c:s -0.03514457 # e 0.88224787 summary(mod6) # Random effects: # Groups Name Variance Std.Dev. # subject (Intercept) 7.078e-02 2.660e-01 # subject.1 A 6.176e-02 2.485e-01 # subject.2 B 0.000e+00 0.000e+00 # subject.3 C 6.979e-02 2.642e-01 # subject.4 A:B 1.549e-16 1.245e-08 # subject.5 A:C 4.566e-03 6.757e-02 # subject.6 B:C 0.000e+00 0.000e+00 # Residual 6.587e-01 8.116e-01 # Number of obs: 96, groups: subject, 12  Now we can see what is pathological about the original example. The best-fitting model is one that implies that several of the random variance components are negative. But lmer() (and most other mixed model programs) constrains the estimates of variance components to be non-negative. This is generally considered a sensible constraint, since variances can of course never truly be negative. However, a consequence of this constraint is that mixed models are unable to accurately represent datasets that feature negative intraclass correlations, that is, datasets where observations from the same cluster are less (rather than more) similar on average than observations drawn randomly from the dataset, and consequently where within-cluster variance substantially exceeds between-cluster variance. Such datasets are perfectly reasonable datasets that one will occasionally come across in the real world (or accidentally simulate!), but they cannot be sensibly described by a variance-components model, because they imply negative variance components. They can however be "non-sensibly" described by such models, if the software will allow it. aov() allows it. lmer() does not. • +1. Re I am not sure why these two methods give different results, as again I think they are "in principle" equivalent, but maybe it is for some numerical/computational reasons - do you perhaps understand this better now (two years later)? I tried to figure out what is the difference, but don't get it either... – amoeba Oct 10 '16 at 19:05 • @amoeba My current thinking is still pretty much the same as then: AFAIK, the two models are statistically equivalent (in the sense that they make the same predictions about the data and imply the same standard errors), even though the random effects are parameterized differently. I think the observed differences -- which seem to only occur sometimes -- are just due to computational issues. In particular, I suspect that you could fiddle around with optimizer settings (like varying starting points or using more strict convergence criteria) until the two models returned exactly the same answer. – Jake Westfall Oct 10 '16 at 21:18 • Thanks for your reply. I am rather unconvinced: I tried to fiddle with optimizer settings and could not change anything in the outcomes; my impression is that both models are well converged. I might ask this is a separate question at some moment. – amoeba Oct 12 '16 at 9:43 • I keep investigating this issue and realized the following: as soon as repeated-measures factors have more than two levels, these two methods become entirely different! If A has$k$levels then (1|A:sub) estimates only one variance parameter, whereas (0+A|sub) estimates$k-1$variance parameters and$k(k-1)/2$correlations between them. My understanding is that classical ANOVA estimates only one variance parameter and so should be equivalent to the first method, not the second. Right? But for$k=2\$ both methods estimate one parameter, and I am still not quite sure why they disagree. – amoeba Oct 12 '16 at 15:59
• Getting back to this issue... I noticed that for the two-factor case where two lmer calls produce identical anova() output, the random effect variances are nevertheless quite different: see VarCorr(mod1) and VarCorr(mod2). I don't quite understand why this happens; do you? For mod3 and mod4, one can see that four out of seven variances for mod3 are actually equal to zero (for mod4 all seven are non-zero); this "singularity" in mod3 is probably why the anova tables differ. Apart from that, how would you use your "preferred way" if a and b had more than two levels? – amoeba Jun 28 '17 at 11:17

Are a, b, c fixed or random effects? If they are fixed, your syntax will simply be

summary(aov(y~a*b*c+Error(subject), d))
Error: subject
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 11  13.47   1.225

Error: Within
Df Sum Sq Mean Sq F value  Pr(>F)
a          1   1.41   1.408   1.730 0.19235
b          1   3.12   3.118   3.829 0.05399 .
c          1   6.38   6.381   7.836 0.00647 **
a:b        1   1.57   1.571   1.929 0.16889
a:c        1   1.09   1.091   1.339 0.25072
b:c        1   1.41   1.413   1.735 0.19168
a:b:c      1   0.10   0.101   0.124 0.72518
Residuals 77  62.70   0.814

library(lmerTest)
anova(lmer(y ~ a*b*c+(1|subject), data=d))
Analysis of Variance Table of type 3  with  Satterthwaite
approximation for degrees of freedom
Sum Sq Mean Sq NumDF  DenDF F.value   Pr(>F)
a     1.4085  1.4085     1 76.991  1.7297 0.192349
b     3.1180  3.1180     1 76.991  3.8291 0.053995 .
c     6.3809  6.3809     1 76.991  7.8363 0.006469 **
a:b   1.5706  1.5706     1 76.991  1.9289 0.168888
a:c   1.0907  1.0907     1 76.991  1.3394 0.250716
b:c   1.4128  1.4128     1 76.991  1.7350 0.191680
a:b:c 0.1014  0.1014     1 76.991  0.1245 0.725183

• They are fixed effects. However, the ANOVA model you fit is not the model that seems to be the classical repeated measures ANOVA model, see e.g., here. See the Error strata in your and my case. – Henrik Oct 2 '14 at 20:02
• Actually how they do it is incorrect. If you have a fully crossed factorial repeated measures design (or randomized block factorial design), you should get only 1 error term, aside from subject, for all effects (i.e., Within). See Experimental Design: Procedures for Behavioral Sciences (2013) by Kirk, chapter 10 (p.458) or my post here – Masato Nakazawa Oct 2 '14 at 20:12
• Let's sidestep this question for the moment and assume that the model I fitted would be the correct model. How would you fit this using lmer? I will nevertheless get my copy of Kirk (2nd edition only) and see what it says. – Henrik Oct 2 '14 at 20:27
• I'm curious to know what you think of Kirk's chapter. I think the chapter number in the 2nd ed. is different. Meanwhile I'll try to fit different lmer models. The best way to check the model fit is to check their dfs using lmerTest because the K-R approximation is supposed to give you exact dfs and hence p-values. – Masato Nakazawa Oct 2 '14 at 20:39
• I have the 2nd edition of Kirk, where I believe the relevant discussion is on pp. 443-449, which discusses a two-way (not three-way) example. Expected mean squares, either assuming additivity of A and B or not, are given on p. 447. Assuming A and B are fixed and subjects/blocks are random, we can see from the expected mean squares listed by Kirk under "non-additive model" that tests of A, B, and AB each involve different error terms, namely, the relevant interactions with block/subject. The same principle extends to the present three-way example. – Jake Westfall Nov 4 '14 at 3:31