Recreating three-way interaction plot from sjplot in ggplot

Question

I used sjPlot to plot a significant three-way interaction from a glm model using a Gamma distribution. The DV is continuous, two of the IVs are continuous, and one is categorical. I'd love to have the plot be faceted by the categorical variable, with lines on each that show mean +/- 1 sd (shown below). I got a very nice facet-wrapped plot and am able to modify the axis limits, titles, etc. However, I don't know how to change certain aesthetics, such as removing the axis lines/grids, making the background color white, changing the color of the bands at the top of the facets, etc. I can do this in ggplot, but I'd really love to use the specific predicted values from the model because it contains covariates in addition to the IV, DV, and moderators.

Here's an example using the iris dataset:

iris.glm<-glm(Sepal.Length ~ Sepal.Width*Petal.Width*Species + Petal.Length, 
              family = Gamma(link="inverse"), data=iris)
summary(iris.glm)
plot_model(iris.glm, type="int", mdrt.values="meansd", axis.lim=c(4,10), 
           title="", axis.title=c("Sepal Width", "Sepal Length"), 
           legend.title="Petal Width", colors=c("#66c2a5","#fc8d62", "#8da0cb"))

However, I'd like to do more refining of the aesthetics. When I try to recreate this plot in ggplot, I'm doing a terrible job! In particular I don't know how to essentially convert the continuous moderator to a three-level category (-1 SD, Mean, 1 SD). Do I need to create a separate data frame with predicted values from the model and create a plot from that? If so, could someone help with that? Apologies in advance, still a novice with ggplot.

iris %>% 
  ggplot(aes(x=Sepal.Width, y=Sepal.Length, color=Petal.Width, group=Petal.Width)) + 
  facet_grid(~Species) +    
  stat_smooth(method = "glm", se=T) + 
  xlab("Sepal Width") + 
  ylab("Sepal Length") + 
  coord_cartesian(xlim =c(2, 4.5), ylim = c(4, 10))

Answer 1

The sjPlot author is very helpful and he lurks here from time to time, so you might want to be a bit patient and wait for a direct answer.

In the meantime, you can also try the plot_cap() function from marginaleffects (disclaimer: I am the author). The nice thing is that this function returns “normal” ggplot2 objects, so you can customize using all the normal aesthetics and functions like theme_classic() or ylim():

library(ggplot2)
library(marginaleffects)

iris.glm <- glm(
    Sepal.Length ~ Sepal.Width * Petal.Width * Species + Petal.Length,
    family = Gamma(link = "inverse"), data = iris)

plot_cap(
    iris.glm,
    condition = list(
        "Sepal.Width",
        "Petal.Width" = "threenum",
        "Species")) +
    ylim(c(4, 10)) +
    theme_classic()

Also, you can return the underlying data too by setting draw=FALSE. Then you can plot that yourself if you need more customization:

p <- plot_cap(
    iris.glm,
    draw = FALSE,
    condition = list(
        "Sepal.Width",
        "Petal.Width" = "threenum",
        "Species"))

head(p)
#>   rowid     type estimate std.error statistic      p.value  conf.low conf.high
#> 1     1 response 5.534242 0.4372848 12.655920 1.037525e-36  6.548356  4.792110
#> 2     2 response 5.122283 0.3618975 14.153961 1.765283e-45  5.945597  4.499251
#> 3     3 response 4.466816 0.4016798 11.120340 9.989003e-29  5.422541  3.797504
#> 4     4 response 6.827409 2.8194809  2.421513 1.545605e-02 35.819951  3.773308
#> 5     5 response 5.282150 0.1117132 47.283120 0.000000e+00  5.510572  5.071911
#> 6     6 response 4.739130 0.2387601 19.848918 1.125985e-87  5.258362  4.313224
#>   Sepal.Length Petal.Length condition1 condition2 condition3
#> 1     5.843333        3.758          2        -SD     setosa
#> 2     5.843333        3.758          2        -SD versicolor
#> 3     5.843333        3.758          2        -SD  virginica
#> 4     5.843333        3.758          2       Mean     setosa
#> 5     5.843333        3.758          2       Mean versicolor
#> 6     5.843333        3.758          2       Mean  virginica