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.
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()