Specify an lmer model for fully crossed data with replications

Question

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")

statmerkur · Accepted Answer · 2022-09-14T22:29:35.923

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. — sharoz, Sep 15 '22 at 16:57

score 2 · Answer 2 · answered Sep 08 '22 at 16:08

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. — sharoz, Sep 12 '22 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). — statmerkur, Sep 14 '22 at 23:01

Specify an lmer model for fully crossed data with replications

ANOVA with afex

Equivalent with lmer?

Edit:

2 Answers2

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

[...]

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

ANOVA with `afex`

Equivalent with `lmer`?