blog/2021/01/15/msm-gee-multilevel/index

[giscus[bot]]

blog/2021/01/15/msm-gee-multilevel/index

Use R to correctly close backdoor confounding in panel data with marginal structural models and inverse probability weights with both GEE and multilevel models

https://www.andrewheiss.com/blog/2021/01/15/msm-gee-multilevel/index.html

[MichaelWillox]

Hi Andrew,

If I hadn't seen your blog, I would have taken a pass on MSMs. Everything is incredibly clear and easy to follow. After playing around with your code for a while and puzzling about why MSM might produce results that differ from a multi-level model, I tried treating the covariate c as a time-varying variable. To do that, I lagged it one period as you did in your data generating function for the treatment t and outcome d variables.

# The only change is that c is included here.
mutate(across(c(t, d, c), list(lag = lag)))

Then, to calculate the IPTW, I added c_lag to the denominator. I also played with dropping c from the numerator as it appears below, but it has almost no effect.

# Automatically with ipwtm()
panel_weights_auto <- ipwtm(
 exposure = t,
 family = "gaussian", 
 numerator = ~ t_lag  , 
 denominator = ~ t_lag + d_lag + c_lag , 
 id = id,
 timevar = y,
 type = "all",
 corstr = "ar1",
 data = as.data.frame(df_panel_sans_na))

Then, as you did, running the MSM produced noticeably better results for the coefficient on t, while the coefficient for t_lag is just slightly better.

> m_auto_gee <- geeglm(d ~ t + t_lag, 
+                      data = df_panel_weights_auto, 
+                      id = id, waves = y, 
+                      weights = ipw)
> 
> summary(m_auto_gee)

Call:
geeglm(formula = d ~ t + t_lag, data = df_panel_weights_auto, 
    weights = ipw, id = id, waves = y)

 Coefficients:
            Estimate  Std.err   Wald Pr(>|W|)    
(Intercept) -0.00144  0.01534   0.01     0.93    
t               0.39413  0.01652 568.98   <2e-16 ***
t_lag        0.21681  0.01548 196.23   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

The results for the multi-level model are about the same for t and slightly better for t_lag.

> m_lmer_auto_gee_wt <- lmer(d ~ t + t_lag + (1 | id), 
+                            data = df_panel_weights_auto, 
+                            weights = ipw)
> summary(m_lmer_auto_gee_wt)
Linear mixed model fit by REML ['lmerMod']
Formula: d ~ t + t_lag + (1 | id)
   Data: df_panel_weights_auto
Weights: ipw

REML criterion at convergence: 18786

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-5.489 -0.573  0.001  0.552  6.813 

Random effects:
 Groups   Name        Variance Std.Dev.
 id       (Intercept) 0.176    0.419   
 Residual             0.374    0.611   
Number of obs: 8000, groups:  id, 2000

Fixed effects:
            Estimate Std. Error t value
(Intercept) -0.00635    0.01211   -0.52
t               0.39063    0.01024   38.14
t_lag        0.18221    0.01037   17.57

I'm curious to know if what I've described here makes sense to you.

Thanks again,
Michael