Interpretation of relationship between residuals and predicted values in a mixed model

Question

I'm using a linear mixed effects model (varying slope + intercept, with glmer from "lme4") to look at the association of the interaction between two categorical predictors (time + X1 - where X1 is a grouping variable) with a continuous outcome (Y1, has both positive and negative values)

The model looks like this:

glmer(Y1~ Time*X1+ (1 + Time+ X1 | Subject))

[Subject is to account for random intercepts in the individual subjects)

The problem is, my residuals vs predictors plot looks like this:

The way I interpret (based on this) this is that there's little heteroscedasticity here, but the model is biased. (in the sense that there's a relationship between the residuals + predictors).

This source tells me I might be missing a predictor, but I can't think of any that I'm missing.

So my questions are:

1. How big of a problem is this in terms of interpreting the model output?

2. Could it have something to do with the glm family that I'm using ("Gaussian")? Is this appropriate for such a continuous outcome variable?

3. What should/can I do to improve the model?

Thanks for the tip, Florian. This is what the plot looks like when I use "re.form = NULL" in predict():

Looks exactly the same. Any ideas why?

Answer 1

Per default, fitted() uses the full model predictions, including the REs. Plotting residuals against these conditional predictions often creates spurious patterns that have nothing to do with a model misspecification, see a bit more detail here https://github.com/florianhartig/DHARMa/issues/43.

Solution is to plot residuals against the fixed effects predictions only, adding re.form = NULL to the predict function - I would only worry if you still have a pattern.

Additionally, you should check REs against predictors though.