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.
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