# Specify an lmer model for fully crossed data with replications

I have experiment data with multiple independent variables:

• shape - 2 factors
• color - 2 factors
• size - 3 factors

Each subject provides a response (the dependent variable) for each combination of those independent variables many times.

Here is example data:

library(tidyverse)

data = expand_grid(
subject = paste("S", 1:11),
shape = c("circle", "square"),
color = c("red", "blue"),
size = c("small", "medium", "large"),
replicate = 1:50 # many replicate responses
) %>%
mutate(response = rnorm(n()))


# ANOVA with afex

The afex library lets you describe the variables and provides the results of an F-test:

afex::aov_ez(
id = "subject",
dv = "response",
data = data,
within = c("shape", "color", "size"),
anova_table=list(correction = "none") # turn off sphericity correction
)

Response: response
Effect    df  MSE       F   ges p.value
1            shape 1, 10 0.01    0.09 <.001    .768
2            color 1, 10 0.01    1.43  .008    .260
3             size 2, 20 0.01 7.88 **  .054    .003
4      shape:color 1, 10 0.03    1.19  .015    .302
5       shape:size 2, 20 0.02    1.68  .030    .211
6       color:size 2, 20 0.02  3.32 +  .064    .057
7 shape:color:size 2, 20 0.03    2.31  .050    .125
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1


# Equivalent with lmer?

Is it possible to specify an equivalent model via lmer() and get a similar output? If so, how?

These attempts don't seem right, as the degrees of freedom seem way off:

lmer(response ~ shape*color*size + (1 | subject), data) %>%
report::report()

lmer(response ~ shape*color*size + (1 | subject) + (1 | subject:shape) + (1 | subject:color) + (1 | subject:size), data) %>%
report::report()


# Edit:

It might be easier to match the afex output by aggregating first

data_aggregated = data %>%
group_by(subject, shape, color, size) %>%
summarise(response = mean(response), .groups = "drop")


Each observation in the data you simulate is independent from all the other observations. This doesn't match your $$2 \times 2 \times 3$$ within-subjects design with $$50$$ observations per cell and subject. In what follows I'll ignore that and act as if data and design would match.

The repeated measures ANOVA you specify via afex::aov_ez() works with data that are mean-aggregated to $$2 \times 2 \times 3$$ values per subject, one for each cell of your design.

library(tidyverse)
library(lmerTest)
set.seed(42)

data = expand_grid(
subject = paste("S", 1:11),
shape = c("circle", "square"),
color = c("red", "blue"),
size = c("small", "medium", "large"),
replicate = 1:50 # many replicate responses
) %>%
mutate(response = rnorm(n()))

rm_anova <- afex::aov_ez(
id = "subject",
dv = "response",
data = data,
within = c("shape", "color", "size"),
anova_table = list(correction = "none") # turn off sphericity correction
)

summary(rm_anova$$data$$long)
# subject      shape     color        size       response
# S 1    :12   circle:66   red :66   small :44   Min.   :-0.398478
# S 10   :12   square:66   blue:66   medium:44   1st Qu.:-0.111595
# S 11   :12                         large :44   Median :-0.010960
# S 2    :12                                     Mean   :-0.007078
# S 3    :12                                     3rd Qu.: 0.092962
# S 4    :12                                     Max.   : 0.411217
# (Other):60


To mimic the afex::aov_ez()-ANOVA by a linear mixed model (which might not be the most sensible thing to do) you should therefore work with the aggregated data.

Let's first state the repeated measures ANOVA explicitly in syntax of stats::aov(), which shows the involved error strata:

summary(aov(response ~ shape * color * size + Error(subject / (shape * color * size)),
data = rm_anova$$data$$long))
# [...]
# Error: subject:shape
#           Df  Sum Sq  Mean Sq F value Pr(>F)
# shape      1 0.00279 0.002794   0.155  0.702
# Residuals 10 0.18025 0.018025
#
# Error: subject:color
#           Df  Sum Sq  Mean Sq F value Pr(>F)
# color      1 0.00835 0.008346   0.697  0.423
# Residuals 10 0.11976 0.011976
#
# Error: subject:size
#           Df Sum Sq  Mean Sq F value Pr(>F)
# size       2 0.0061 0.003049   0.109  0.897
# Residuals 20 0.5570 0.027849
#
# Error: subject:shape:color
#             Df  Sum Sq Mean Sq F value Pr(>F)
# shape:color  1 0.01731 0.01731   1.077  0.324
# Residuals   10 0.16068 0.01607
#
# Error: subject:shape:size
#             Df Sum Sq  Mean Sq F value Pr(>F)
# shape:size  2 0.0041 0.002049   0.118   0.89
# Residuals  20 0.3480 0.017399
#
# Error: subject:color:size
#             Df Sum Sq Mean Sq F value Pr(>F)
# color:size  2 0.0713 0.03563   1.887  0.178
# Residuals  20 0.3776 0.01888
#
# Error: subject:shape:color:size
#                  Df Sum Sq Mean Sq F value Pr(>F)
# shape:color:size  2 0.0642 0.03210   1.474  0.253
# Residuals        20 0.4355 0.02178


The corresponding linear mixed model in lme4::lmer()/lmerTest::lmer() syntax is

lmm <- lmer(response ~ shape * color * size +
(1 | subject) +
(1 | subject:shape) + (1 | subject:color) + (1 | subject:size) +
(1 | subject:shape:color) + (1 | subject:shape:size) + (1 | subject:color:size),
data = rm_anova$$data$$long)


However, the variance components in this linear mixed model are estimated by a different method (REML by default) and constrained to be non-negative. As a consequence, the corresponding ANOVA-table might differ from the one based on afex::aov_ez()/stats::aov(), see also this answer.

