Causal Inference X — Causal Survival Analysis¶
Treatment effects on time-to-event outcomes: IPW, g-computation, RMST, and the hazards of hazard ratios¶
The arc closes where causal inference meets survival analysis — the setting of most clinical and many economic questions, where the outcome is a time to an event (death, relapse, unemployment exit, default) and the data are censored: for many units we only know they had not yet had the event when observation stopped. Everything from the identification subsections still applies — we need to defeat confounding — but two survival-specific complications appear, and one deeply-rooted habit needs to be unlearned.
The estimands are counterfactual survival curves $S_a(t)=P(T^a>t)$: the probability of surviving past $t$ if everyone were assigned treatment $a$. From them come the honest causal contrasts — the survival-probability difference $S_1(t)-S_0(t)$ and the restricted mean survival time (RMST) difference, the difference in average event-free time up to a horizon $\tau$. Notably absent from that list is the quantity almost every survival paper reports: the hazard ratio.
This notebook makes four points, each on a simulation with a known counterfactual truth:
- The hazards of hazard ratios (Hernán 2010) — the HR is non-collapsible and generally time-varying, so a single Cox HR is a misleading summary of a causal effect even in a clean experiment;
- Inverse-probability weighting (IPTW) — reweight by the propensity of treatment to build confounding-adjusted survival curves;
- G-computation (standardization) — fit an outcome model, predict everyone's survival under each treatment, and average;
- RMST as the interpretable, collapsible causal summary — and a note on marginal structural models for time-varying confounding (Robins).
Python-lead (from-scratch Kaplan-Meier / IPW / g-computation + lifelines); R companion uses survival, survRM2, and ipw. This is the final example of the Causal Inference arc, and it cross-links directly to the Survival catalog (Weibull PH, Cox-via-Poisson, frailty, interval-censored).
1. The data, and why the hazard ratio misleads¶
We simulate an observational study of $n$ patients. A confounder $L$ (disease severity) drives both treatment ($A$: sicker patients are more likely treated) and survival (sicker patients die sooner) — classic confounding by indication. Survival times follow a log-normal model, $\log T^a=\mu(L)+\gamma\,a+\sigma W$ with $\gamma>0$ (treatment prolongs survival), and observation is censored administratively at $\tau$ and by random dropout. Because we generate both potential times $T^0,T^1$ for everyone, we know the truth: the RMST difference is about 1.6 years.
Two problems surface immediately. First, confounding: the naive Kaplan-Meier comparison of treated vs untreated understates the benefit, because the treated were sicker to begin with. Second, and more subtly, the hazard ratio is a treacherous summary. In this (realistic, non-proportional-hazards) data the true ratio of hazards $h_1(t)/h_0(t)$ changes with time — a large early benefit that fades — so no single number describes it; yet a Cox model forces exactly one. The HR is also non-collapsible (the marginal HR differs from the covariate-conditional HR even with no confounding), so it is not a clean causal contrast. This is Hernán's "hazards of hazard ratios".
import numpy as np, pandas as pd, matplotlib.pyplot as plt, warnings
warnings.filterwarnings("ignore")
from lifelines import CoxPHFitter
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; ORANGE="#dd6b20"; GREY="#a0aec0"; PURP="#6b46c1"
rng=np.random.default_rng(0); n=4000; tau=8.0
L=rng.uniform(0,1,n) # severity confounder
e=1/(1+np.exp(-(-0.3+2.2*(L-0.5)))) # P(treat|L): sicker -> more likely treated
A=(rng.uniform(size=n)<e).astype(int)
mu=1.6-1.2*L; gamma=0.5; sig=0.5; W=rng.normal(0,1,n)
T0=np.exp(mu+sig*W); T1=np.exp(mu+gamma+sig*W) # counterfactual survival times (log-normal AFT)
rmst_true=np.mean(np.minimum(T1,tau))-np.mean(np.minimum(T0,tau))
TA=np.where(A==1,T1,T0)
C=np.minimum(rng.exponential(12,n),tau) # censoring + admin cutoff at tau
T=np.clip(np.minimum(TA,C),1e-3,None); Dr=(TA<=C).astype(int)
d=pd.DataFrame(dict(T=T,D=Dr,A=A,L=L))
print(f"n={n}, treated fraction={A.mean():.2f}, event rate={Dr.mean():.2f}, horizon tau={tau}")
print(f"TRUE RMST difference (treatment − control) = {rmst_true:.2f} years")
# naive Cox HR
hr=np.exp(CoxPHFitter().fit(d[['T','D','A']],'T','D').params_['A'])
hr_adj=np.exp(CoxPHFitter().fit(d[['T','D','A','L']],'T','D').params_['A']) # same model, L included
# true time-varying hazard ratio from the counterfactual populations
def hz(Tx,t,h=0.3):
S=np.mean(Tx>t); f=np.mean((Tx>t)&(Tx<=t+h))/h; return f/max(S,1e-9)
tg=np.linspace(0.5,6,25); hrt=[hz(T1,t)/hz(T0,t) for t in tg]
fig,ax=plt.subplots(1,2,figsize=(13,4.4))
ax[0].hist(L[A==1],bins=25,alpha=.55,color=RED,density=True,label="treated"); ax[0].hist(L[A==0],bins=25,alpha=.55,color=BLUE,density=True,label="untreated")
ax[0].set_xlabel("severity L"); ax[0].set_ylabel("density"); ax[0].set_title("Confounding by indication: the treated are sicker"); ax[0].legend()
ax[1].plot(tg,hrt,"o-",color=PURP,lw=2,label="true hazard ratio h₁(t)/h₀(t)"); ax[1].axhline(hr,color=RED,ls="--",label=f"single Cox HR = {hr:.2f}")
ax[1].axhline(1,color="k",lw=.6); ax[1].set_xlabel("time t"); ax[1].set_ylabel("hazard ratio"); ax[1].set_title("The true HR varies over time — one Cox number can't capture it"); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.show()
print(f"The Cox model reports a single HR = {hr:.2f}, but the true hazard ratio sweeps from ~{hrt[0]:.2f} (large early benefit)")
print(f"to ~{hrt[-1]:.2f} (benefit nearly gone). A time-averaged HR hides that story -- Hernán's 'hazards of hazard ratios'.")
print()
print("SEPARATE THE TWO COMPLAINTS, THOUGH, BECAUSE THAT HR HAS *TWO* THINGS WRONG WITH IT.")
print(f" unadjusted Cox HR {hr:.3f} the number above")
print(f" the same Cox, L included {hr_adj:.3f} confounding removed")
print(f" true HR over follow-up {min(hrt):.3f} to {max(hrt):.3f}")
print()
print(f"Confounding accounts for {hr-hr_adj:.3f} of the distance -- treated patients were sicker, so the raw")
print("comparison understates the benefit, exactly as the naive RMST does later. That part is fixable by")
print("adjustment, and adjusting fixes it.")
print()
print(f"What adjustment does NOT fix is the rest. The L-adjusted {hr_adj:.2f} is a perfectly respectable")
print(f"confounding-free number, and it is still a single constant standing in for something that runs")
print(f"from {min(hrt):.2f} to {max(hrt):.2f}. It lands somewhere inside that range with nothing to say about where, and")
print("no diagnostic on the adjusted model reveals the problem, because the problem is not bias in the")
print("usual sense -- it is that the estimand itself is a summary of a curve that changed shape.")
print("That is why the rest of this notebook estimates RMST instead of arguing about which HR to report.")
n=4000, treated fraction=0.43, event rate=0.69, horizon tau=8.0 TRUE RMST difference (treatment − control) = 1.62 years
The Cox model reports a single HR = 0.60, but the true hazard ratio sweeps from ~0.14 (large early benefit) to ~0.75 (benefit nearly gone). A time-averaged HR hides that story -- Hernán's 'hazards of hazard ratios'. SEPARATE THE TWO COMPLAINTS, THOUGH, BECAUSE THAT HR HAS *TWO* THINGS WRONG WITH IT. unadjusted Cox HR 0.601 the number above the same Cox, L included 0.352 confounding removed true HR over follow-up 0.139 to 0.750 Confounding accounts for 0.249 of the distance -- treated patients were sicker, so the raw comparison understates the benefit, exactly as the naive RMST does later. That part is fixable by adjustment, and adjusting fixes it. What adjustment does NOT fix is the rest. The L-adjusted 0.35 is a perfectly respectable confounding-free number, and it is still a single constant standing in for something that runs from 0.14 to 0.75. It lands somewhere inside that range with nothing to say about where, and no diagnostic on the adjusted model reveals the problem, because the problem is not bias in the usual sense -- it is that the estimand itself is a summary of a curve that changed shape. That is why the rest of this notebook estimates RMST instead of arguing about which HR to report.
2. Inverse-probability weighting — confounding-adjusted survival curves¶
The first job is to remove the confounding. Inverse-probability-of-treatment weighting (IPTW) reweights each unit by the inverse of its probability of receiving the treatment it got, creating a pseudo-population in which treatment is independent of the confounders $L$. We use stabilized weights $sw_i=\frac{P(A=a_i)}{P(A=a_i\mid L_i)}$ (lower variance than raw weights), then estimate a weighted Kaplan-Meier curve in each arm. From scratch, the weighted KM simply replaces counts of at-risk and events with sums of weights.
The adjusted curves separate correctly: the naive KM understated the treatment benefit (treated patients were sicker), and IPTW restores it. Reading the RMST difference (area between the curves up to $\tau$) off the weighted curves recovers the true ~1.6 years, while the naive curves badly under-report it.
def km(T,D,w):
et=np.unique(T[D==1]); S=1.0; ts=[0.0]; sv=[1.0]
for t in et:
atrisk=w[T>=t].sum(); dead=w[(T==t)&(D==1)].sum(); S*=(1-dead/atrisk); ts.append(t); sv.append(S)
return np.array(ts),np.array(sv)
def step(ts,sv,grid): return sv[np.searchsorted(ts,grid,side="right")-1]
def rmst(T,D,w,tau):
ts,sv=km(T,D,w); ts=np.append(ts,tau); sv=np.append(sv,sv[-1]); m=ts<=tau
return np.trapezoid(sv[m],ts[m])
w1=np.ones(n); sw=A*(A.mean()/e)+(1-A)*((1-A.mean())/(1-e)) # stabilized IPTW weights
grid=np.linspace(0,tau,200)
naive_diff=rmst(T[A==1],Dr[A==1],w1[A==1],tau)-rmst(T[A==0],Dr[A==0],w1[A==0],tau)
ipw_diff =rmst(T[A==1],Dr[A==1],sw[A==1],tau)-rmst(T[A==0],Dr[A==0],sw[A==0],tau)
fig,ax=plt.subplots(1,2,figsize=(13,4.6))
for lbl,w,a in [("naive",w1,ax[0]),("IPTW-adjusted",sw,ax[1])]:
for arm,c,nm in [(1,RED,"treated"),(0,BLUE,"untreated")]:
ts,sv=km(T[A==arm],Dr[A==arm],w[A==arm]); a.step(grid,step(ts,sv,grid),where="post",color=c,lw=2.2,label=nm)
a.set_xlabel("time t (years)"); a.set_ylabel("survival S(t)"); a.set_ylim(0,1); a.legend()
a.set_title(f"{lbl} Kaplan-Meier (RMST diff {rmst(T[A==1],Dr[A==1],w[A==1],tau)-rmst(T[A==0],Dr[A==0],w[A==0],tau):.2f})")
plt.tight_layout(); plt.show()
print(f"TRUE RMST difference = {rmst_true:.2f} years")
print(f"naive KM RMST diff = {naive_diff:.2f} (confounded: treated were sicker -> benefit understated)")
print(f"IPTW RMST diff = {ipw_diff:.2f} (pseudo-population removes confounding -> recovers the truth)")
TRUE RMST difference = 1.62 years naive KM RMST diff = 0.97 (confounded: treated were sicker -> benefit understated) IPTW RMST diff = 1.69 (pseudo-population removes confounding -> recovers the truth)
3. G-computation (standardization) — the outcome-model route¶
IPW models the treatment; g-computation models the outcome. We fit a survival regression $S(t\mid A,L)$ (a Weibull/log-normal AFT), then form each counterfactual survival curve by standardization — predict every patient's survival curve setting $A=1$, average over the observed distribution of $L$ to get $\hat S_1(t)$; repeat with $A=0$ for $\hat S_0(t)$:
$$\hat S_a(t)=\frac1n\sum_{i=1}^{n}\hat S\big(t\mid A=a,\,L_i\big).$$
This is the survival version of the g-formula (Robins). Like IPW it recovers the true RMST difference, and the two are complementary — IPW is consistent if the treatment model is right, g-computation if the outcome model is right (combining them gives a doubly-robust estimator, the survival analogue of the AIPW/DML of earlier notebooks). We use lifelines' AFT fitter for the outcome model and confirm the standardized curves against the truth.
from lifelines import WeibullAFTFitter
aft=WeibullAFTFitter().fit(d[['T','D','A','L']],'T','D')
tg=np.linspace(0.05,tau,200)
d1=d.copy(); d1['A']=1; d0=d.copy(); d0['A']=0
S1=aft.predict_survival_function(d1,times=tg).mean(axis=1).values # standardized S_1(t)
S0=aft.predict_survival_function(d0,times=tg).mean(axis=1).values # standardized S_0(t)
gcomp_diff=np.trapezoid(S1,tg)-np.trapezoid(S0,tg)
# truth curves for overlay
S1_true=np.array([np.mean(T1>t) for t in tg]); S0_true=np.array([np.mean(T0>t) for t in tg])
fig,ax=plt.subplots(1,2,figsize=(13,4.6))
ax[0].plot(tg,S1,color=RED,lw=2.5,label="g-comp $\\hat S_1(t)$"); ax[0].plot(tg,S0,color=BLUE,lw=2.5,label="g-comp $\\hat S_0(t)$")
ax[0].plot(tg,S1_true,color=RED,lw=1,ls=":",label="true $S_1$"); ax[0].plot(tg,S0_true,color=BLUE,lw=1,ls=":",label="true $S_0$")
ax[0].set_xlabel("time t"); ax[0].set_ylabel("survival"); ax[0].set_ylim(0,1); ax[0].set_title("Standardized counterfactual survival curves"); ax[0].legend(fontsize=8)
ax[1].bar(["true","naive KM","IPTW","g-comp"],[rmst_true,naive_diff,ipw_diff,gcomp_diff],color=[GREY,ORANGE,GREEN,BLUE])
ax[1].axhline(rmst_true,color="k",ls="--");
for i,v in enumerate([rmst_true,naive_diff,ipw_diff,gcomp_diff]): ax[1].text(i,v+0.03,f"{v:.2f}",ha="center")
ax[1].set_ylabel("RMST difference (years)"); ax[1].set_title("Both adjustment methods recover the causal RMST")
plt.tight_layout(); plt.show()
print(f"g-computation RMST difference = {gcomp_diff:.2f} (true {rmst_true:.2f}). The standardized curves track the truth,")
print("and IPW (treatment model) and g-computation (outcome model) agree -- two routes to the same adjusted effect.")
g-computation RMST difference = 1.57 (true 1.62). The standardized curves track the truth, and IPW (treatment model) and g-computation (outcome model) agree -- two routes to the same adjusted effect.
4. RMST as the estimand, and marginal structural models¶
Putting the summaries side by side makes the case for RMST. The hazard ratio compressed a time-varying, non-collapsible effect into one slippery number. The RMST difference — "treatment adds about 1.6 event-free years up to the 8-year horizon" — is a collapsible, absolute, directly interpretable causal quantity on the time scale patients and policymakers care about, and it is exactly what IPW and g-computation delivered. The survival-probability difference at a chosen time is its equally-honest pointwise cousin. Modern practice (and regulators) increasingly favour these over the reflexive Cox HR, precisely for the reasons Hernán laid out.
One frontier remains, and it is where survival g-methods truly earn their keep: time-varying treatment with time-varying confounding. When a confounder is affected by past treatment and also predicts future treatment and the outcome (e.g., a lab value that responds to a drug and guides the next dose), standard regression adjustment is biased either way — adjust for it and you block part of the effect and open a collider; don't and confounding remains. Marginal structural models (Robins 2000), fit by inverse-probability-of-treatment weighting over time, are the canonical solution: the same IPW logic used above, applied sequentially, is the only one of these tools that handles it. Our point-treatment example is the one-period special case; the full machinery is the reason IPW/g-methods are indispensable in longitudinal causal inference.
# summary comparison of estimands / methods
summ=pd.DataFrame({
"estimate":[f"{hr:.2f}","—",f"{naive_diff:.2f}",f"{ipw_diff:.2f}",f"{gcomp_diff:.2f}",f"{rmst_true:.2f}"],
"note":["single number for a time-varying, non-collapsible effect",
"(hazard ratio is not a clean causal contrast)",
"confounded (treated were sicker)",
"confounding removed via treatment model",
"confounding removed via outcome model",
"the known truth"]},
index=["Cox hazard ratio","—","naive RMST diff","IPTW RMST diff","g-comp RMST diff","TRUE RMST diff"])
print(summ.to_string())
# survival-probability difference at t=4, adjusted
ts1,sv1=km(T[A==1],Dr[A==1],sw[A==1]); ts0,sv0=km(T[A==0],Dr[A==0],sw[A==0])
sd=step(ts1,sv1,4)-step(ts0,sv0,4); sd_true=np.mean(T1>4)-np.mean(T0>4)
print(f"\nIPTW survival-probability difference at t=4: {sd:+.3f} (true {sd_true:+.3f})")
print("RMST and survival-difference are the interpretable, collapsible causal contrasts; the hazard ratio is not.")
estimate note Cox hazard ratio 0.60 single number for a time-varying, non-collapsible effect — — (hazard ratio is not a clean causal contrast) naive RMST diff 0.97 confounded (treated were sicker) IPTW RMST diff 1.69 confounding removed via treatment model g-comp RMST diff 1.57 confounding removed via outcome model TRUE RMST diff 1.62 the known truth IPTW survival-probability difference at t=4: +0.323 (true +0.317) RMST and survival-difference are the interpretable, collapsible causal contrasts; the hazard ratio is not.
5. Summary — and the close of the arc¶
Causal survival analysis brings the whole arc to time-to-event outcomes under censoring. On a confounded simulation with a known truth of ~1.6 added event-free years:
- the hazard ratio — the field's default — proved a poor causal summary: a single Cox HR papered over a genuinely time-varying effect (large early, fading later) and is non-collapsible, exactly Hernán's warning;
- inverse-probability weighting built a confounding-free pseudo-population and adjusted survival curves that recovered the truth;
- g-computation reached the same answer through an outcome model and standardization (the survival g-formula);
- RMST and the survival-probability difference are the interpretable, collapsible causal estimands to report — and marginal structural models extend the IPW logic to the time-varying confounding that defeats ordinary adjustment.
Cross-links. This notebook is the causal face of the Survival catalog — the same Weibull/Cox/Kaplan-Meier machinery (from the Weibull-PH, Cox-via-Poisson, frailty, and interval-censored notebooks), now asked to answer a causal question rather than a descriptive one. The IPW here is the survival version of the matching/weighting of subsection 2; g-computation is the standardization / back-door adjustment of the DAG notebook; the doubly-robust combination is the AIPW/DML idea of subsection 9 carried to censored data; and the untestable assumption underneath it all is again unconfoundedness (now "no unmeasured confounding, plus non-informative censoring").
With this, the Causal Inference arc is complete across all ten subsections — from Fisher's randomized experiments, through the observational identification toolkit (matching, IV, RD, panel, difference-in-differences, synthetic control), to the graphical language of the SCM, machine-learning estimation of heterogeneous and debiased effects, and finally causal inference for censored survival. Two languages — potential outcomes and structural graphs — one disciplined question throughout: what would have happened otherwise?