Welcome to Barely Informative Content (or BIC for short); this is a corner of the internet where I ramble about statistics and other topics that take my fancy. Expect niche takes and mild existential dread. Reader discretion is advised.
Date: 23rd September 2026.
If you're anything like me, you may have found that understanding the "maximal model" structure in mixed-effects modelling isn't always straightforward, especially once experimental designs start getting more complex. However, I have written this post to try and make it a little easier.
Maximal models refer to models where the random-effects structure incorporates all elements justified by the experimental design. A bit opaque, I know, but stick with me... we will work through some examples below!
The push for maximal random-effects structures largely comes from work by Barr et al. (2013), who showed that when things aren't quite maximal, you can inflate Type I error rates. In other words, when you don't specify the random-effects structure correctly by incorporating all of the right elements, you are more likely to find an effect when one shouldn't actually be there... **dramatic music**
Note: Maximal isn't necessarily always better. If you're really interested in understanding why, see Matuschek et al. (2017).
But if we want to try fitting the maximal random-effects structure justified by our experimental design, we first need to know what that structure actually is. So, how do we figure that out?
Let's say we have a study interested in how schema congruence impacts later memory. Participants complete two within-subject conditions: congruent (items appear in a context you would expect, such as a toilet in a bathroom) and incongruent (items appear somewhere unexpected, such as a microwave in the living room). Each item-location association is unique and belongs to one of these two conditions. We later test participants' memory for these associations to see if people remember more or less depending on whether the association is congruent or incongruent.
To assess this, we decide to average across all of the items and calculate two proportions per participant, one for the congruent items and one for the incongruent items. From there, we fit a linear mixed-effects model with the following formula:
Recall ~ Congruency + (1 | Participant)
What is our model saying?
First, Congruency is our fixed effect. We are asking whether, on average, memory differs between congruent and incongruent events. We then have (1 | Participant). This is our random intercept; it allows participants to differ in their overall memory ability. In other words, we are telling the model to account for the fact that some people might simply have better memory than others. Importantly, including Participant as a random effect also accounts for the fact that each participant provides multiple responses. Responses from the same person will be more alike than responses from completely different people, so we don't want to treat them as entirely independent observations!
However, the model above is NOT maximal!
Why? Well, we know from our experimental design that Congruency varies within Participant. What I mean is, every participant completes both congruent and incongruent trials. Because participants experience both levels of Congruency, it is possible that the effect differs between participants. Congruency might have a huge effect on memory for Participant 1, a small effect for Participant 2, and even the opposite effect for Participant 3. At present, our model doesn't allow for this... it simply allows participants to have different overall levels of memory, but assumes that the effect of Congruency is the same across participants.
So, because Congruency varies within Participant, we should allow the effect of Congruency to vary across Participants. In other words, we add a random slope for Congruency:
Recall ~ Congruency + (1 + Congruency | Participant)
Our model now allows participants to differ in two ways:
Participants can differ in their overall memory ability (our Participant random intercept).
Participants can differ in the size and direction of the Congruency effect (our Congruency random slope).
In the context of this experimental design and analytical approach, we have now fit the maximal random-effects structure.
The example above was a very simple one. In reality, we often have more complex experimental designs, which means we need to think more carefully about which variables vary within and between each of our random effects. So, let's expand our schema example by adding a between-subjects variable: Age Group (AgeGrp). Age Group has two levels: younger adults and older adults. Each participant belongs to only one of these groups. That is, a younger adult cannot simultaneously be an older adult [unless they have somehow figured out time travel...]! However, like before, each participant still completes both Congruency conditions.
Our model now looks like this:
Recall ~ AgeGrp * Congruency + (1 + Congruency | Participant)
Interestingly, this model is still maximal! You may be a little confused about how when we've added an entirely new variable? Well, it's because Age Group does not vary within Participant. A given participant is either a younger or an older adult; they cannot be both. So, unlike Congruency, there is no within-participant variation in Age Group. We cannot ask whether the effect of being younger versus older varies within an individual participant because no participant experiences both levels! Therefore, we can't add a random slope for Age Group within Participant.
What I am trying to emphasise here is that we shouldn't simply be asking: "What predictors are in my model?", but instead: "Which predictors vary within each of my random effects?"... Hold onto that question; it's going to become useful in just a second!
The models above assumed that we carried out a linear mixed-effects model. To do this, we averaged across items and analysed the resulting proportions. In reality, we would probably want to work with the trial-level data, where each trial is coded as either correct (1) or incorrect (0). But why would we want to do this? Well, averaging throws away information! Trial-level analysis retains every response and allows us to model the probability of a correct response directly. It also allows us to account for the fact that some participants simply have better memory and that some items may simply be easier to remember than others. Of course, because our outcome is now binary, we would use a generalised linear mixed-effects model rather than a linear mixed-effects model... but don't worry about this for now! Instead, let's focus on what our model formula will look like...
Recall ~ AgeGrp * Congruency + (1 + Congruency | Participant) + (1 | Item)
You should immediately notice that we now have an additional random intercept: (1 | Item). This works similarly to our random intercept for Participant. Just as some participants may have better memory than others, some items may simply be easier to remember than others. By including the random intercept for Item, we allow each item to have its own baseline level of memorability. It also accounts for the fact that the same items are seen by multiple participants. Responses involving the same item may therefore be more similar than responses involving a completely different item.
However, by moving to an item-level model, we no longer have a maximal model structure!
So, what would the maximal structure be now that we have Items in the mix too? Your initial instinct might be: "Well, Congruency is one of our predictors, we know it is within-subjects, so surely we also need a random slope for Congruency within Item?" Something like this:
Recall ~ AgeGrp * Congruency + (1 + Congruency | Participant) + (1 + Congruency | Item)
But remember what we explicitly stated at the beginning:
"Each item-location association is unique and belongs to one of these two conditions."
In other words, each Item is either congruent or incongruent. An Item does not appear in both conditions. This means that Congruency does not vary within Item. So, including (1 + Congruency | Item) doesn't make sense for this particular design. To estimate a random slope for Congruency within Item, we would need to observe the same Item under different levels of Congruency. For example, having an experiment where we use the bathroom context with both a congruent item (e.g., toilet) and an incongruent item (e.g., toaster). But that's not how we've designed our experiment! Each Item has only one level of Congruency.
Okay, so does that mean we did have the maximal model before...? Nope! There is another factor in our model that does vary within Item: Age Group.
Think about how the experiment is structured. The same Item is seen by both younger and older adults. So, although Age Group DOES NOT vary within Participant, it DOES vary within Item. For each Item, we have responses from younger adults and responses from older adults. This means the effect of Age Group could differ across Items. For some Items, younger and older adults might perform very similarly, while for other Items, there might be a much larger difference between the two groups! So, we can add a random slope for Age Group within Item:
Recall ~ AgeGrp * Congruency + (1 + Congruency | Participant) + (1 + AgeGrp | Item)
And that gives us the maximal random-effects structure justified by this particular experimental design. By incorporating Age Group as a random slope within Item, we allow the effect of Age Group to vary across Items.
You should notice a nice symmetry here:
Congruency varies within Participant, so Participant gets a Congruency slope.
Age Group varies within Item, so Item gets an Age Group slope.
This gives us a useful question to keep coming back to: What varies within each random effect?
Yes! Funnily enough, one of my favourite R packages, performance (Lüdecke et al., 2021), has a handy little function called check_group_variation().
This function can help with identifying the random-effect structure for your model. That's not to say you don't have to think about your experimental design; you absolutely do! But it provides a useful sanity check when things start getting complicated. For our schema example, we could run:
performance::check_group_variation(schema_data,
select = c("Congruency", "AgeGrp"),
by = c("Participant", "Item"))
Here, we're essentially asking: How do Congruency and AgeGrp vary with respect to Participant and Item? We should then get a printout like this:
Group Variable Variation
Participant Congruency within
Participant AgeGrp between
Item Congruency between
Item AgeGrp within
... which, thankfully, confirms everything we just spent all that time figuring out! Congruency varies within Participant, but between Items. Age Group varies between Participants, but within Items. Woohoo!
Note: check_group_variation() isn't some magical function that can tell you what random effects to put into your model; it describes how variables vary with respect to your grouping factors. As such, you still need to understand your experimental design. But if you ever find yourself staring at your dataset thinking... "Wait... does this vary within Participant? Within Item? Between them? Both? What is going on!?" ...then check_group_variation() is a very useful starting point!
Putting my "teaching hat" on a little bit, I thought it would be good to also do some exercises where you get to apply what you've learned to some other scenarios where you need to figure out the maximal model for each. I have done it in the form of multiple-choice questions, with the correct answers listed below the Reference List of this post.
Exercise 1: A group of psychology students investigate how quickly adults identify facial emotions. Each participant sees many different faces, which are either presented upright or inverted. The same face stimulus can be presented upright or inverted, but not to the same participant. The outcome is reaction time in milliseconds for correctly identified emotions. The authors fit a generalised linear mixed-effects model, but are not sure what the maximal model structure should be.
A. RT ~ Orientation + (1 + Orientation | Participant) + (1 + Orientation | Face)
B. RT ~ Orientation + (1 + Orientation | Participant) + (1 | Face)
C. RT ~ Orientation + (1 | Participant) + (1 | Face)
D. RT ~ Orientation + (1 | Participant) + (1 + Orientation | Face)
Exercise 2: Some clinical psychologists are investigating therapists' judgements about patients potentially requiring further assessment. Each therapist reads multiple clinical vignettes, with vignettes manipulating symptom severity as either low or high. Therapists belong to one of two groups based on level of clinical experience: trainee or qualified. Notably, each vignette is presented in both severity versions, but across different therapists. The outcome is whether the therapist recommends further assessment. The authors fit a generalised linear mixed-effects model, but are not sure what the maximal model structure should be.
A. Judgement ~ Experience * Severity + (1 + Experience | Therapist) + (1 + Experience + Severity | Vignette)
B. Judgement ~ Experience * Severity + (1 + Experience * Severity | Therapist) + (1 + Experience * Severity | Vignette)
C. Judgement ~ Experience * Severity + (1 + Severity | Therapist) + (1 + Experience + Severity | Vignette)
D. Judgement ~ Experience * Severity + (1 + Severity | Therapist) + (1 + Experience * Severity | Vignette)
Exercise 3: A cognitive neuroscientist investigates whether reward influences episodic memory through activity in the striatum. During fMRI, participants encode a series of images, but before each image appears, participants are told that they can earn either a high, low, or no monetary reward for successfully remembering it later. Each image is assigned to only one reward condition across all participants. The authors fit a linear mixed-effects model but are not sure what the maximal model structure should be.
A. Activity ~ Reward + (1 + Reward | Participant) + (1 + Reward | Image)
B. Activity ~ Reward + (1 | Participant) + (1 + Reward | Image)
C. Activity ~ Reward + (1 + Reward | Participant) + (1 | Image)
D. Activity ~ Reward + (1 | Participant) + (1 | Image)
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
Lüdecke, D., Ben-Shachar, M., Patil, I., Waggoner, P., & Makowski, D. (2021). performance: An R Package for Assessment, Comparison and Testing of Statistical Models. Journal of Open Source Software, 6(60), 3139. https://doi.org/10.21105/joss.03139
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
Exercise 1: The correct answer was A. Orientation varies within Participant because each participant sees both upright and inverted faces. We therefore include a random slope for Orientation within Participant. Orientation also varies within Face because, across participants, the same face can be presented both upright and inverted. Therefore, we include a random slope for Orientation within Face.
Exercise 2: The correct answer was D. Severity varies within Therapist because each therapist judges both low- and high-severity vignettes. However, therapist experience does not vary within a given Therapist because each therapist is either a trainee or qualified. Of course, within a particular Vignette, both Experience and Severity vary because, across participants, each vignette is judged by trainee and qualified therapists and appears at both levels of Severity. Their interaction therefore also varies within Vignette. We therefore include random slopes for Therapist, Severity, and their interaction within Vignette.
Exercise 3: The correct answer was C. Reward varies within Participant because each participant experiences high-, low- and no-reward trials. However, somewhat strangely I must admit, reward does not vary within Image because each image is assigned to only one reward condition across everyone. Therefore, we include a random intercept, but no random slope for Reward, within Image.
Date: 23rd September 2026.
In my very first post on this blog, I discussed how categorical random slopes had to be correlated when using lme4 because using || simply didn't work. I then discussed a workaround using afex.
Well, I have NEWS...
As of version 2.0-1, lme4 now allows us to specify uncorrelated categorical random slopes directly. No workaround required! You simply need to add the following at the beginning of your script:
options(lme4.doublevert.default = "diag_special")
Then, using || will allow you to specify uncorrelated random slopes, including for categorical predictors. A tiny line of code!!
This is, of course, an exciting development for those of us who apparently have very specific ideas about what constitutes exciting news.
You're all welcome... I contributed absolutely nothing to this development.
Date: 2nd August 2025.
Recently, I ran into three quirks when fitting mixed-effects models in R. First up: how lme4 (Bates et al., 2015) handles uncorrelated (categorical) random slopes. Then I stumbled into a factor level fiasco courtesy of afex (Singmann et al., 2024). And just when I thought I had it all working... performance (Lüdecke et al., 2021) decides I no longer need all my usual diagnostics.
Just a quick reminder: If you're testing within-subject effects, it's best practice to include random slopes for those effects as they help control Type I error by preventing artificially low p-values that can occur when ignoring this variability (Barr et al., 2013; Barr 2013). That said, it doesn't mean it is easy to fit these models, especially when trying to avoid overparameterisation/fit issues! 😭
When you specify a model like: Y ~ X + (1 + X | Participant)... you’re (implicitly) telling lme4 to include both random intercepts and slopes, and to estimate their correlation. However, if the correlation between intercepts and slopes is extremely high (e.g., > +.80 or < –.80), your model may be overfitting or become unstable even in the absence of a warning/error message (see Bates et al., 2018; Matuschek et al., 2017).
A common solution is to remove the correlation but keep both the intercept and slope: Y ~ X + (1 + X || Participant). Here, the || syntax tells lme4 to model the slope and intercept as uncorrelated, helping reduce complexity without discarding important structure. However, this little trick only works if your slope is numeric! If your random slope is categorical? Tough luck. lme4 digs in its heels and refuses to drop the correlation, no matter how nicely you ask or how many tears you shed.
Turns out others have picked up on this quirk (see: Stack Exchange post from 2018). One suggested workaround is through the afex package using the afex::lmer_alt function. Through this function, you can model uncorrelated random slopes even when the predictors are categorical:
GLM = afex::lmer_alt(Y ~ 1 + X + (X || Participant) + (1 | Stimulus),
data = DT,
family = binomial(link = "logit"),
control = glmerControl(optimizer = "bobyqa"))
This allows full use of glmer functionality, but lets you get around lme4's issues with categorical predictors and ||.
However, there’s a subtle issue with afex::lmer_alt... it doesn’t automatically ignore unused factor levels. So, if you're anything like me (don't be like me!) and try to make code that works across experiments like so:
DT = read_csv(File_Location) %>%
mutate(X = factor(X,
levels = c("A", "B", "C", "D")))
... you're going to run into some issues if your data doesn’t contain data for all levels of a factor (e.g., no “C” trials in this dataset). In those cases, afex will still try to fit a slope for that level, even if there's not a single data point in sight. It's like trying to model a ghost! When this happens, you will get warnings like:
Warning messages: In vcov.merMod(object, use.hessian = use.hessian): variance-covariance matrix computed from finite-difference Hessian is not positive definite or contains NA values: falling back to var-cov estimated from RX
To fix this, make sure unused factor levels are dropped using the droplevels base R function (e.g., DT$X = droplevels(DT$X)) before you fit the model via afex. Basically, when preparing your data, be mindful of how you factorise your variables!
So, you have finally wrangled a model that fits with uncorrelated slopes and ghost levels removed, but you now want to check the model meets the necessary statistical assumptions, so you turn to your trusted diagnostic tool: performance::check_model. This usually gives a lovely set of plots to help assess key assumptions, including whether the random effects are normally distributed (a core assumption in mixed-effects models, whether you use a Gaussian distribution or not).
But plot twist! When you ran performance::check_model(GLM), the random effects distribution plots were missing. Nada. Gone. Like it never existed. To make sure you weren't imagining things, you ran the same model using lme4::glmer directly (with or without random slopes), and voilà, those helpful diagnostic plots returned. So it turns out that performance::check_model doesn't support the afex::lmer_alt workaround.
[I would insert a pretty picture here showing you the difference in outputs for the function, but it won't let me!]
Now I hear you asking: “But how do I check the random effects if performance::check_model won't do it?” Great question. Below is a workaround that gets the job done... but it isn't the most efficient way of doing things.
To assess whether your random effects are normally distributed, you’ll need the help of the lattice package (Sarkar, 2008). Specifically, we’ll use the lattice::qqmath function to create Q-Q plots manually. Here’s the basic idea:
Use lme4::ranef to extract the random effects from your model; this gives you a list of data frames, one for each grouping factor.
Pick the random effect you want to inspect (e.g., the slope for condition B might show up as re1.B. in the output of lme4::ranef.
These names can be a little cryptic, so use names(REs[["Participant"]]) to reveal them!
Use qqmath to plot the quantiles and check if they roughly follow a normal distribution.
Below is code that I used to do exactly that!
# Extract the random effects from the model (will be a list)
REs = lme4::ranef(GLM)
# Get the variable names for Participant effects
RE_Names = names(REs[["Participant"]])
# Generate Q-Q plot for random effects
lattice::qqmath(
# Random slope for level "B" of a categorical predictor (from ranef summary)
~ `re1.B.`,
# Grab it from the random effects list generated (specify which random effect!)
data = REs[["Participant"]],
# Make plot
panel = function(x, ...) {
# Plots sample vs. theoretical quantiles
panel.qqmath(x, ...)
# Adds reference line for normality
panel.qqmathline(x, col = "red", lty = 2)
}
)
To change the random effect you're looking at, just swap the name of the first argument used when calling lattice::qqmath. For example, you could put ~ `(Intercept)` to look only at the random intercept per participant. If you have another grouping factor (e.g., Stimulus), just change the list you are pulling from. For example, data = REs[["Stimulus"]], but make sure you have the correct variable names for the random effects in that list!
Mixed-effects modelling: come for the statistical power, stay for the weird late-night debugging! I hope that something in this post will be helpful to someone else besides me, even if only a little bit helpful.
Update: The bug in performance is now under investigation. You can find out more at the issue posted on GitHub: Normality of random effect not returned using ||.
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
Barr, D. J. (2013). Random effects structure for testing interactions in linear mixed-effects models. Frontiers in Psychology, 4. https://doi.org/10.3389/fpsyg.2013.00328
Bates, D., Mächler, M., Bolker, B., & Walker, S. (2015). Fitting Linear Mixed-Effects Models Using lme4. Journal of Statistical Software, 67(1), 1–48. https://doi.org/10.18637/jss.v067.i01
Bates, D., Kliegl, R., Vasishth, S., & Baayen, H. (2018, May 26). Parsimonious Mixed Models. ArXiv. https://doi.org/10.48550/arXiv.1506.04967
Lüdecke, D., Ben-Shachar, M., Patil, I., Waggoner, P., & Makowski, D. (2021). performance: An R Package for Assessment, Comparison and Testing of Statistical Models. Journal of Open Source Software, 6(60), 3139. https://doi.org/10.21105/joss.03139
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
Sarkar, D. (2008). Lattice: Multivariate Data Visualization with R. Springer: New York. https://doi.org/10.1007/978-0-387-75969-2
Singmann, H., Bolker, B., Westfall, J., Aust, F., & Ben-Shachar, M. (2024). afex: Analysis of Factorial Experiments. R package version 1.3-1. https://CRAN.R-project.org/package=afex