How to plot a two-way interaction with two categorical predictors from a logistic regression model in R?

Question

I am trying to plot in R a two-way interaction with two dichotomous predictors from a logistic regression model (i.e., the DV is dichotomous as well) such that the y-axis will present the probabilities and the predicted values with have SE bars.

I've tried using afex_plot:

glm <- glm(Y ~ X*Z, data = Data,
            family = "binomial")
library("afex")  
afex_plot(glm, X", "Y")

and received the following error:

Error in tbl[, vars, drop = FALSE] : incorrect number of dimensions

Then I tried to use ggpredict and plot functions and succeed to generate the desired plot. However, ggpredict with a binomial model allows you to present on the plot only the CIs and not the SE bars:

mydf <- ggpredict(glm, c("X","Z"))
plot(mydf)

Then I tried to use the SEs that ggpredict saves to calculate the error bars manually but then realized that the standard errors are always on the link-scale, and not back-transformed for non-Gaussian models!

How do I plot the desired graph?

Thanks in advance, Amit

Answer 1

Not sure why ggpredict isn't giving you what you want. Here's what I get when using it and note, everything is automatically on the response scale.

set.seed(32043)
dat <- tibble(
X = sample(c("L1", "L2"), 250, replace=TRUE),
Z = sample(c("L1", "L2"), 250, replace=TRUE)
)
b <- c(-1, 1, -2, 3)
mat <- model.matrix(~X*Z, data=dat)
p <- plogis(mat %*% b)
dat$y <- rbinom(250, 1, p)
mod <- glm(y ~ X*Z, data = dat, family="binomial")

library(ggeffects)
g <- ggpredict(mod, c("X", "Z"))
g
# # Z = L2
# 
# x  | Predicted |   SE |       95% CI
# ------------------------------------
# L2 |      0.84 | 0.34 | [0.73, 0.91]
# L1 |      0.07 | 0.46 | [0.03, 0.16]
# 
# # Z = L1
# 
# x  | Predicted |   SE |       95% CI
# ------------------------------------
# L2 |      0.47 | 0.25 | [0.35, 0.59]
# L1 |      0.34 | 0.29 | [0.23, 0.48]

Then, you could just use the built-in plot method for these objects.

plot(g)

If you wanted the lines like you get from afex_plot(), you could add them with the appropriate geometry.

plot(g) + 
  geom_line(aes(linetype=group), 
            position=position_dodge(width=.25))

---

Edit: SE bars instead of CI

I misunderstood the original post - where specifically SE bars are desired. The way that I would normally do this is with a parametric bootstrap. First do a parametric bootstrap of the coefficients

b <- MASS::mvrnorm(1000, coef(mod), vcov(mod))

Then, make the data that defines the four predictions - all of the combinations of the two values of X and Z:

newdat <- expand.grid(
  X = factor(c(1,2), labels=c("L1", "L2")),
  Z = factor(c(1,2), labels=c("L1", "L2")), 
  y=0
)

Then, make the model matrix and generate the predicted probabilities for all of the bootstrapped coefficients:

Xnew <- model.matrix(mod, data=newdat)

probs <- plogis(Xnew %*% t(b))

Next, put the predictions and the standard errors (calculated across all of the bootstrapped probabilities for each condition) into the prediction data frame.

newdat$fit <- predict(mod, newdata=newdat, type="response")
newdat$se <- apply(probs, 1, sd)

Then, you could make the plot adding and subtracting a standard error from each fitted probability:

ggplot(newdat, aes(x=X, y=fit, colour=Z)) + 
  geom_errorbar(aes(ymin=fit-se, ymax=fit+se), 
                width=.05, 
                position=position_dodge(width=.25)) + 
  geom_point(position=position_dodge(width=.25)) + 
  geom_line(aes(x=as.numeric(X), linetype=Z), 
            position=position_dodge(width=.25)) + 
  theme_bw()