Causal Inference X — Causal Survival Analysis (R companion)¶

survival, survRM2, and ipw — the hazard-ratio critique made formal, and adjusted RMST¶

R is the native home of survival analysis, and its packages make the causal story crisp:

  • survival — coxph for the hazard ratio and, crucially, cox.zph, a formal test of the proportional-hazards assumption that exposes when a single HR is misleading;
  • survRM2::rmst2 — the reference implementation of the restricted mean survival time comparison;
  • ipw::ipwpoint — inverse-probability-of-treatment weights for confounding adjustment.

On the same simulated data as the Python notebook (true RMST difference ≈ 1.6 years), R confirms every result — and cox.zph adds a formal rejection of proportional hazards, turning "the hazard ratio is treacherous" from an assertion into a test.

1. The hazard ratio, and why cox.zph distrusts it¶

coxph fits the Cox model and reports a single hazard ratio for treatment. But the data are non-proportional (a large early benefit that fades), so that one number is a poor summary. cox.zph tests exactly this: a small p-value means the hazard ratio is not constant over time, so the reported HR is a time-average of a moving target — Hernán's "hazards of hazard ratios," now a formal diagnostic. The HR here is also confounded (sicker patients were treated), which the Kaplan-Meier plot shows.

In [1]:
suppressMessages({library(survival); library(survRM2); library(ipw)})
d<-read.csv("causal_surv_sim.csv"); tau<-8
cx<-coxph(Surv(T,D)~A, d)
cat(sprintf("naive Cox hazard ratio = %.3f\n", exp(coef(cx))))
zph<-cox.zph(cx); cat(sprintf("cox.zph proportional-hazards test: p = %.4f  ", zph$table["A","p"]))
cat(if(zph$table["A","p"]<0.05) "-> REJECT PH: the hazard ratio varies over time (single HR misleads)\n" else "\n")
options(repr.plot.width=8.5, repr.plot.height=4.6)
sf<-survfit(Surv(T,D)~A, d)
plot(sf, col=c("#2b6cb0","#c53030"), lwd=2.5, xlab="time t (years)", ylab="survival S(t)", xlim=c(0,tau),
     main="Naive Kaplan-Meier by treatment (confounded)")
legend("topright", c("untreated","treated"), col=c("#2b6cb0","#c53030"), lwd=2.5, bty="n")

Warning message:
"package 'survRM2' was built under R version 4.6.1"

Warning message:
"package 'ipw' was built under R version 4.6.1"

naive Cox hazard ratio = 0.601

cox.zph proportional-hazards test: p = 0.0003  
-> REJECT PH: the hazard ratio varies over time (single HR misleads)
No description has been provided for this image

2. Naive vs IPTW-adjusted RMST — survRM2 and ipw¶

survRM2::rmst2 gives the unadjusted RMST difference between arms — confounded here, understating the benefit because the treated were sicker. ipw::ipwpoint then estimates stabilized inverse-probability-of-treatment weights (numerator $P(A)$, denominator $P(A\mid L)$), and a weighted Kaplan-Meier in each arm yields the confounding-adjusted survival curves. Integrating them recovers the true RMST difference of ~1.6 years — the adjusted answer survRM2 alone cannot produce, supplied by the IPW weights.

In [2]:
r0<-rmst2(d$T, d$D, d$A, tau=tau)
naive_rmst<-r0$unadjusted.result[1,1]
w<-ipwpoint(exposure=A, family="binomial", link="logit", numerator=~1, denominator=~L, data=d)
d$sw<-w$ipw.weights
wrmst<-function(sub){ s<-survfit(Surv(T,D)~1, data=sub, weights=sub$sw)
  tt<-c(0,s$time); ss<-c(1,s$surv); keep<-tt<=tau; tt<-c(tt[keep],tau); ss<-ss[1:length(tt)]; sum(diff(tt)*head(ss,-1)) }
