How to plot the marginal effect of a standard deviation increase and a standard deviation decrease both on the same plot?

Question

I would like to plot the marginal effect of a standard deviation increase and a standard deviation decrease of a model on the same plot. I am using the marginaleffects package, but I am open to other solutions.

This is an example, suppose I have the following model:

library(marginaleffects)
library(ggplot2)

mod <- lm(mpg ~ hp * disp + factor(am), data = mtcars)

I would like to get the standard deviation increase and standard deviation decrease of hp. The standard deviation is:

sd(mtcars$hp)
> sd(mtcars$hp)
[1] 68.56287

Then, I can plot individually each standard deviation using the codes:

plot1 <- plot_comparisons(mod, variables = list(hp = 68.56287), condition = "disp")
plot1

plot2 <- plot_comparisons(mod, variables = list(hp = -68.56287), condition = "disp")
plot2

Instead of having these values in two plots, I would like to get them in one plot. How can I do that?

EDIT: Following Vincent suggestion, I think I did it:

EDIT2: There was a mistake in previous code I posted, in which I had already standardized the data then I was calculating the plot with the standard deviation again (Thank you for pointing that out, Daniel). I corrected the codes.

plot1 <- plot_comparisons(mod, variables = list(hp = 68.56287), condition = "disp", draw = F)
plot1
plot2 <- plot_comparisons(mod, variables = list(hp = -68.56287), condition = "disp", draw = F)
plot2

plot1$example <- 1
plot2$example <- 2

plots <- rbind(plot1, plot2)

ggplot(plots, aes(disp, estimate, fill=factor(example), color = factor(example)))+
  geom_line()+
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high, color=NULL), alpha=0.2)

Now I have to figure out how to fix the labels, then I am good.

Answer 1

One solution is to use draw=FALSE, combine the datasets, and use ggplot2:

library(marginaleffects)
library(ggplot2)

mtcars2 <- within(mtcars, hp <- scale(hp)[,1])
mod <- lm(mpg ~ hp * disp + factor(am), data = mtcars2)
plot1 <- plot_comparisons(mod, variables = list(hp = 68.56287), condition = "disp", draw = FALSE)
plot2 <- plot_comparisons(mod, variables = list(hp = -68.56287), condition = "disp", draw = FALSE)
plot1$hp <- 69
plot2$hp <- -69
plots <- rbind(plot1, plot2)

ggplot(plots, aes(x = disp)) +
    geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = factor(hp)), alpha = 0.2) +
    geom_line(aes(y = estimate, color = factor(hp))) +
    labs(
        x = "Displacement",
        y = "Effect of a change in horsepower\non predicted miles per gallon",
        color = "Horsepower",
        fill = "Horsepower") +
    scale_color_manual(values = palette("Okabe-Ito")[2:3]) +
    scale_fill_manual(values = palette("Okabe-Ito")[2:3]) +
    theme_bw()

Answer 2

You could use the ggeffects-package, which let's you easily create plots of predictions or marginal effects at specific values:

library(ggeffects)

mtcars2 <- within(mtcars, hp <- scale(hp)[,1])
mod <- lm(mpg ~ hp * disp + factor(am), data = mtcars2)

# focal terms specified as strings
ggpredict(mod, c("disp", "hp[-68.5, 68.5]")) |> plot()

# or as list
ggpredict(mod, list(disp = seq(100, 500, 10), hp = c(-68.5, 68.5))) |> plot()