Plotting two(!) regression lines from a linear regression model ( > 3 predictors) into the same plot with confidence intervals (sjPlot or ggplot)

Question

here's my problem:

I have a linear regression model with 5 predictors, all of which use the same 6-point scale. I would like to plot two of the slopes (predictions) into one plot, including confidence intervals. Thus far, I've been trying to use sjPlot (which is based on ggpredict I believe) to achieve this. I managed to create two separate plots, but cannot figure out how to combine them into one plot -- ideally with a legend naming the two predictors.

Here's an example code:

library(sjPlot)
dat = mtcars
fit = lm(mpg ~ disp + hp + drat + wt + qsec, data = dat)

plot_model(fit, type = "eff", terms = c("disp"))
plot_model(fit, type = "eff", terms = c("hp"))

I tried terms = c("disp", "hp") but this makes plot_model try to plot the interaction -- the effect of disp at different levels of hp.

Thanks!

Answer 1

The plots in plot_model are made by getting a model's predictions and standard errors across the range of the variable in question, while holding all other covariates at their mean value. In your example, that means we can effectively recreate your two plots in vanilla ggplot by doing:

dat = mtcars
fit = lm(mpg ~ disp + hp + drat + wt + qsec, data = dat)

hp_preds <- predict(fit, newdata = data.frame(hp = 52:335, 
                                              disp = mean(mtcars$disp), 
                                              drat = mean(mtcars$drat), 
                                              wt = mean(mtcars$wt), 
                                              qsec = mean(mtcars$qsec)), 
                    se = TRUE)

ggplot(cbind(hp = 52:335, as.data.frame(hp_preds[1:2])), aes(hp, fit)) +
  geom_ribbon(aes(ymin = fit - 1.96*se.fit, ymax = fit + 1.96*se.fit),
              alpha = 0.15) +
  geom_line(size = 1) +
  ggtitle('Predicted values of mpg')

And

disp_preds<- predict(fit, newdata = data.frame(hp = mean(mtcars$hp), 
                                               disp = 71:472, 
                                               drat = mean(mtcars$drat), 
                                               wt = mean(mtcars$wt), 
                                               qsec = mean(mtcars$qsec)), 
                    se = TRUE)

ggplot(cbind(disp = 71:472, as.data.frame(disp_preds[1:2])),
       aes(disp, fit)) +
  geom_ribbon(aes(ymin = fit - 1.96*se.fit, ymax = fit + 1.96*se.fit),
              alpha = 0.15) +
  geom_line(size = 1) +
  scale_y_continuous(breaks = c(15, 20, 25), limits = c(15, 27.8)) +
  ggtitle('Predicted values of mpg')

This shows that we can simply overlay the two plots, as long as the x axis is specific about the fact that the units are different for the variables:

ggplot(cbind(disp = 71:472, as.data.frame(disp_preds[1:2])), aes(disp, fit)) +
  geom_ribbon(aes(ymin = fit - 1.96*se.fit, ymax = fit + 1.96*se.fit,
                  fill = 'disp'), alpha = 0.15) +
  geom_line(size = 1, aes(color = 'disp')) +
  geom_ribbon(data = cbind(hp = 52:335, as.data.frame(hp_preds[1:2])),
              aes(x = hp, ymin = fit - 1.96*se.fit, ymax = fit + 1.96*se.fit,
                  fill = 'hp'), alpha = 0.15) +
  geom_line(data = cbind(hp = 52:335, as.data.frame(hp_preds[1:2])),
            size = 1, aes(x = hp, color = 'hp')) +
  labs(fill = 'Predictor', color = 'Predictor', 
       x = 'disp or hp value', y = 'mpg') +
  ggtitle('Predicted values of mpg')