A randomized experiment is assumed. Y is the binary indicator of the poor outcome, E is the exposure indicator, and M is the effect modifier. As it is a randomized experiment M is blanced across levels of the assigned exposure (no E-M association). M also does not affect the outcome among the unexposed. The effect of E is beneficial in direction among patients in the M = 0 subgroup (non-significant trend), it is harmful in direction among patients in the M = 1 subgroup (non-significant trend). The interaction term is significant suggesting a statistically significant effect modification. What should we conclude in this case?
dat <- expand.grid(Y = 0:1,
E = 0:1,
M = 0:1)
dat$count <- c(5,5,6,4, 5,5,4,6)*10
xtabs(count ~ E + Y + M, dat)[2:1,2:1,]
## , , M = 0
##
## Y
## E 1 0
## 1 40 60
## 0 50 50
##
## , , M = 1
##
## Y
## E 1 0
## 1 60 40
## 0 50 50
## No M->E association
xtabs(count ~ M + E, dat)[2:1,2:1]
## E
## M 1 0
## 1 100 100
## 0 100 100
## Among unexposed, no M->Y association
xtabs(count ~ E + Y + M, subset(dat, E == 0))[2:1,2:1,]
## Error in xtabs(count ~ E + Y + M, subset(dat, E == 0))[2:1, 2:1, ]: subscript out of bounds
## Reflevel of M = 0
model1 <- glm(formula = Y ~ E*factor(M, 0:1),
weights = count,
family = binomial(link = "logit"),
data = dat)
tableone::ShowRegTable(model1)
## exp(beta) [confint] p
## (Intercept) 1.00 [0.67, 1.48] 1.000
## E 0.67 [0.38, 1.17] 0.156
## factor(M, 0:1)1 1.00 [0.57, 1.74] 1.000
## E:factor(M, 0:1)1 2.25 [1.02, 4.99] 0.045
## Reflevel of M = 1 (to obtain the E effect and 95% CI among M = 1)
model2 <- glm(formula = Y ~ E*factor(M, 1:0),
weights = count,
family = binomial(link = "logit"),
data = dat)
tableone::ShowRegTable(model2)
## exp(beta) [confint] p
## (Intercept) 1.00 [0.67, 1.48] 1.000
## E 1.50 [0.86, 2.64] 0.156
## factor(M, 1:0)0 1.00 [0.57, 1.74] 1.000
## E:factor(M, 1:0)0 0.44 [0.20, 0.98] 0.045
library(ggplot2)
library(magrittr)
datPlot <- rbind(exp(cbind(coef(model1), suppressMessages(confint(model1))))[2,],
exp(cbind(coef(model2), suppressMessages(confint(model2))))[2,]) %>%
as.data.frame
names(datPlot) <- c("OR","lower","upper")
datPlot$M <- 0:1
ggplot(data = datPlot, mapping = aes(x = factor(M), y = OR, ymin = lower, ymax = upper)) +
layer(geom = "errorbar", width = 0.1) +
layer(geom = "point") +
layer(geom = "hline", yintercept = log10(1), size = 0.1) +
labs(title = "p for interaction is significant") +
scale_y_log10(limit = c(0.1,10)) +
coord_flip() +
theme_bw() + theme(legend.key = element_blank())