Multilevel Model Residuals Scatterplot Assumptions

Question

I am conducting multilevel modelling (MLM) in SPSS (mixed modeling) to analyze cross-sectional repeated measures data. One of my dependent variables is a survey question scaled 1 to 10, which produced the attached scatterplot of the residuals, which is angled downward but still balanced around the zero line. Could this "angling" cause any concern for interpreting that the MLM assumptions are being met?

Answer 1

You only have a fixed number of outcomes, so the pattern you are observing reflects that.

In fact, you can quite closely reproduce this pattern by plotting residuals of a binomial process:

...or the Pearson residuals of a multilevel model:

What to make of this?

What you're really observing is that:

• If your model assumes a normally distributed error term, this assumption is not reasonable, and you should use a GLMM instead.
• For models that do allow discrete probability distributions, a simple residuals vs fitted plot is not very informative.
• For multilevel models, you have to think about what assumptions are actually being made, and how you might go about diagnosing those.
• If you are willing to use R, all of these considerations can be taken into account using the amazing package DHARMa.

---

Simulated data to reproduce the figures:

# Binomial example
set.seed(2024)
n    <- 1000
size <- 10
x    <- rnorm(n, 10, 3)
ID   <- factor(rep(1:100, each = 10))
Z    <- model.matrix(~ 0 + ID)
eta  <- -1.8 + 0.2 * x + Z %*% rnorm(ncol(Z))
y    <- rbinom(n, size, binomial()$linkinv(eta)) + 1
DF   <- data.frame(x, y)
rm(n, size, x, eta, y)
LM <- lm(y ~ x, DF)

plot(LM, which = 1, pch = 21, bg = "steelblue1", xlim = c(2, 10))

# Binomial multilevel example
set.seed(2024)
n    <- 1000
size <- 10
x    <- rnorm(n, 6, 3)
ID   <- factor(rep(1:100, each = 10))
Z    <- model.matrix(~ 0 + ID)
eta  <- -0.75 + 0.15 * x + Z %*% rnorm(ncol(Z), 0, 0.5)
y    <- rbinom(n, size, binomial()$linkinv(eta)) + 1
DF   <- data.frame(x, y, ID)
rm(n, size, x, eta, y)

LM <- lm(y ~ x, DF)

require("lme4")
LMM <- lmer(y ~ x + (1 | ID), DF)

plot(residuals(LMM) ~ fitted(LMM), pch = 21, bg = "steelblue1", xlim = c(2, 10),
     cex = 0.7)
mtext("Pearson Residuals vs Fitted", 3, 0.5)