Think about collapsibility of the odds ratio

[LucyMcGowan]

hmmm @malcolmbarrett this is reminding me of simulations I ran in my doctorate program, but looks like collapsibility is not fixed by IPW and in fact is worse?

library(tidyverse)
library(propensity)
library(survey)

set.seed(1)
n <- 10000

dat <- tibble(
  x_1 = rnorm(n),
  x_2 = rnorm(n),
  x_3 = rnorm(n),
  trt_lin = - 0.5 + 0.25 * x_1 + 0.75 * x_2,
  trt_p = exp(trt_lin) / (1 + exp(trt_lin)),
  treatment = rbinom(n, 1, trt_p),
  y_lin = - 3 + treatment + x_1 + x_2 + 2 * x_3,
  y_p = exp(y_lin) / (1 + exp(y_lin)),
  outcome = rbinom(n, 1, y_p)
)


dat %>%
  mutate(
    p = predict(glm(treatment ~ x_1 + x_2, data = dat, family = binomial()), type = "response"),
    wt = wt_atm(p, treatment)
  ) -> dat
#> ℹ Setting treatment to `1`

## collapsibility is bad (removing x_3 hurts even though not related to treatment)
glm(outcome ~ treatment + x_1 + x_2, data = dat, family = binomial)
#> 
#> Call:  glm(formula = outcome ~ treatment + x_1 + x_2, family = binomial, 
#>     data = dat)
#> 
#> Coefficients:
#> (Intercept)    treatment          x_1          x_2  
#>     -1.9825       0.7342       0.6274       0.6727  
#> 
#> Degrees of Freedom: 9999 Total (i.e. Null);  9996 Residual
#> Null Deviance:       9980 
#> Residual Deviance: 8465  AIC: 8473
## effect should be 1, but is 0.7342

## when we put x_3 in this is fixed
glm(outcome ~ treatment + x_1 + x_2 + x_3, data = dat, family = binomial)
#> 
#> Call:  glm(formula = outcome ~ treatment + x_1 + x_2 + x_3, family = binomial, 
#>     data = dat)
#> 
#> Coefficients:
#> (Intercept)    treatment          x_1          x_2          x_3  
#>     -3.0349       1.0268       0.9537       1.0201       1.9572  
#> 
#> Degrees of Freedom: 9999 Total (i.e. Null);  9995 Residual
#> Null Deviance:       9980 
#> Residual Deviance: 5710  AIC: 5720
## yay!

## let's try to fit the model with weighting
d <- svydesign(
  ~1,
  weights = ~wt,
  data = dat
)
svyglm(outcome ~ treatment, design = d, family = binomial)
#> Warning in eval(family$initialize): non-integer #successes in a binomial glm!
#> Independent Sampling design (with replacement)
#> svydesign(~1, weights = ~wt, data = dat)
#> 
#> Call:  svyglm(formula = outcome ~ treatment, design = d, family = binomial)
#> 
#> Coefficients:
#> (Intercept)    treatment  
#>     -1.6251       0.6821  
#> 
#> Degrees of Freedom: 9999 Total (i.e. Null);  9998 Residual
#> Null Deviance:       10600 
#> Residual Deviance: 10400     AIC: 12060
## uh oh, it's even worse! (0.6821 should be 1)