Section 3: Mixed-Effects Models for Factorial Designs
Keith Lohse, PhD, PStat
2026-07-30
Note that this page is in the process of being updated 2022-10-04
For these examples, we are going to work with a fictional data that that has two fully nested (i.e., between-subjects) and two fully crossed (i.e., within-subject) factors. We have two hypothetical groups of older adults (OA, >65 y) and younger adults (YA, <35 y). Within each of those groups, half of the participants were in the control group (C) and half were in the treatment group (T). Each participant was also measured on four different trials (1 to 4) in three different conditions (A, B, and C).
First, we will open the relevant R packages that we will use and then we will download the data from the web. (If you don’t have these packages installed, you will need install the packages before running the library() function.) Once the data are downloaded we can use the head() function to view the first ten rows of the data.
library(tidyverse); library(RCurl); library(ez); library(lme4); library(car); library(lmerTest)
DATA <- read.csv("https://raw.githubusercontent.com/keithlohse/mixed_effects_models/master/data_AGING_example.csv",
stringsAsFactors = TRUE)head(DATA, 10)For the examples we will walk through, we will ignore different variables at different times. For instance, we will want to average across trials to create a data set with one observation per condition (which we will call data_COND), and we will want to average across conditions to create a data set with one observation for each trial (which we will call data_TIME).
Averaging across trials, here are the first ten rows of data_COND.
data_COND <- DATA %>% group_by(subID, condition, age_group, group) %>% # Note we ignore trial to average across it
summarize(speed = mean(speed, na.rm=TRUE)) %>% # I have included na.rm=TRUE even though there are no missing data
arrange(age_group, subID, condition) # finally, we can resort our data by subIDs within groups## `summarise()` has grouped output by 'subID', 'condition', 'age_group'. You can
## override using the `.groups` argument.head(data_COND, 10)Averaging across conditions, here are the first ten rows of data_TIME.
data_TIME <- DATA %>% group_by(subID, time, age_group, group) %>% # Note we ignore condition to average across it
summarize(speed = mean(speed, na.rm=TRUE)) %>% # I have included na.rm=TRUE even though there are no missing data
arrange(age_group, subID, time) # finally, we can resort our data by subIDs within groups## `summarise()` has grouped output by 'subID', 'time', 'age_group'. You can
## override using the `.groups` argument.head(data_TIME, 10)With these data in place, we can now test all of our different models. We are going to explore the appropriate mixed-effects regression (MER) models for these different situations:
- MER as One-Way Repeated Measures ANOVA
- When we have a single crossed factor (i.e., one within-subjects factor) in our mixed-effects model.
- MER as Two-Way Repeated Measures ANOVA
- When we have multiple crossed factors (i.e., two within-subjects factors) in our mixed-effects model.
- MER as Mixed-Factorial ANOVA with a Single Within-Subjects
Factor
- When we have a single crossed factor (i.e., one within-subjects factor) and a single nested factor (i.e., between-subjects factors) in our mixed-effects model.
- MER as Mixed-Factorial ANOVA with Multiple Within-Subjects
Factors
- When we have multiple crossed factors (i.e., two within-subjects factor) and a single nested factor (i.e., between-subjects factors) in our mixed-effects model.
- MER as Three-Way Repeated-Measures ANOVA (No
Between-Subjects Factor)
- When we have three crossed factors (i.e., three within-subjects factors) and no nested, between-subjects factor at all in our mixed-effects model. This situation is a natural extension of Sections 1 and 2 above (which handled one and two crossed factors, respectively) and it is also exactly the situation that Lohse, Kozlowski, and Strube (2023; the companion paper to this repository) describe when they discuss adding a third within-subject factor to a fully factorial design. We work through a fully-verified example below, including what happens when the “textbook correct” random-effects structure runs into trouble.
1. MER as One-Way Repeated Measures ANOVA
(A Single Crossed Factor)
For this example, we will focus on only the effect of condition, so we will use the data_COND dataset to average across different trials. First, let’s plot the data to get a better sense of what the data look like.
ggplot(data_COND, aes(x = condition, y = speed)) +
geom_point(aes(fill=condition), pch=21, size=2,
position=position_jitter(w=0.2, h=0))+
geom_boxplot(aes(fill=condition), col="black",
alpha=0.4, width=0.5, outlier.shape = NA)+
scale_x_discrete(name = "Condition") +
scale_y_continuous(name = "Speed (m/s)") +
theme(axis.text=element_text(size=16, color="black"),
axis.title=element_text(size=16, face="bold"),
plot.title=element_text(size=16, face="bold", hjust=0.5),
panel.grid.minor = element_blank(),
strip.text = element_text(size=16, face="bold"),
legend.position = "none")
1.1. As an ANOVA…
To implement a simple one-way repeated measures ANOVA, we have a few options. We could directly code our ANOVA using the aov() function in R:
summary(aov(speed ~ condition + Error(subID/condition), data=data_COND))##
## Error: subID
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 39 2.751 0.07054
##
## Error: subID:condition
## Df Sum Sq Mean Sq F value Pr(>F)
## condition 2 0.2820 0.14098 24.46 5.67e-09 ***
## Residuals 78 0.4496 0.00576
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Or we can use the ezANOVA() function from the “ez” package.
ezANOVA(data = data_COND,
dv = .(speed),
wid = .(subID),
within = .(condition)
)## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 2 condition 2 78 24.46039 5.674295e-09 * 0.08095761
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 2 condition 0.8955986 0.1230706
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF] p[HF]<.05
## 2 condition 0.9054678 2.494269e-08 * 0.947018 1.300658e-08 *Although these two different functions present the results in slightly different ways, note that the f-values for the two different omnibus tests match. By either method, the f-observed for the main-effect of condition is F(2,78)=24.46, p<0.001.
If you are familiar with issues of contrast versus treatment coding, Type I versus Type III Sums of Squared Errors, and collinearity, then directly controlling your own ANOVA using base R is probably a safe bet.
If you are not familiar with those terms/issues, then the ezANOVA code is probably the best option for you. However, I would strongly encourage you to develop a more detailed understanding of general-linear models before jumping into mixed-effect models.
By default the aov() function provides a test of your statistical model using Type I Sums of Squared Errors, whereas the ezANOVA() function provides a test of your statistical model using Type III Sums of Squared Errors. In a completely orthogonal design where all factors are statistically independent of each other (i.e., there is no collinearity), the Type I and Type III Sums of Squared Errors will agree. By default, R also uses treatment coding (or “dummy” coding) for categorical factors rather than orthogonal contrast codes. The omnibus F-tests that we see in the ANOVA output will be the same regardless of the types of codes used, but if we dig into individual regression coefficients, it is important to remember how these variables were coded so that we can interpret them correctly.
1.2. Getting the same result with a mixed-effect model…
Because we have a single within-subject factor, we will need to add a random-effect of subject to account for individual differences between subjects. By partitioning the between-subjects variance out of our model, we can fairly test the effect of condition, because our residuals will now be independent of each other.
# First we will define our model
mod1 <- lmer(speed ~
# Fixed Effects:
condition +
# Random Effects:
(1|subID),
# Define the data:
data=data_COND, REML = TRUE)
# We can then get the ANOVA results for our model:
anova(mod1)Most critically, note that the F-value, F(2,78)=24.46, p<0.001 is identical in both the RM ANOVA and in the mixed-effects regression model.
If we want to delve deeper into our model, we can also use the summary() function to get more information about model fit statistics, parameter estimates, random-effects and residuals.
summary(mod1)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: speed ~ condition + (1 | subID)
## Data: data_COND
##
## REML criterion at convergence: -162.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.5766 -0.3736 -0.0165 0.3682 3.8369
##
## Random effects:
## Groups Name Variance Std.Dev.
## subID (Intercept) 0.021594 0.14695
## Residual 0.005764 0.07592
## Number of obs: 120, groups: subID, 40
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.84250 0.02615 52.09109 32.215 <2e-16 ***
## conditionB 0.11831 0.01698 78.00000 6.970 9e-10 ***
## conditionC 0.05050 0.01698 78.00000 2.975 0.0039 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cndtnB
## conditionB -0.325
## conditionC -0.325 0.5002. MER as Two-Way Repeated Measures ANOVA
(Multiple Crossed Factors)
For this example, we will be analyzing both Time and Condition, so we will use our full data frame, DATA. First, let’s plot the data to get a better sense of what our data look like.
DATA$time <- factor(DATA$time)
ggplot(DATA, aes(x = time, y = speed)) +
geom_point(aes(fill=condition), pch=21, size=2,
position=position_jitterdodge(dodge.width = 0.5, jitter.width = 0.1))+
geom_boxplot(aes(fill=condition), col="black",
alpha=0.4, width=0.5, outlier.shape = NA)+
labs(fill="Condition")+
scale_x_discrete(name = "Trial") +
scale_y_continuous(name = "Speed (m/s)", limits = c(0,3)) +
theme(axis.text=element_text(size=16, color="black"),
axis.title=element_text(size=16, face="bold"),
legend.text=element_text(size=16, color="black"),
legend.title=element_text(size=16, face="bold"),
plot.title=element_text(size=16, face="bold", hjust=0.5),
panel.grid.minor = element_blank(),
strip.text = element_text(size=16, face="bold"),
legend.position = "bottom")
2.1. As an ANOVA…
Before we run our models, we want to convert trial to a factor so that our model is treating time categorically rather than continuously. (In later modules, we will discuss how to mix continuous and categorical factors.) Additionally, remember that we are now using our larger data set DATA rather than our aggregated data set data_COND. To implement a two-way repeated measures ANOVA, we have the same options as before. We can directly code our ANOVA using the aov() function in R:
DATA$time <- factor(DATA$time)
summary(aov(speed ~ condition*time + Error(subID/(condition*time)), data=DATA))##
## Error: subID
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 39 11.01 0.2822
##
## Error: subID:condition
## Df Sum Sq Mean Sq F value Pr(>F)
## condition 2 1.128 0.5639 24.46 5.67e-09 ***
## Residuals 78 1.798 0.0231
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: subID:time
## Df Sum Sq Mean Sq F value Pr(>F)
## time 3 10.618 3.539 106 <2e-16 ***
## Residuals 117 3.909 0.033
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: subID:condition:time
## Df Sum Sq Mean Sq F value Pr(>F)
## condition:time 6 0.7851 0.13084 13.53 3.59e-13 ***
## Residuals 234 2.2626 0.00967
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Or we can use the ezANOVA() function from the “ez” package.
ezANOVA(data = DATA,
dv = .(speed),
wid = .(subID),
within = .(time, condition)
)## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 2 time 3 117 105.95095 3.295544e-33 * 0.35881815
## 3 condition 2 78 24.46039 5.674295e-09 * 0.05610418
## 4 time:condition 6 234 13.53149 3.588475e-13 * 0.03973035
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 2 time 0.3262678 5.356910e-08 *
## 3 condition 0.8955986 1.230706e-01
## 4 time:condition 0.2380896 9.170350e-05 *
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF]
## 2 time 0.6917461 9.817867e-24 * 0.7313027 5.956088e-25
## 3 condition 0.9054678 2.494269e-08 * 0.9470180 1.300658e-08
## 4 time:condition 0.6852099 1.068940e-09 * 0.7760402 1.056951e-10
## p[HF]<.05
## 2 *
## 3 *
## 4 *Again, regardless of the coding approach that you use, these two different functions produce the same f-observed for the main-effects of time, F(3,117)=105.95, p<0.001, and Condition, F(2,78)=24.46, p<0.001, and the Time x Condition interaction, F(6,234)=13.53, p<0.001.
2.2. Two-way RM ANOVA as a mixed-effect model…
Because we now have two crossed factors, we need to account not only for the fact that we have multiple observations for each participant, but we have multiple observations for each person at each level of each factor.
For instance, for the effect of Condition, we actually have Condition effects at Trial 1, Trial 2, Trial 3, and Trial 4. The converse is also true for the effect of Time, we have four different trials in Condition A, Condition B, and Condition C. Thus, in order to appropriately account for the statistical dependencies in our data, we need to add random-effects of “time:subID” and “condition:subID” to the model.
The colon operator (“:”) means that we are crossing or multiplying these factors. That is, if we have subject ID’s A, B, and C and Trials 1, 2, 3, and 4, then we end up with A1, A2, A3, A4, B1, B2, B3, B4, etc.
For more information on how random-effects are specified and what they mean in R, I recommend looking at this discussion on Stack Exchange: https://stats.stackexchange.com/questions/228800/crossed-vs-nested-random-effects-how-do-they-differ-and-how-are-they-specified
# First we will define our model
mod2 <- lmer(speed ~
# Fixed Effects:
time*condition +
# Random Effects
(1|subID)+ (1|time:subID) + (1|condition:subID),
# Define your data,
data=DATA, REML=TRUE)
# We can then get the ANOVA results for our model:
anova(mod2)Again, for our purposes, the critical thing to note is that the F-values for the main-effects and interactions are the same between our RM ANOVAs and the mixed-effects model. Without delving into the mathematical details, this is a good demonstration that the appropriate random-effects our regression model make it analogous to the factorial ANOVA. This allows us to capitalize on the benefits of mixed-effects regression for designs that we would normally analyze using factorial ANOVA.
If we want to delve deeper into our model, we can also use the summary() function to get more information about model fit statistics, parameter estimates, random-effects and residuals.
summary(mod2)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: speed ~ time * condition + (1 | subID) + (1 | time:subID) + (1 |
## condition:subID)
## Data: DATA
##
## REML criterion at convergence: -454.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9509 -0.4451 -0.0673 0.3844 3.2827
##
## Random effects:
## Groups Name Variance Std.Dev.
## time:subID (Intercept) 0.007912 0.08895
## condition:subID (Intercept) 0.003346 0.05785
## subID (Intercept) 0.019616 0.14006
## Residual 0.009669 0.09833
## Number of obs: 480, groups: time:subID, 160; condition:subID, 120; subID, 40
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.00975 0.03184 109.12304 31.716 < 2e-16 ***
## time2 -0.15300 0.02965 249.81033 -5.160 5.03e-07 ***
## time3 -0.25350 0.02965 249.81033 -8.550 1.25e-15 ***
## time4 -0.26250 0.02965 249.81033 -8.853 < 2e-16 ***
## conditionB 0.24625 0.02551 260.37316 9.653 < 2e-16 ***
## conditionC 0.10125 0.02551 260.37316 3.969 9.34e-05 ***
## time2:conditionB -0.07025 0.03110 234.00008 -2.259 0.02479 *
## time3:conditionB -0.19800 0.03110 234.00008 -6.367 1.01e-09 ***
## time4:conditionB -0.24350 0.03110 234.00008 -7.831 1.68e-13 ***
## time2:conditionC -0.03825 0.03110 234.00008 -1.230 0.21990
## time3:conditionC -0.06425 0.03110 234.00008 -2.066 0.03991 *
## time4:conditionC -0.10050 0.03110 234.00008 -3.232 0.00141 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) time2 time3 time4 cndtnB cndtnC tm2:cB tm3:cB tm4:cB
## time2 -0.466
## time3 -0.466 0.500
## time4 -0.466 0.500 0.500
## conditionB -0.401 0.320 0.320 0.320
## conditionC -0.401 0.320 0.320 0.320 0.500
## tim2:cndtnB 0.244 -0.524 -0.262 -0.262 -0.609 -0.305
## tim3:cndtnB 0.244 -0.262 -0.524 -0.262 -0.609 -0.305 0.500
## tim4:cndtnB 0.244 -0.262 -0.262 -0.524 -0.609 -0.305 0.500 0.500
## tim2:cndtnC 0.244 -0.524 -0.262 -0.262 -0.305 -0.609 0.500 0.250 0.250
## tim3:cndtnC 0.244 -0.262 -0.524 -0.262 -0.305 -0.609 0.250 0.500 0.250
## tim4:cndtnC 0.244 -0.262 -0.262 -0.524 -0.305 -0.609 0.250 0.250 0.500
## tm2:cC tm3:cC
## time2
## time3
## time4
## conditionB
## conditionC
## tim2:cndtnB
## tim3:cndtnB
## tim4:cndtnB
## tim2:cndtnC
## tim3:cndtnC 0.500
## tim4:cndtnC 0.500 0.5003. MER as Mixed-Factorial ANOVA with One Crossed Factor
2 (Age) x 3 (Condition) ANOVA Example
For this example, we will now consider both the between-subject factor of age, as well as the within-subject factor of condition. Do to this, we will average over the different trials (1, 2, 3, and 4) to get one observation in each condition. This design is sometimes referred as a “mixed-factorial” design, because we have a mix of between-subjects and within-subject factors. Depending on your background, you might be more familiar with this as a “split plot” design. Split plot comes from agronomy (where a lot of statistics where developed!) and refers to the fact that some plots are assigned to different levels of Factor A (e.g., Treatment versus Control), but within each plot there is a second level of randomization to different levels of Factor B (e.g., Conditions A, B, or C). Both terms describe the same thing but differ in their unit of analysis. In psychology, a single person is usually the experimental unit (hence “within-subject” variables) whereas in agronomy or biology a physical area might be the unit of analysis (hence “split-plot” variables). The more general way to talk about these terms is as nested- versus crossed-factors.
- For fully nested factors, an experimental unit is represented in only one of the factors (e.g., a person can only be in the treatment group or the control group).
- For fully crossed factors, an experimental unit is represented at all levels of the factor (e.g., e.g., a person is tested in conditions A, B, and C).
In our example, participants are nested within age-group, because each person is represented at only one level of that factor (i.e., you are only a younger adult or older adult). However, Condition is a crossed factor, because each person is represented at all levels of Condition (i.e., each person was measured in all three conditions). We can see this more clearly if we plot all of our data.
ggplot(data_COND, aes(x = condition, y = speed)) +
geom_point(aes(fill=age_group), pch=21, size=2,
position=position_jitterdodge(dodge.width = 0.5, jitter.width = 0.1))+
geom_boxplot(aes(fill=age_group), col="black",
alpha=0.4, width=0.5, outlier.shape=NA)+
labs(fill="Age Group")+
scale_x_discrete(name = "Condition") +
scale_y_continuous(name = "Speed (m/s)") +
theme(axis.text=element_text(size=16, color="black"),
axis.title=element_text(size=16, face="bold"),
legend.text=element_text(size=16, color="black"),
legend.title=element_text(size=16, face="bold"),
plot.title=element_text(size=16, face="bold", hjust=0.5),
panel.grid.minor = element_blank(),
strip.text = element_text(size=16, face="bold"),
legend.position = "bottom")
3.1. As an ANOVA…
For this mixed factorial ANOVA, we have one factor of Condition that varies within-subjects, but we also have a factor of Age Group that varies between subjects. As before, we can directly code this into our analysis of variance using the aov() function or using the ezANOVA() function from the “ez” package.
summary(aov(speed ~ age_group*condition + Error(subID/condition), data=data_COND))##
## Error: subID
## Df Sum Sq Mean Sq F value Pr(>F)
## age_group 1 1.524 1.5238 47.18 3.71e-08 ***
## Residuals 38 1.227 0.0323
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: subID:condition
## Df Sum Sq Mean Sq F value Pr(>F)
## condition 2 0.2820 0.14098 27.271 1.18e-09 ***
## age_group:condition 2 0.0567 0.02834 5.482 0.00597 **
## Residuals 76 0.3929 0.00517
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Or we can use the ezANOVA() function from the “ez” package.
ezANOVA(data = data_COND,
dv = .(speed),
wid = .(subID),
within = .(condition),
between = .(age_group)
)## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 2 age_group 1 38 47.176031 3.706667e-08 * 0.48465601
## 3 condition 2 76 27.271263 1.181245e-09 * 0.14822123
## 4 age_group:condition 2 76 5.481689 5.972183e-03 * 0.03379572
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 3 condition 0.9008222 0.1448182
## 4 age_group:condition 0.9008222 0.1448182
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF]
## 3 condition 0.9097709 5.566488e-09 * 0.9530301 2.646471e-09
## 4 age_group:condition 0.9097709 7.698688e-03 * 0.9530301 6.815470e-03
## p[HF]<.05
## 3 *
## 4 *The ezANOVA() function provides us with more detail, for instance automatically providing Mauchly’s Test for Sphericity and both the Greenhouse-Geisser and Hyunh-Feldt corrections in the event that the sphericity assumption is violated. Most importantly, however, the outputs of these two functions agree when we look at the f-statistics for all main-effects and interactions when sphericity is assumed.
3.2. Mixed Factorial ANOVA as a mixed-effect model…
For this mixed-factorial design, we need to account for the fact that we have multiple observations coming from each person, so we will add a random-effect of “subID”. After accounting for this statistical dependence in our data, we can now fairly test the effects of Age Group, and Condition with residuals that are independent of each other.
# First we will define our model
mod3 <- lmer(speed ~
# Fixed Effects:
age_group*condition +
# Random Effects
(1|subID),
# Define your data,
data=data_COND, REML=TRUE)
# We can then get the ANOVA results for our model:
anova(mod3)Critically, note that the f-statistics and degrees of freedom all match our ANOVA outputs (when sphericity is assumed), with a main-effect of Age Group, F(1,38)=47.18, Condition, F(2,76)=27.27, and the Age Group x Condition interaction, F(2,76)=5.48.
If we want to delve deeper into our model, we can also use the summary() function to get more information about model fit statistics, parameter estimates, random-effects, and residuals.
summary(mod3)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: speed ~ age_group * condition + (1 | subID)
## Data: data_COND
##
## REML criterion at convergence: -189.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6329 -0.4286 -0.0582 0.3525 3.7641
##
## Random effects:
## Groups Name Variance Std.Dev.
## subID (Intercept) 0.009044 0.0951
## Residual 0.005169 0.0719
## Number of obs: 120, groups: subID, 40
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.70000 0.02666 62.99267 26.258 < 2e-16 ***
## age_groupYA 0.28500 0.03770 62.99267 7.560 2.11e-10 ***
## conditionB 0.16950 0.02274 76.00000 7.455 1.21e-10 ***
## conditionC 0.08875 0.02274 76.00000 3.903 0.000204 ***
## age_groupYA:conditionB -0.10238 0.03215 76.00000 -3.184 0.002107 **
## age_groupYA:conditionC -0.07650 0.03215 76.00000 -2.379 0.019862 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) ag_gYA cndtnB cndtnC a_YA:B
## age_groupYA -0.707
## conditionB -0.426 0.302
## conditionC -0.426 0.302 0.500
## ag_grpYA:cB 0.302 -0.426 -0.707 -0.354
## ag_grpYA:cC 0.302 -0.426 -0.354 -0.707 0.5004. MER as Mixed-Factorial ANOVA with Multiple Crossed Factors
2 (Age Group) x 3 (Condition) x 4 (Trial) ANOVA
Next, we will consider the more complicated example of our fully factorial design with a nested factor of Age-Group, and crossed factors of Condition and Time.
DATA$time <- factor(DATA$time)
ggplot(DATA, aes(x = time, y = speed)) +
geom_point(aes(fill=age_group), size=2, shape=21,
position=position_jitterdodge(dodge.width = 0.5, jitter.width = 0.1))+
geom_boxplot(aes(fill=age_group), col="black",
alpha=0.4, width=0.5, outlier.shape = NA)+
labs(fill="Age Group")+
facet_wrap(~condition)+
scale_x_discrete(name = "Trial") +
scale_y_continuous(name = "Speed (m/s)", limits = c(0,3)) +
theme(axis.text=element_text(size=16, color="black"),
axis.title=element_text(size=16, face="bold"),
legend.text=element_text(size=16, color="black"),
legend.title=element_text(size=16, face="bold"),
plot.title=element_text(size=16, face="bold", hjust=0.5),
panel.grid.minor = element_blank(),
strip.text = element_text(size=16, face="bold"),
legend.position = "bottom")
4.1. As an ANOVA…
For this mixed factorial ANOVA, we have two factors that vary within subjects, Condition and Time, and we have one factor that varies between subjects, Age Group. As before, we can directly code this into our analysis of variance using the aov() function or using the ezANOVA() function from the “ez” package.
summary(aov(speed ~ age_group*condition*time + Error(subID/(condition*time)), data=DATA))##
## Error: subID
## Df Sum Sq Mean Sq F value Pr(>F)
## age_group 1 6.095 6.095 47.18 3.71e-08 ***
## Residuals 38 4.910 0.129
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: subID:condition
## Df Sum Sq Mean Sq F value Pr(>F)
## condition 2 1.1278 0.5639 27.271 1.18e-09 ***
## age_group:condition 2 0.2267 0.1133 5.482 0.00597 **
## Residuals 76 1.5715 0.0207
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: subID:time
## Df Sum Sq Mean Sq F value Pr(>F)
## time 3 10.618 3.539 113.88 <2e-16 ***
## age_group:time 3 0.365 0.122 3.92 0.0105 *
## Residuals 114 3.543 0.031
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: subID:condition:time
## Df Sum Sq Mean Sq F value Pr(>F)
## condition:time 6 0.7851 0.13084 13.718 2.71e-13 ***
## age_group:condition:time 6 0.0880 0.01467 1.539 0.166
## Residuals 228 2.1746 0.00954
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ezANOVA(data = DATA,
dv = .(speed),
wid = .(subID),
within = .(condition, time),
between = .(age_group)
)## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 2 age_group 1 38 47.176031 3.706667e-08 * 0.33318007
## 3 condition 2 76 27.271263 1.181245e-09 * 0.08462822
## 5 time 3 114 113.882930 3.746328e-34 * 0.46536693
## 4 age_group:condition 2 76 5.481689 5.972183e-03 * 0.01824443
## 6 age_group:time 3 114 3.919720 1.049688e-02 * 0.02908814
## 7 condition:time 6 228 13.718345 2.705162e-13 * 0.06046299
## 8 age_group:condition:time 6 228 1.538560 1.664920e-01 0.00716581
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 3 condition 0.9008222 1.448182e-01
## 4 age_group:condition 0.9008222 1.448182e-01
## 5 time 0.3554491 3.872111e-07 *
## 6 age_group:time 0.3554491 3.872111e-07 *
## 7 condition:time 0.2268549 8.327332e-05 *
## 8 age_group:condition:time 0.2268549 8.327332e-05 *
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe
## 3 condition 0.9097709 5.566488e-09 * 0.9530301
## 4 age_group:condition 0.9097709 7.698688e-03 * 0.9530301
## 5 time 0.7021534 1.020427e-24 * 0.7443488
## 6 age_group:time 0.7021534 2.194748e-02 * 0.7443488
## 7 condition:time 0.6718683 1.241027e-09 * 0.7615676
## 8 age_group:condition:time 0.6718683 1.933581e-01 0.7615676
## p[HF] p[HF]<.05
## 3 2.646471e-09 *
## 4 6.815470e-03 *
## 5 4.686178e-26 *
## 6 1.975578e-02 *
## 7 1.230714e-10 *
## 8 1.856659e-014.2. Multi-Way Mixed Factorial ANOVA as a mixed-effect model…
For this mixed-factorial design, we need to account for the fact that we have multiple observations coming from each person, but we also need to acocunt for the fact that we multiple observations for each Condition and on each Trial. In order to account for this dependence in our data, we need to include random-effects of subject, subject:condition, and subject:time. Adding these random-effects to our model will make our mixed-effects model statistically equivalent to the mixed-factorial ANOVAs that we ran above.
# First we will define our model
mod4 <- lmer(speed ~
# Fixed Effects:
age_group*condition*time +
# Random Effects
(1|subID)+(1|condition:subID)+(1|time:subID),
# Define your data,
data=DATA, REML=TRUE)
# We can then get the ANOVA results for our model:
anova(mod4)Note that the f-statistics and the degrees of freedom match what we got from the ANOVA tables above (when sphericity is assumed). There are a lot of effects so I wont list them all, but note the main-effects of Age Group, F(1,38)=47.18, Condition, F(2,76)=27.27, and Time, F(3,114)=113.88.
Readers familiar with repeated measures ANOVA wont be surprised to see that that the output of the aov() function partitions the effects into different sections. The between-subject effect of Age Group is evaluated using residual error at the level of “Error:subID”. The effects of Condition and Age Group x Condition are evaluated using residuals of “Error: subID:condition”. This is precisely what we are doing when we include the three different random-intercepts in the mixed-effects regression. The \((1|subID)\) partitions out the average individual differences between subjects, \((1|condition:subID)\) and \((1|time:subID)\) partitions out the variance at the different levels of our repeated measures.
If we want to delve deeper into our model, we can also use the summary() function to get more information about model fit statistics, parameter estimates, random-effects and residuals.
summary(mod4)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## speed ~ age_group * condition * time + (1 | subID) + (1 | condition:subID) +
## (1 | time:subID)
## Data: DATA
##
## REML criterion at convergence: -463.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.8077 -0.4203 -0.0405 0.3229 3.2902
##
## Random effects:
## Groups Name Variance Std.Dev.
## time:subID (Intercept) 0.007181 0.08474
## condition:subID (Intercept) 0.002785 0.05277
## subID (Intercept) 0.007249 0.08514
## Residual 0.009538 0.09766
## Number of obs: 480, groups: time:subID, 160; condition:subID, 120; subID, 40
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.79400 0.03657 186.55256 21.710 < 2e-16
## age_groupYA 0.43150 0.05172 186.55256 8.343 1.58e-14
## conditionB 0.34850 0.03510 263.60569 9.928 < 2e-16
## conditionC 0.17750 0.03510 263.60569 5.056 8.01e-07
## time2 -0.06700 0.04089 249.82467 -1.639 0.102553
## time3 -0.13950 0.04089 249.82467 -3.412 0.000753
## time4 -0.16950 0.04089 249.82467 -4.145 4.65e-05
## age_groupYA:conditionB -0.20450 0.04964 263.60569 -4.119 5.09e-05
## age_groupYA:conditionC -0.15250 0.04964 263.60569 -3.072 0.002350
## age_groupYA:time2 -0.17200 0.05782 249.82467 -2.975 0.003222
## age_groupYA:time3 -0.22800 0.05782 249.82467 -3.943 0.000105
## age_groupYA:time4 -0.18600 0.05782 249.82467 -3.217 0.001468
## conditionB:time2 -0.15450 0.04368 228.00014 -3.537 0.000489
## conditionC:time2 -0.08950 0.04368 228.00014 -2.049 0.041587
## conditionB:time3 -0.26350 0.04368 228.00014 -6.033 6.44e-09
## conditionC:time3 -0.12300 0.04368 228.00014 -2.816 0.005285
## conditionB:time4 -0.29800 0.04368 228.00014 -6.823 7.94e-11
## conditionC:time4 -0.14250 0.04368 228.00014 -3.263 0.001273
## age_groupYA:conditionB:time2 0.16850 0.06177 228.00014 2.728 0.006868
## age_groupYA:conditionC:time2 0.10250 0.06177 228.00014 1.659 0.098394
## age_groupYA:conditionB:time3 0.13100 0.06177 228.00014 2.121 0.035011
## age_groupYA:conditionC:time3 0.11750 0.06177 228.00014 1.902 0.058388
## age_groupYA:conditionB:time4 0.10900 0.06177 228.00014 1.765 0.078951
## age_groupYA:conditionC:time4 0.08400 0.06177 228.00014 1.360 0.175185
##
## (Intercept) ***
## age_groupYA ***
## conditionB ***
## conditionC ***
## time2
## time3 ***
## time4 ***
## age_groupYA:conditionB ***
## age_groupYA:conditionC **
## age_groupYA:time2 **
## age_groupYA:time3 ***
## age_groupYA:time4 **
## conditionB:time2 ***
## conditionC:time2 *
## conditionB:time3 ***
## conditionC:time3 **
## conditionB:time4 ***
## conditionC:time4 **
## age_groupYA:conditionB:time2 **
## age_groupYA:conditionC:time2 .
## age_groupYA:conditionB:time3 *
## age_groupYA:conditionC:time3 .
## age_groupYA:conditionB:time4 .
## age_groupYA:conditionC:time4
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 24 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it5. MER as Three-Way Repeated-Measures ANOVA
(Three Crossed Factors, No Between-Subjects Factor)
In Sections 1 and 2 above, we built up the mixed-effects analogue of a repeated-measures ANOVA one crossed factor at a time (one-way, then two-way). In Sections 3 and 4, we then added a between-subjects factor (Age Group) alongside one or two within-subjects factors. But we have not yet pushed the purely within-subjects case (no between-subjects factor at all) beyond two crossed factors. This is worth doing for its own sake, because a fully within-subjects, three-way design is common in rehabilitation and movement science (e.g., every participant experiences every combination of Load x Terrain x Fatigue, or every combination of Medication x Time-of-Day x Task), and because it is exactly the situation that our companion paper (Lohse, Kozlowski, & Strube, 2023, Communications in Kinesiology) flags as getting “more complicated” once a third within-subject factor is added, but only sketches in a single sentence and footnote. Below, we build out a fully worked, numerically-verified example, including a very important complication that shows up in real (not just hypothetical) data: what happens when the “textbook correct” random-effects structure runs into estimation trouble.
For this example, we simulate a new hypothetical data set
(data_GAIT_3way_example.csv) in which 32 participants
walked an obstacle course under every combination of three fully crossed
within-subject factors:
- Load: carrying a Light vs. a Heavy pack.
- Terrain: walking on Flat ground, Uneven ground, or over Stairs.
- Fatigue: tested Pre- vs. Post-fatigue protocol.
Every participant is measured exactly once under all \(2 \times 3 \times 2 = 12\) combinations of
these factors (384 observations total), and there is no between-subjects
factor at all (e.g., no separate age or treatment groups). The R script
used to simulate this data set (simulate_data_GAIT_3way.R)
is included in the repository alongside the data, so that instructors
can regenerate it, inspect the true, known effects that were simulated,
or adapt the script to build their own examples with different numbers
of factors, levels, or subjects.
ggplot(DAT3, aes(x = terrain, y = speed, fill = load)) +
geom_boxplot(col="black", alpha=0.6, width=0.6, outlier.shape = NA,
position = position_dodge(width=0.7)) +
geom_point(pch=21, size=1.5, alpha=0.6,
position = position_jitterdodge(dodge.width = 0.7, jitter.width = 0.1)) +
facet_wrap(~fatigue) +
labs(fill="Load") +
scale_x_discrete(name = "Terrain") +
scale_y_continuous(name = "Speed (m/s)") +
theme(axis.text=element_text(size=14, color="black"),
axis.title=element_text(size=14, face="bold"),
legend.text=element_text(size=14, color="black"),
legend.title=element_text(size=14, face="bold"),
strip.text = element_text(size=14, face="bold"),
panel.grid.minor = element_blank(),
legend.position = "bottom")
5.1. As an ANOVA…
Because every factor here is within-subjects, we account for
the statistical dependence of every observation coming from the same
person by nesting the entire factorial design under subject in
our Error() term:
aov_3way <- aov(speed ~ load*terrain*fatigue + Error(subID/(load*terrain*fatigue)), data=DAT3)
summary(aov_3way)##
## Error: subID
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 31 2.572 0.08298
##
## Error: subID:load
## Df Sum Sq Mean Sq F value Pr(>F)
## load 1 3.778 3.778 375.5 <2e-16 ***
## Residuals 31 0.312 0.010
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: subID:terrain
## Df Sum Sq Mean Sq F value Pr(>F)
## terrain 2 6.619 3.310 207.1 <2e-16 ***
## Residuals 62 0.991 0.016
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: subID:fatigue
## Df Sum Sq Mean Sq F value Pr(>F)
## fatigue 1 1.3799 1.3799 122.8 2.62e-12 ***
## Residuals 31 0.3484 0.0112
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: subID:load:terrain
## Df Sum Sq Mean Sq F value Pr(>F)
## load:terrain 2 0.1650 0.08252 15.76 2.93e-06 ***
## Residuals 62 0.3247 0.00524
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: subID:load:fatigue
## Df Sum Sq Mean Sq F value Pr(>F)
## load:fatigue 1 0.1155 0.11554 14.89 0.000539 ***
## Residuals 31 0.2405 0.00776
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: subID:terrain:fatigue
## Df Sum Sq Mean Sq F value Pr(>F)
## terrain:fatigue 2 0.0018 0.000879 0.145 0.866
## Residuals 62 0.3766 0.006074
##
## Error: subID:load:terrain:fatigue
## Df Sum Sq Mean Sq F value Pr(>F)
## load:terrain:fatigue 2 0.0014 0.000701 0.213 0.809
## Residuals 62 0.2039 0.003289Notice that this single Error() call produces
eight different strata: the (unused, since we have no
between-subjects factor) Error: subID stratum, one stratum
for each of the three main effects (subID:load,
subID:terrain, subID:fatigue), one stratum for
each of the three two-way interactions (subID:load:terrain,
subID:load:fatigue, subID:terrain:fatigue),
and finally the three-way interaction stratum
(subID:load:terrain:fatigue). Because each subject
contributes exactly one observation to every cell of the design, this
final three-way stratum is also serving as our “leftover” residual error
– there is no deeper partition of the data left to make.
5.2. Three-Way RM ANOVA as a mixed-effect model…
To reproduce these eight strata in a single lmer() call,
we can extend the same logic we used for the two-way, purely
within-subjects case in Section 2 above. There, to account for two
crossed factors we needed random intercepts for subject,
subject:A, and subject:B, and the model’s own
residual automatically played the role of the finest
(subject:A:B) stratum. For three crossed
factors, we push that logic one level deeper: we need a random intercept
for subject, one for each main effect nested in
subject (subject:load,
subject:terrain, subject:fatigue), and one for
each two-way interaction nested in subject
(subject:load:terrain, subject:load:fatigue,
subject:terrain:fatigue). The model’s residual then
automatically absorbs the finest, three-way stratum
(subject:load:terrain:fatigue), which is exactly the
stratum used to test the three-way interaction in the classical ANOVA
above.
# First we will define our model
mod5 <- lmer(speed ~
# Fixed Effects:
load*terrain*fatigue +
# Random Effects: one term for subject, one for each factor
# nested in subject, and one for each two-way interaction
# nested in subject:
(1|subID) +
(1|subID:load) + (1|subID:terrain) + (1|subID:fatigue) +
(1|subID:load:terrain) + (1|subID:load:fatigue) + (1|subID:terrain:fatigue),
data=DAT3, REML=TRUE)
# We can then get the ANOVA results for our model:
anova(mod5)# It is always worth checking whether the model actually converged cleanly:
isSingular(mod5)## [1] FALSEWith isSingular(mod5) returning FALSE, this
model converged without hitting a boundary, and – as promised – every
F-value and every denominator degrees of freedom matches the classical
repeated-measures ANOVA above exactly: Load, F(1,31)=375.52; Terrain,
F(2,62)=207.06; Fatigue, F(1,31)=122.77; Load x Terrain, F(2,62)=15.76;
Load x Fatigue, F(1,31)=14.89; Terrain x Fatigue, F(2,62)=0.14; and Load
x Terrain x Fatigue, F(2,62)=0.21. As in Sections 1, 2, and 4, this is a
good demonstration that the appropriately-specified mixed-effects model
is statistically equivalent to the classical repeated-measures ANOVA for
a complete, balanced, fully-crossed design.
A general recipe. Generalizing across sections: for
\(k\) fully crossed, purely
within-subjects factors with no missing cells, you need a random
intercept of subject, plus one random intercept for
subject crossed with every combination of 1 up to \(k-1\) of your within-subjects factors. The
model’s residual then automatically plays the role of the finest, \(k\)-way stratum. For \(k=1\): (1|subject). For \(k=2\):
(1|subject) + (1|subject:A) + (1|subject:B). For \(k=3\): the seven-term structure used in
mod5 above. Notice how quickly this grows – a fourth
crossed factor would require \(1 + 4 + 6 =
11\) random-effect terms before the residual even gets involved!
This rapid growth in complexity is itself an important, practical reason
to think carefully before adding more and more crossed within-subject
factors to a single design (see Section 5.3 below).
5.3. A Note on Boundary (Singular) Fits in Saturated, Three-Way Designs
The random-effects structure above is the theoretically “correct”
one, in the sense that it reproduces the classical ANOVA’s error strata
exactly. But it asks REML to estimate seven separate variance
components from data that, ultimately, come from only 32 subjects – and
some of those variance components (like
subject:load:fatigue, which is the interaction of two
factors that each only have two levels) reflect only a
single contrast per subject. This is a very easy structure to over-ask
of your data, and it is worth seeing exactly what happens when you
do.
To illustrate, let’s refit the identical model, but using only the first 12 (of 32) subjects from the very same simulated population:
DAT3_small <- DAT3 %>% filter(subID %in% levels(subID)[1:12]) %>% droplevels()
mod5_small <- lmer(speed ~
load*terrain*fatigue +
(1|subID) +
(1|subID:load) + (1|subID:terrain) + (1|subID:fatigue) +
(1|subID:load:terrain) + (1|subID:load:fatigue) + (1|subID:terrain:fatigue),
data=DAT3_small, REML=TRUE)## boundary (singular) fit: see help('isSingular')isSingular(mod5_small)## [1] TRUEVarCorr(mod5_small)## Groups Name Std.Dev.
## subID:terrain:fatigue (Intercept) 4.6103e-02
## subID:load:terrain (Intercept) 2.8232e-02
## subID:load:fatigue (Intercept) 3.2771e-02
## subID:terrain (Intercept) 2.8384e-02
## subID:fatigue (Intercept) 1.6710e-02
## subID:load (Intercept) 7.0654e-07
## subID (Intercept) 8.8057e-02
## Residual 6.1938e-02anova(mod5_small)With only 12 subjects, R warns us of a “boundary (singular) fit,” and
looking at VarCorr(mod5_small) shows why: the variance
estimate for subID:load has been pushed all the way down to
0.00. Load only has two levels, so subject:load reflects a
single per-subject contrast – exactly the kind of variance component
that is hardest to distinguish from residual noise with a modest sample
size.
Critically, this is not just a cosmetic warning. Because
subID:load collapsed into the residual, every fixed
effect that depends on that stratum now uses a different,
non-integer denominator degrees of freedom, and no longer matches the
classical ANOVA exactly: Load itself (F=160.32, df=1,21.1, versus the
ANOVA’s F=188.2, df=1,11), Load x Terrain (F=3.08, df=2,26.1, versus the
ANOVA’s F=2.94, df=2,22 – note this even flips from non-significant to a
different p-value), and Load x Fatigue (F=8.93, df=1,17.9, versus the
ANOVA’s F=7.97, df=1,11) have all drifted noticeably from their ANOVA
counterparts. Terrain, Fatigue, and Terrain x Fatigue – none of which
depended on the collapsed subID:load term – still match
almost exactly (e.g., Terrain: F=108.75 vs. the ANOVA’s 108.7, both with
df=2,22). The three-way interaction shows only a small drift (F=0.40,
df=2,21.6, versus the ANOVA’s F=0.41, df=2,22), since it is one step
further removed from the collapsed term. In other words, a singular fit
does not just add a scary warning message; in a design like this one, it
can measurably erode the very ANOVA-equivalence property that was the
whole reason for building this random-effects structure in the first
place, and the terms most directly tied to the collapsed variance
component are affected the most.
This is precisely the tension that Barr, Levy, Scheepers, & Tily (2013) and Bates, Kliegl, Vasishth, & Baayen (2015) debate under the labels “keep it maximal” versus “parsimonious mixed models” (the latter, along with a closely related paper by Barr alone, is already in the reference list of the companion Communications in Kinesiology paper for this repository), and that Matuschek, Kliegl, Vasishth, Baayen, & Bates (2017) resolve empirically: keeping every theoretically-justified random effect can protect against inflated Type I error, but it can also cost real statistical power once a variance component is not well-supported by the data, and (as shown here) it does not even guarantee the exact ANOVA-equivalence it was meant to provide. In practice, we would suggest:
- If your design is complete and balanced (as it is
here – every subject has exactly one observation in every cell, with no
missing data), the classical repeated-measures ANOVA (
aov()orezANOVA()) will always give you the “gold standard” answer directly, without any risk of a singular fit. The main reasons to reach for the mixed-effects version of a fully-crossed, purely within-subjects design are (a) you do have some missing data or an unbalanced design, where the classical ANOVA’sError()syntax breaks down butlmer()still works, or (b) you specifically want to report and interpret the variance components themselves (e.g., “how much do people vary in their Load x Fatigue effect?”). - If you need the mixed-effects version, consider a
more parsimonious structure: dropping the correlations among random
effects with the
||operator, or dropping the specific variance components that are estimated at (or very near) the zero boundary, and refitting. This will not recover the exact ANOVA F-values, but it can resolve the convergence warning and typically improves statistical power relative to the maximal model. - Plan for this ahead of time. If you know your design will have three or more fully crossed within-subjects factors, simulating data first (as we did to build this section) is a cheap way to check whether your planned sample size can support a maximal random-effects structure at all, before you collect any real data.
References for this Section
Barr, D. J., Levy, R., Scheepers, C., & Tily, H. J. (2013). Random effects structure for confirmatory hypothesis testing: Keep it maximal. Journal of Memory and Language, 68(3), 255-278. https://doi.org/10.1016/j.jml.2012.11.001
Bates, D., Kliegl, R., Vasishth, S., & Baayen, H. (2015). Parsimonious Mixed Models (Version 2). arXiv. https://doi.org/10.48550/ARXIV.1506.04967
Lohse, K. R., Kozlowski, A. J., & Strube, M. J. (2023). Model Specification in Mixed-Effects Models: A Focus on Random Effects. Communications in Kinesiology. https://doi.org/10.51224/cik.2023.52
Matuschek, H., Kliegl, R., Vasishth, S., Baayen, H., & Bates, D. (2017). Balancing Type I error and power in linear mixed models. Journal of Memory and Language, 94, 305-315. https://doi.org/10.1016/j.jml.2017.01.001