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

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?