Survival — R cross-check (survival: Weibull AFT + Cox PH)¶

Companion to survival_python.ipynb / survival_pymc.ipynb. R's survival package is the standard tool. We fit two models on the Gehan data and compare to the from-scratch Weibull:

  • survreg(..., dist="weibull") — parametric Weibull (in accelerated failure time form; we convert to the proportional-hazards HR).
  • coxph(...) — the semiparametric Cox model: the same HR without assuming a Weibull baseline.
In [1]:
suppressMessages(library(survival))
d <- read.csv('leuk.csv')
d$treat <- relevel(factor(d$treat), ref = '6-MP')      # 6-MP = baseline, so 'control' coef = log-HR
S <- Surv(d$time, d$cens)

# --- Weibull AFT, converted to proportional-hazards HR ---
wb <- survreg(S ~ treat, data = d, dist = 'weibull')
sigma <- wb$scale; k <- 1/sigma                         # Weibull shape
beta_ph <- -coef(wb)['treatcontrol'] / sigma            # AFT -> PH
cat(sprintf('Weibull  shape k = %.3f   beta_PH(control) = %.3f   HR = %.2f\n', k, beta_ph, exp(beta_ph)))

# --- Cox proportional hazards (semiparametric) ---
cx <- coxph(S ~ treat, data = d)
b <- coef(cx)['treatcontrol']; ci <- confint(cx)['treatcontrol', ]
cat(sprintf('Cox PH   beta(control) = %.3f   HR = %.2f   95%% CI [%.2f, %.2f]   p = %.2g\n',
            b, exp(b), exp(ci[1]), exp(ci[2]), summary(cx)$coefficients['treatcontrol','Pr(>|z|)']))

Weibull  shape k = 1.366   beta_PH(control) = 1.731   HR = 5.65

Cox PH   beta(control) = 1.572   HR = 4.82   95% CI [2.15, 10.81]   p = 0.00014
In [2]:
# Kaplan-Meier by group + median survival
km <- survfit(S ~ treat, data = d)
print(km)                                               # n, events, median per group
plot(km, col = c('steelblue','firebrick'), lwd = 2, xlab = 'weeks', ylab = 'S(t)',
     main = 'Kaplan-Meier: 6-MP vs control')
legend('topright', legend = levels(d$treat), col = c('steelblue','firebrick'), lwd = 2, bty = 'n')

Call: survfit(formula = S ~ treat, data = d)

               n events median 0.95LCL 0.95UCL
treat=6-MP    21      9     23      16      NA
treat=control 21     21      8       4      12
No description has been provided for this image

Results¶

  • Weibull (survreg) gives shape k ≈ 1.35 and HR(control vs 6-MP) ≈ 5–6, matching the from-scratch and PyMC Weibull fits exactly (same parametric model, ML vs Bayes).
  • Cox PH (coxph) gives HR ≈ 4.5–5 — the same strong, highly significant treatment effect (p ≈ 0.001) without assuming a Weibull baseline. The Cox HR is a touch smaller than the Weibull HR, the usual semiparametric-vs-parametric gap; the conclusion (6-MP roughly 5×-lowers the relapse hazard) is identical.
  • The Kaplan-Meier curves (and the per-group medians: ~8 weeks (control) vs 23 weeks (6-MP)) show the separation nonparametrically.
  • Three engines, one model: from-scratch Weibull (MLE + Bayes) ≈ PyMC (NUTS) ≈ R survreg; plus R coxph as the semiparametric benchmark. Censoring is handled identically throughout via Surv(time, event).