binomial baseline control in RCT.md

library(marginaleffects)
library(emmeans)

Make some data

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

Fit unadjusted model

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

Fit adjusted model

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:

library(brms)

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))


m <- brm(outcome ~ rx + baseline, data = data,
         family = bernoulli(),
         backend = "cmdstanr")


diff <- 
  # Predicted value if all were in the treatment group
  data |> transform(rx = TRUE) |> 
  posterior_epred(m, newdata = _) - 
  
  # Predicted value if all were in the control group
  data |> transform(rx = FALSE) |> 
  posterior_epred(m, newdata = _)

diff |> 
  apply(1, mean) |> # average across observations
  bayestestR::describe_posterior()
#> Summary of Posterior Distribution
#> 
#> Parameter | Median |       95% CI |   pd |          ROPE | % in ROPE
#> --------------------------------------------------------------------
#> Posterior |   0.39 | [0.34, 0.44] | 100% | [-0.10, 0.10] |        0%

# Same result
marginaleffects::marginaleffects(m, variables = "rx") |> 
  summary()
#> Average marginal effects 
#>   Term     Contrast Effect  2.5 % 97.5 %
#> 1   rx TRUE - FALSE 0.3864 0.3354  0.439
#> 
#> Model type:  brmsfit 
#> Prediction type:  response