Constrasts in R after logistic regression with interaction

Question

I would like to calculate contrasts of predictive probabilities. In a nutshell, I would like to compute the difference in differences of probabilities. The code below illustrates exactly my problem.

Among men, I calculate the predictive probability of the event by smoking:

p1 = prob. of event|sex==1 and smoking == 1;

p2 = prob. of event|sex==1 and smoking == 0

Among women, I calculate the predictive probability of the event by smoking.

p3 = prob. of event|sex==0 and smoking == 1;

p4 = prob. of event|sex==0 and smoking == 0

Now, what I want to calculate is the (p4-p3)-(p2-p1) along with 95% CI. In the example below, the difference in differences would be around -0.06 [e.g., (0.44-0.5)-(0.52-0.52)]. How do we estimate the SE/95% CI for that difference in R?

library(ggeffects)
N = 100
set.seed(123456)
hypertension <- round(runif(N))
age <- rnorm(N,45,5)
sex <- round(runif(N))
smoking <- round(runif(N))
model <- glm(hypertension ~ age + sex + smoking +sex*smoking, family = "binomial")
summary(model)
ggpredict(model, c("smoking", "sex")) 

Thank you so much for any tips/suggestions.

All the best,

Jacob

Answer 1

The bootstrap will probably be your easiest route here. Something like

df = data.frame(hypertension, age, sex, smoking)

stat = function(dat, j) {
  dat = dat[j,]
  model <- glm(hypertension ~ age + sex + smoking +sex*smoking,
    family = "binomial", data = dat)
  pred <- ggpredict(model, c("smoking", "sex"))
  pred <- pred$predicted    
  (pred[4] - pred[2]) - (pred[3] - pred[1]) # double-check the indexing here!
}

library(boot)
br <- boot(df, stat, R = 1000)
boot.ci(br)

would work. By the way, I think the ggpredict call is something of a bottleneck in the above code, so I would write your own function to compute the predictions (preferably using predict(model, ...) to handle the interactions).