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

Question

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.

Answer 1

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)

Answer 2

In regards to your question, and also to clarify the answer provided by Frans, I want to point out that a res ~ fitted plot has 2 axis, a) the x axis where you have the model predictions, and b) the y axis, where you have the residuals. For both axis, you can condition on the REs or not.

1. When you condition on the REs when calculating the predictions on x, you can/will indeed get get the pattern that you highlight, in particular when the residuals are unconditional (but I believe there is a light pattern as well if you calculate conditional residuals). I show this here https://github.com/florianhartig/DHARMa/issues/43. For this reason, predicted values (on x) in DHARMa are always marginal predictions (= unconditional, not including the REs), and I believe it is good advice to NOT use res ~ FullPredictions, and rather use res ~ MarginalPredictions for mixed models.

2. The re.form argument that Frans mentions determines in DHARMa how the residuals are calculated (i.e. the values on y, not on x). As shown in the GH issue above, this affects the correlation against conditions predictions, but I think this point is secondary to your question.