Understanding residuals vs. fitted plot for a linear mixed model

Question

I am modeling body mass (y var) according to indices of dysregulation for different physiologic systems (x vars). I did a likelihood ratio test, which supported using a linear mixed model, with a random effect of individual, over a simple linear model.

library(lme4)
require(RLRsim)
mod1 <- lm(mass_kg_zscore ~ metabolic_index + immune_index + neuro_index, data = pb)
mod2 <- lmer(mass_kg_zscore ~ metabolic_index + immune_index + neuro_index + (1|id), data = pb, REML = FALSE)
exactLRT(mod2,mod1) # individual, null - tests individual 
# P =0.0026, so there is support for keeping the RE of individual

However, the residuals vs. fitted plot for the linear mixed model looks pathologic:

plot(fitted(mod2),residuals(mod2))
abline(h = 0)

The QQ plot for mod2 looks relatively okay:

Now, if I remove the random effect of individual, and just do a simple linear model, my residuals vs. fitted plot looks much better:

My question is, why does the inclusion of the random effect cause such a wonky residuals vs. fitted plot? Does this mean I should exclude the RE, even though the likelihood ratio test supported keeping the RE in the model?

Thank you!

Answer 1

My question is, why does the inclusion of the random effect cause such a wonky residuals vs. fitted plot?

It's quite hard to say, without having the data itself.

First note that the bulk of range of the residuals for the mixed model is around -1 to +1 (ignoring a few outliers). However, for the linear model, the bulk of the range of the residuals is -2 to +2. This is telling you that the random effects are reducing the residual variance, which is to be expected.

Your main concern is that there appears to be an upward trend in the residuals of the mixed model. This might just be a random artefact. Do you have more data ? If so, it would be very interesting to see it was reproduced with further data. Note that in the linear model (lm) you have not adjusted for the effects of id. There is clearly some variation attributable to it, so it would be a good idea to fit id in the linear model, as a fixed effect and then compare the plots again.

Does this mean I should exclude the RE, even though the likelihood ratio test supported keeping the RE in the model?

No, I don't think so. I would be inclined to retain the random intercepts, and look into fitting a robust linear mixed model using package robustlmm