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!

If you want two slopes, then you need to remember that you need two x axes; the slopes represent a change in the outcome variable per unit change in the predictor variable on the x axis, and two slopes means two predictor variables, which means two x axes. This would make for a pretty confusing plot. Wouldn't it be better to use multiple panels here? — Allan Cameron, Sep 02 '22 at 13:24
Oh, right! Unlike in the mtcar example, my predictors are all on the same 6-point scale. So the x-axes look identical, just the slopes are at a differnt angle — format_0, Sep 02 '22 at 13:28

score 0 · Answer 1 · answered Sep 02 '22 at 14:03

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')

Awesome! Thank you so much for this elegant solution. I wonder if this could be made even easier by extracting the predicted values and the SEs out of the plot_model object. I mean, not necessarily easier, as your solution does not even require the plot_model package, so ... — format_0, Sep 02 '22 at 14:31

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

1 Answers1