rm_anova
# [...]
#             Effect    df  MSE    F  ges p.value
# 1            shape 1, 10 0.02 0.16 .001    .702
# 2            color 1, 10 0.01 0.70 .004    .423
# 3             size 2, 20 0.03 0.11 .003    .897
# 4      shape:color 1, 10 0.02 1.08 .008    .324
# 5       shape:size 2, 20 0.02 0.12 .002    .890
# 6       color:size 2, 20 0.02 1.89 .031    .178
# 7 shape:color:size 2, 20 0.02 1.47 .028    .253
# [...]

anova(lmm, ddf = "Kenward-Roger")
# [...]
#                    Sum Sq  Mean Sq NumDF DenDF F value Pr(>F)
# shape            0.002794 0.002794     1    10  0.1550 0.7020
# color            0.008346 0.008346     1    10  0.4631 0.5116
# size             0.005130 0.002565     2    20  0.1424 0.8682
# shape:color      0.017306 0.017306     1    10  0.9603 0.3502
# shape:size       0.004098 0.002049     2    20  0.1137 0.8931
# color:size       0.071255 0.035628     2    20  1.9771 0.1646
# shape:color:size 0.064199 0.032099     2    20  1.7813 0.1941

• Thank you! Yeah, I tried to keep the example really simple, but simulating a DV that was actually a function of the variables would have, in hind sight, been more clear. Commented Sep 15, 2022 at 16:57

This can be done by using lmerTest, which is an extension for lmer. The correct formula is response ~ shape*color*size + (1 | subject).

Data setup:

library(tidyverse)
library(lme4)
library(lmerTest)

set.seed(12345)
data = expand_grid(
subject = paste("S", 1:11),
shape = c("circle", "square"),
color = c("red", "blue"),
size = c("small", "medium", "large"),
replicate = 1:50 # many replicate responses
) %>%
mutate(response = rnorm(n()))


afex::aov_ez

afex::aov_ez(
id = "subject",
dv = "response",
data = data,
within = c("shape", "color", "size"),
anova_table=list(correction = "none") # turn off sphericity correction
)

Anova Table (Type 3 tests)

Response: response
Effect    df  MSE      F   ges p.value
1            shape 1, 10 0.03   0.02 <.001    .888
2            color 1, 10 0.02 3.51 +  .030    .090
3             size 2, 20 0.02   0.30  .005    .747
4      shape:color 1, 10 0.01   0.38  .001    .551
5       shape:size 2, 20 0.03   0.52  .011    .603
6       color:size 2, 20 0.02   0.17  .002    .848
7 shape:color:size 2, 20 0.02   0.64  .009    .537
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1
Warning message:
More than one observation per cell, aggregating the data using mean (i.e, fun_aggregate = mean)!


lmer/lmerTest

fit <- lmer(response ~ shape*color*size + (1 | subject), data)
anova(fit)
> anova(fit)
Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)
shape            0.0312  0.0312     1  6578  0.0314 0.85928
color            3.8021  3.8021     1  6578  3.8336 0.05028 .
size             0.6400  0.3200     2  6578  0.3226 0.72424
shape:color      0.1783  0.1783     1  6578  0.1797 0.67161
shape:size       1.3370  0.6685     2  6578  0.6740 0.50969
color:size       0.2897  0.1448     2  6578  0.1460 0.86414
shape:color:size 1.1174  0.5587     2  6578  0.5633 0.56934
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1



In general, the p-value for aov_ez and lmer/lmerTest are pretty close. However, I am worry about MSE/F in the aov_ez and ddf in lmerTest. Therefore, I do a bit more investigation mentioned in this post:

Using nlme::lme

model <- nlme::lme(response ~ shape*color*size, random=~1|subject, data = data)
anova(model)

numDF denDF  F-value p-value
(Intercept)          1  6578 0.071599  0.7890
shape                1  6578 0.031435  0.8593
color                1  6578 3.833618  0.0503
size                 2  6578 0.322648  0.7242
shape:color          1  6578 0.179745  0.6716
shape:size           2  6578 0.674027  0.5097
color:size           2  6578 0.146024  0.8641
shape:color:size     2  6578 0.563332  0.5693


Using lm

fit.lm <- lm(response ~ shape*color*size +  subject, data)
anova(fit.lm)

Analysis of Variance Table

Response: response
Df Sum Sq Mean Sq F value  Pr(>F)
shape               1    0.0  0.0312  0.0314 0.85928
color               1    3.8  3.8021  3.8336 0.05028 .
size                2    0.6  0.3200  0.3226 0.72424
subject            10   10.4  1.0445  1.0531 0.39535
shape:color         1    0.2  0.1783  0.1797 0.67161
shape:size          2    1.3  0.6685  0.6740 0.50969
color:size          2    0.3  0.1448  0.1460 0.86414
shape:color:size    2    1.1  0.5587  0.5633 0.56934
Residuals        6578 6524.0  0.9918
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


It seems like lme and lm agrees with lme4/lmerTest. As shown in the results in lm, ddf 6578 comes from the degree of freedom for residuals. I believe this is correct since this is a balance design. Therefore, I would go with the results from lmerTest.

• Thanks. The different degrees of freedom suggest a very different analysis approach, possibly one that fails to account for the fully crossed variables. The post you linked to indeed uses both within and between-subject variables. Commented Sep 12, 2022 at 13:30
• A linear mixed model that only includes a single by-subject random intercept is not equivalent to the model the OP specifies via afex::aov_ez(), and is very unlikely to be appropriate for a $2 \times 2 \times 3$ within-subjects design with $50$ observations per cell and subject. This might be more obvious with simulated data that match the design (instead of $6600$ i.i.d draws from a standard normal distribution). Commented Sep 14, 2022 at 23:01