I’m trying to visualise my adjusted lmer model showing an interaction of education status using the predict function and ggplot, but I get a zig-zag output when I add in my covariates. I want to include 95% CI bands as well.

I managed to plot my unadjusted model using the following code:

Unadjusted Model

lmer3<- lmer(outcome~time_years*educ+ before_after*educ+ time_since_job*educ + (1 | pidp), data = dat)
summary(lmer3) 

predictions <- predict(lmer3, newdata = dat, re.form = NA, type = "response", se.fit = TRUE)
dat$predicted_values <- predictions$fit
# Getting 95% CI bands
dat$lower_bound <- dat$predicted_values - 1.96 * dat$se_predicted
dat$upper_bound <- dat$predicted_values + 1.96 * dat$se_predicted

GGplot Code

ggplot(dat, aes(x = time_years, y = predicted_values, color = as.factor(educ))) +
  theme_classic()+
    geom_line() + ##line showing predicted values from model 
    labs(x = "Time (years)", y = "Daily MET-min", color = "Education ") +
    geom_ribbon(aes(ymin = lower_bound , ymax = upper_bound), alpha = 0.3, linetype= 0,fill = "grey")+
    geom_vline(xintercept = 0, linetype = "dashed", color = "red") +
    scale_color_manual(values = c("blue", "red"), label=c("No Degree", "Degree or higher")) + 
    geom_vline(xintercept = 0,linetype = "dashed", color = "red") +
      labs(x = "Time relative to job start (years)", y = "Daily Exercise", color = "Education Status")

unadjusted graph

But when I add in my covariates, I get a messy and nonsensical graph:

Adjusted Model

lmer4<- lmer(outcome~time_years*educ+ before_after*educ+ time_since_job*educ + sex_dv + ethnicity + parentEducation + (1 | pidp), data = dat)
summary(lmer4)

predictions<-predict(lmer4, newdata = dat, re.form = NA, type = "response", se.fit = TRUE)

Adjusted Graph

I’m not sure why this is happening, or how to fix it. Thanks in advance for any suggestions!