I am running a linear mixed lmer of the following form:
Y ~ (1|ID) + Age + Sex + Education + Time*Factor + Time*Diagnosis

My control term is control <- lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000))

I have about 8800 observations on about 2100 individuals. Variable amounts of missing data.

I use predict_response from the ggeffects package to obtain marginal mean predictions for Time, Factor, Diagnosis. To visualize the model fit, I plotted the actual mean ± standard error from the data against the predicted R with some homebrew ggplot code to graphically assess the model results.

I have about 15 different Y variables. All my models converged easily. But some dependent variables have biased predictions. See figure below:

Any help addressing the issue or model terms is appreciated. I cannot use a random slope because of sampling limitations.