ipw_rmst<-wrmst(d[d$A==1,])-wrmst(d[d$A==0,])
cat(sprintf("survRM2 naive RMST difference = %.3f  (confounded)\n", naive_rmst))
cat(sprintf("IPTW-adjusted RMST difference  = %.3f  (true ~1.62)\n", ipw_rmst))
options(repr.plot.width=8.5, repr.plot.height=4.6)
grid<-seq(0,tau,length.out=200); stepf<-function(s,g){ i<-findInterval(g,c(0,s$time)); c(1,s$surv)[i+1] }
s1<-survfit(Surv(T,D)~1,data=d[d$A==1,],weights=d$sw[d$A==1]); s0<-survfit(Surv(T,D)~1,data=d[d$A==0,],weights=d$sw[d$A==0])
plot(grid, stepf(s1,grid), type="s", col="#c53030", lwd=2.5, ylim=c(0,1), xlab="time t", ylab="survival S(t)",
     main=sprintf("IPTW-adjusted survival curves (RMST diff %.2f)", ipw_rmst))
lines(grid, stepf(s0,grid), type="s", col="#2b6cb0", lwd=2.5)
legend("topright", c("treated (adjusted)","untreated (adjusted)"), col=c("#c53030","#2b6cb0"), lwd=2.5, bty="n")

survRM2 naive RMST difference = 0.974  (confounded)

IPTW-adjusted RMST difference  = 1.618  (true ~1.62)
No description has been provided for this image

3. G-computation — standardization with survreg¶

The outcome-model route: fit a Weibull AFT with survreg(Surv(T,D) ~ A + L), then standardize — predict every patient's survival curve setting $A=1$ and averaging over the observed $L$ to form $\hat S_1(t)$, likewise $\hat S_0(t)$ — and take the RMST difference. This is the survival g-formula; it recovers the same ~1.6-year effect through the outcome model that IPW reached through the treatment model. Reporting both is the honest robustness check (and combining them is doubly robust).

In [3]:
fit<-survreg(Surv(T,D)~A+L, data=d, dist="weibull")
scale<-fit$scale; grid<-seq(0.05,tau,length.out=200)
Sfun<-function(a){ lp<-predict(fit, newdata=transform(d,A=a), type="lp")   # log-scale linear predictor
  sapply(grid, function(t) mean(exp(-(t/exp(lp))^(1/scale)))) }            # avg Weibull survival over L
S1<-Sfun(1); S0<-Sfun(0)
gcomp<-sum(diff(grid)*head(S1,-1))-sum(diff(grid)*head(S0,-1))
cat(sprintf("g-computation (survreg standardization) RMST difference = %.3f  (true ~1.62)\n", gcomp))
options(repr.plot.width=7.5, repr.plot.height=4.4)
plot(grid,S1,type="l",col="#c53030",lwd=2.5,ylim=c(0,1),xlab="time t",ylab="survival",main="Standardized counterfactual survival curves (g-computation)")
lines(grid,S0,col="#2b6cb0",lwd=2.5); legend("topright",c("S1(t) treat-all","S0(t) treat-none"),col=c("#c53030","#2b6cb0"),lwd=2.5,bty="n")

g-computation (survreg standardization) RMST difference = 1.565  (true ~1.62)
No description has been provided for this image

4. Summary¶

R's survival packages reproduced the causal analysis and sharpened its central message: cox.zph formally rejected proportional hazards (p ≈ 0.0003), proving the single Cox hazard ratio is a time-average of a changing effect; survRM2 gave the interpretable (confounded) RMST difference; ipw supplied the weights that adjusted it back to the truth; and survreg standardization (g-computation) reached the same ~1.6-year effect through the outcome model. IPW and g-computation agreeing is the two-model robustness check.

The workflow to carry: prefer RMST or survival-probability differences to the hazard ratio (and test PH with cox.zph before ever quoting an HR); adjust for confounding with IPTW (ipw) or standardization (survreg g-formula); and reach for marginal structural models when treatment and confounders both vary over time. This closes the Causal Inference arc — the same survival/survRM2 tools of the Survival catalog, now answering a causal question. Cross-links: IPTW ↔ the matching/weighting of subsection 2; g-computation ↔ back-door standardization of the DAG notebook; doubly-robust survival ↔ the AIPW/DML of subsection 9.