# Isn't it normal that residual plots for mixed effect models will show a trend?

Whenever we include a random effect in a mixed model, since the estimates are shrunk, isn't it normal at the fundamental level that the residual plot will show a positive correlation between fitted values and residuals?

Here's a super simple model with no fixed effect and, as random effect, 500 subjects observed 4 times:

set.seed(1)
N = 500
R = 4

subjects = rep(factor(1:N),each = R)
mu = rep(runif(N),each = R)
y = mu + 0.2*rnorm(N*R)
df = data.frame(subjects,mu,y)

mdl = lme4::lmer(y ~ 1 + (1|subjects),df)
performance::check_model(mdl)


It feels intuitive to me that estimates being shrunk, subjects with higher (lower) $$\mu$$ will show a positive (negative) residual.

Is it a problem? How can we diagnose (the linearity of) our models then?

Note: the trend exists also with a normal distribution of $$\mu$$s.

The random effect is assumed to follow a normal distribution, so if you simulate it as uniformly distributed, it's hard to pinpoint why you observe something unexpected. Nevertheless, as you indicated, the positive slope persists if you use a normal distribution, which I will use here to suggest a solution:

library("lme4")
set.seed(1)
N <- 500
R <- 4
subjects <- rep(factor(1:N),each = R)
mu <- rep(rnorm(N), each = R)
y <- mu + 0.2 * rnorm(N * R)
df <- data.frame(subjects, mu, y)
mdl <- lme4::lmer(y ~ 1 + (1|subjects), df)
plot(resid(mdl) ~ fitted(mdl), col = "steelblue", pch = 16)
abline(coef(lm(resid(mdl) ~ fitted(mdl))), col = 2, lwd = 2)


One solution, provided by the DHARMa package, is to use simulated residuals conditional only on the fixed effects:

library("DHARMa")
simres <- simulateResiduals(mdl, re.form = ~ 0)


Compare this to simulated residuals conditional on the fixed and random effects:

simres2 <- simulateResiduals(mdl, re.form = ~ 1 | subjects)


• I am unsure how this alleviates the issue of correlation between random effects and residuals, which are assumed to be independent in a random effects model. Commented Aug 2 at 18:47
• @larmor The plot in the question shows the residuals against the fitted values, which are not only comprised of the fixed effects, but also the random effects. In other words, the correlation you observe is between $\hat{\mathbf{\epsilon}}$ and $\mathbf{X\hat\beta} + \mathbf{Z \hat \upsilon}$. If you draw any conclusions about the fixed effects, the random effects and $\epsilon$ are assumed zero. The plot shown here only compares $\hat\epsilon$ to $\mathbf{X\hat\beta}$ (using a simulation-based approach) and is therefore a better way to assess issues with the fixed effects, like non-linearity. Commented Aug 2 at 19:01
• I agree with all your statements, but I'm also worried about the assumption of independence between random effects and residuals. The plot in my question is displaying a correlation between them (since there is no fixed effect to correlate with). Commented Aug 3 at 10:59
• After a few days, my understanding is that while it is assumed in a random effect model that μs and ϵs are independent variables, it is not the case that conditional modes (i.e. individual estimates of random effects) and residuals will be uncorrelated. Commented yesterday