library(marginaleffects)
library(emmeans)
set.seed(20220722)
S <- diag(3)
S[1,2] <- S[2,1] <- 0.6
S[1,3] <- S[3,1] <- 0.6
data <- MASS::mvrnorm(1000, rep(0, 3), S) |>
as.data.frame() |>
# make outcome, group, baseline binary
transform(outcome = V1 < 0,
rx = V2 < 0,
baseline = V3 < 0) |>
subset(select = c(outcome, rx, baseline))
cor(data)
#> outcome rx baseline
#> outcome 1.0000000 0.40407111 0.40697644
#> rx 0.4040711 1.00000000 0.04377112
#> baseline 0.4069764 0.04377112 1.00000000
m1 <- glm(outcome ~ rx, data = data,
family = binomial())
exp(coef(m1)["rxTRUE"]) # OR
#> rxTRUE
#> 5.559922
marginaleffects(m1) |> # marginal difference
summary()
#> Average marginal effects
#> Term Contrast Effect Std. Error z value Pr(>|z|) 2.5 % 97.5 %
#> 1 rx TRUE - FALSE 0.4038 0.02891 13.97 < 2.22e-16 0.3471 0.4605
#>
#> Model type: glm
#> Prediction type: response
m2 <- glm(outcome ~ rx + baseline, data = data,
family = binomial())
exp(coef(m2)["rxTRUE"]) # OR is much larger!
#> rxTRUE
#> 7.93829
But marginal differences is about the same
marginaleffects(m2, variables = "rx") |>
summary()
#> Average marginal effects
#> Term Contrast Effect Std. Error z value Pr(>|z|) 2.5 % 97.5 %
#> 1 rx TRUE - FALSE 0.3867 0.02668 14.49 < 2.22e-16 0.3344 0.439
#>
#> Model type: glm
#> Prediction type: response
Created on 2022-07-20 by the reprex package (v2.0.1)
Bayesian equivalent with manual computation of marginal effect: