Causal Inference X(b) — Marginal Structural Models (R companion)¶
ipw::ipwtm — inverse-probability-of-treatment weighting over time, the Robins way¶
R's ipw package (van der Wal & Geskus) is the reference implementation of Robins' marginal structural models. ipwtm builds the time-varying stabilized weights directly from a treatment model; feeding those weights to a weighted pooled regression fits the MSM. This companion reproduces the Python notebook on the same simulated longitudinal survival study with a known effect: the naive analysis is confounded, adjusting for the time-varying confounder is biased toward the null, and the MSM via IPTW recovers the true strong protective effect. Python led with from-scratch IPTW; here ipwtm does the weighting.
1. The treatment-confounder feedback loop¶
Patients are followed over 8 periods. Severity $L_t$ evolves; treatment $A_t$ is more likely when severity is high (confounding), and treatment lowers future severity ($A_{t-1}\to L_t$, the feedback) — which is how it helps, since severity drives the death hazard. So $L_t$ is simultaneously a confounder of $A_t\to Y$ and a mediator of $A_{t-1}\to Y$. We simulate the counterfactuals, so the truth is known: always-treat gives roughly 84% survival versus 45% under never-treat.
suppressMessages({library(ipw); library(survival)})
simulate<-function(n=5000,K=8,cf=NA,seed=0){
set.seed(seed); rows<-list(); alive<-rep(TRUE,n); L<-rnorm(n); Aprev<-rep(0,n); cumA<-rep(0,n); k<-1
for(t in 0:(K-1)){
L<-0.6*L-1.1*Aprev+rnorm(n,0,0.5) # treatment lowers severity (feedback)
pA<-plogis(-0.3+1.3*L) # sicker -> more likely treated (confounding)
A<-if(is.na(cf)) as.numeric(runif(n)<pA) else rep(cf,n)
cumA<-cumA+A
h<-plogis(-2.4+1.1*L-0.15*A) # hazard depends mostly on L
death<-(runif(n)<h)&alive
idx<-which(alive)
rows[[k]]<-data.frame(id=idx,t=t,L=L[idx],A=A[idx],cumA=cumA[idx],death=as.integer(death[idx]),Aprev=Aprev[idx]); k<-k+1
alive<-alive&!death; Aprev<-A; if(sum(alive)==0) break }
do.call(rbind,rows) }
d1<-simulate(cf=1,seed=1); d0<-simulate(cf=0,seed=2)
s1<-1-mean(tapply(d1$death,d1$id,max)); s0<-1-mean(tapply(d0$death,d0$id,max))
cat(sprintf("TRUE counterfactual survival: always-treat %.3f vs never-treat %.3f (difference %+.3f)\n", s1,s0,s1-s0))
cat("L(t) both confounds A(t)->Y and carries the benefit of A(t-1) via A(t-1)->L(t)->Y -- the structure that breaks ordinary regression.\n")
Warning message: "package 'ipw' was built under R version 4.6.1"
TRUE counterfactual survival: always-treat 0.840 vs never-treat 0.453 (difference +0.388)
L(t) both confounds A(t)->Y and carries the benefit of A(t-1) via A(t-1)->L(t)->Y -- the structure that breaks ordinary regression.
2. Both standard analyses fail¶
Pooled logistic regression of death on cumulative treatment, two ways. Naive (death ~ cumA) is confounded — sicker patients accumulate more treatment, so the protective effect is understated. Adjusted (death ~ cumA + L) is biased toward the null — because treatment works through $L$, conditioning on $L$ blocks the benefit. Neither is close to the truth.
d<-simulate(seed=0); d<-d[order(d$id,d$t),]
naive<-coef(glm(death~cumA, d, family=binomial))["cumA"]
adj <-coef(glm(death~cumA+L, d, family=binomial))["cumA"]
cat("Effect of cumulative treatment on log-odds of death (negative = protective):\n")
cat(sprintf(" naive (death ~ cumA) = %+.3f -- confounded: understates the benefit\n", naive))
cat(sprintf(" adjusted (death ~ cumA + L) = %+.3f -- biased toward NULL: L is on the A->L->Y pathway\n", adj))
options(repr.plot.width=7.5, repr.plot.height=4.2)
bp<-barplot(c(naive,adj), names.arg=c("naive\n(death~cumA)","adjusted\n(death~cumA+L)"), col=c("#dd6b20","#c53030"),
ylab="cumulative-treatment coef (log-odds death)", main="Neither standard analysis recovers the strong protective effect")
abline(h=0); text(bp, c(naive,adj)-0.02, sprintf("%+.2f",c(naive,adj)), col="white", font=2)
Effect of cumulative treatment on log-odds of death (negative = protective):
naive (death ~ cumA) = -0.152 -- confounded: understates the benefit
adjusted (death ~ cumA + L) = +0.024 -- biased toward NULL: L is on the A->L->Y pathway
3. The MSM via ipw::ipwtm¶
ipwtm fits the treatment model at each person-period and returns stabilized inverse-probability-of-treatment weights
$$sw_{it}=\prod_{s\le t}\frac{P(A_{is}\mid \bar A_{i,s-1})}{P(A_{is}\mid \bar A_{i,s-1},L_{is})},$$
the numerator omitting the time-varying confounder $L$. Weighting a pooled logistic of death on cumulative treatment by these weights fits the marginal structural model, which builds a pseudo-population where treatment is unconfounded — without conditioning on $L$ — and recovers the true strong protective effect.
w<-ipwtm(exposure=A, family="binomial", link="logit",
numerator=~1, denominator=~L, id=id, timevar=t, type="all", data=as.data.frame(d))
d$sw<-pmin(w$ipw.weights,15)
msm<-coef(glm(death~cumA+t, d, family=binomial, weights=sw))["cumA"]
cat(sprintf("MSM (ipwtm-weighted death ~ cumA): coefficient = %+.3f <-- recovers the strong protective effect\n", msm))
cat(sprintf(" (naive %+.3f, adjusted %+.3f); mean stabilized weight = %.2f (well-behaved)\n", naive, adj, mean(d$sw)))
m<-glm(death~cumA+t, d, family=binomial, weights=sw)
K<-max(d$t); tg<-0:K
surv<-function(strat){ S<-1; out<-numeric(length(tg))
for(i in seq_along(tg)){ cA<-if(strat==1) tg[i]+1 else 0
h<-predict(m, data.frame(cumA=cA,t=tg[i]), type="response"); S<-S*(1-h); out[i]<-S }; out }
S1<-surv(1); S0<-surv(0)
options(repr.plot.width=13, repr.plot.height=4.4); par(mfrow=c(1,2))
barplot(c(naive,adj,msm), names.arg=c("naive","adjusted","MSM-IPTW"), col=c("#dd6b20","#c53030","#2f855a"),
ylab="cumulative-treatment coef", main="Only IPTW recovers the protective effect"); abline(h=0)
plot(tg,S1,type="s",col="#2f855a",lwd=2.5,ylim=c(0,1),xlab="time period",ylab="survival",
main=sprintf("MSM counterfactual survival (true gap %+.2f)",s1-s0))
lines(tg,S0,type="s",col="#c53030",lwd=2.5); abline(h=s1,col="#2f855a",lty=3); abline(h=s0,col="#c53030",lty=3)
legend("bottomleft",c(sprintf("always treat -> %.2f",S1[length(S1)]),sprintf("never treat -> %.2f",S0[length(S0)])),
col=c("#2f855a","#c53030"),lwd=2,bty="n"); par(mfrow=c(1,1))
cat(sprintf("MSM survival gap %+.2f reproduces the true protective effect %+.2f -- by reweighting, not conditioning.\n", S1[length(S1)]-S0[length(S0)], s1-s0))
Warning message in eval(family$initialize): "non-integer #successes in a binomial glm!"
MSM (ipwtm-weighted death ~ cumA): coefficient = -0.354 <-- recovers the strong protective effect
(naive -0.152, adjusted +0.024); mean stabilized weight = 0.94 (well-behaved)
Warning message in eval(family$initialize): "non-integer #successes in a binomial glm!"
MSM survival gap +0.32 reproduces the true protective effect +0.39 -- by reweighting, not conditioning.
4. Summary¶
ipw::ipwtm reproduced Robins' marginal structural model: with a time-varying confounder affected by past treatment, not adjusting left confounding (naive estimate understated a strong protective effect), adjusting introduced new bias (conditioning on the mediator $L$ drove the estimate toward zero), and IPTW — reweighting person-time by the inverse probability of treatment given the past — recovered the true strong protective effect and the counterfactual survival gap, conditioning on nothing.
This is why g-methods exist and why the point-treatment IPW of the causal-survival notebook was only the one-period case: once treatment and a confounder evolve together, only IPTW/MSM, the g-formula, or g-estimation are correct. Practical rule: whenever treatment and a confounder co-evolve over follow-up, do not put the time-varying confounder in an outcome regression — model the treatment process and use stabilized IPTW/MSM. Cross-links: the time-varying generalization of IPW/g-computation (causal survival, subsection 10); the mediator-you-must-not-condition-on is the DAGs collider/mediator lesson over time (subsection 8); the weighting logic is shared with matching/IPW. This completes the depth pass across all ten subsections of the Causal Inference arc.