Survival analysis — Weibull proportional hazards with censoring

From-scratch MLE + Bayes; the BUGS "Leuk" leukemia example

Module: survival_mcmc.py. Time-to-event data (here: weeks to leukemia relapse) bring a genuinely new ingredient — censoring. The Gehan data (BUGS Vol I "Leuk"; PyMC's Bayesian-survival case study) compares the drug 6-MP against placebo in 42 leukemia patients; many 6-MP patients were still in remission when the study ended, so their true relapse time is only known to exceed the observed time.

---

Model & the censored likelihood

A subject with an observed event contributes the density $f(t)=h(t)S(t)$; a censored subject (still event-free at $t$) contributes only the survival $S(t)$. With event indicator $\delta_i$ (1 = relapse seen, 0 = censored), $$L=\prod_i h(t_i)^{\delta_i}\,S(t_i).$$ Weibull proportional hazards: $h(t\mid x)=\lambda_i\,k\,t^{k-1}$, $\lambda_i=e^{x_i'\beta}$, $S(t)=e^{-\lambda_i t^k}$, so $$\log L=\sum_i \delta_i\big[x_i'\beta+\log k+(k-1)\log t_i\big]-\sum_i e^{x_i'\beta}t_i^{k}.$$

• $e^{\beta_j}$ is the hazard ratio for covariate $j$.
• $k$ = Weibull shape: $k>1$ hazard rises over time, $k<1$ falls, $k=1$ = exponential (constant hazard).
• Note this is a Poisson GLM kernel in disguise ($\delta$ as 0/1 "counts", offset $\log t^k$) — the piecewise-exponential = Poisson trick behind BUGS Leuk, linking survival back to the count-model arc.

Algorithm

• MLE: maximise $\log L$ over $(\beta,\log k)$ (Nelder-Mead), SEs from the numerical Hessian.
• Bayes: RW-Metropolis on $(\beta,\log k)$ with a full-covariance proposal from the MLE Hessian; priors $\beta\sim N(0,10^2)$, $\log k\sim N(0,1)$.
• Kaplan-Meier (km_estimator) gives the nonparametric survival curve for comparison.

import numpy as np, pandas as pd, matplotlib.pyplot as plt
from survival_mcmc import simulate_weibull, weibull_mle, weibull_rwm, km_estimator, summary

# --- Section 1: synthetic recovery (Weibull with right-censoring) ---
t, X, delta = simulate_weibull(n=2000, beta=(-1.0, -1.3), k=1.4, seed=1)
print('synthetic: n=%d  censored=%.0f%%' % (len(t), 100*(1-delta.mean())))
mle, se, _ = weibull_mle(X, t, delta)
print('  MLE   beta0 %.3f(%.3f)  beta1 %.3f(%.3f)  k %.3f   (truth -1.0, -1.3, k=1.4)' % (mle[0],se[0],mle[1],se[1],np.exp(mle[2])))
o = weibull_rwm(X, t, delta, seed=1)
for nm,m,s,lo,hi in summary(o['draws'],['beta0','beta1','logk']): print('  Bayes %-6s %.3f (%.3f)'%(nm,m,s))
print('  Bayes HR(treat)=exp(beta1)= %.3f   acceptance %.2f' % (np.exp(o['draws'][:,1].mean()), o['accept']))

synthetic: n=2000  censored=55%
  MLE   beta0 -1.059(0.049)  beta1 -1.380(0.073)  k 1.433   (truth -1.0, -1.3, k=1.4)

  Bayes beta0  -1.061 (0.050)
  Bayes beta1  -1.384 (0.073)
  Bayes logk   0.360 (0.024)
  Bayes HR(treat)=exp(beta1)= 0.251   acceptance 0.46

---

Section 2 — the Gehan leukemia data (BUGS "Leuk")

42 patients, randomised to 6-MP (mercaptopurine) or placebo, time = weeks to leukemia relapse. The placebo arm relapsed fast and fully (all 21 events, median 8 weeks); the 6-MP arm did far better — 12 of 21 were still in remission at the end (censored), median 16 weeks. We code treat so $\beta>0$ means higher hazard for the control group.

d = pd.read_csv('leuk.csv')
trt = (d['treat'] == 'control').astype(float).to_numpy()      # 1 = control (placebo), 0 = 6-MP
t = d['time'].to_numpy(float); delta = d['cens'].to_numpy(float)
X = np.column_stack([np.ones(len(d)), trt])
for g,lab in [(0,'6-MP'),(1,'control')]:
    m = trt==g; print('  %-8s n=%d  events=%d  censored=%d  median.obs=%.0f' % (lab,m.sum(),int(delta[m].sum()),int((1-delta[m]).sum()),np.median(t[m])))

mle, se, _ = weibull_mle(X, t, delta)
exp_mle, _, _ = weibull_mle(X, t, delta, fix_k=True)           # exponential (k=1) for comparison
o = weibull_rwm(X, t, delta, seed=1); dr = o['draws']
hr = np.exp(dr[:,1]); k_post = np.exp(dr[:,2])
print('\nWeibull MLE : beta_control %.3f (%.3f)  k %.3f' % (mle[1], se[1], np.exp(mle[2])))
print('Weibull Bayes: beta_control %.3f [%.2f, %.2f]' % (dr[:,1].mean(), np.percentile(dr[:,1],2.5), np.percentile(dr[:,1],97.5)))
print('  HR (control vs 6-MP) = %.2f [%.2f, %.2f]   => 6-MP HR = %.3f   P(HR>1) = %.3f' % (hr.mean(), np.percentile(hr,2.5), np.percentile(hr,97.5), 1/hr.mean(), (hr>1).mean()))
print('  shape k = %.3f [%.2f, %.2f]  (k>1 => hazard rises with time; CI vs 1 tells Weibull-vs-exponential)' % (k_post.mean(), np.percentile(k_post,2.5), np.percentile(k_post,97.5)))

  6-MP     n=21  events=9  censored=12  median.obs=16
  control  n=21  events=21  censored=0  median.obs=8

Weibull MLE : beta_control 1.731 (0.413)  k 1.366
Weibull Bayes: beta_control 1.745 [0.92, 2.64]
  HR (control vs 6-MP) = 6.31 [2.51, 14.05]   => 6-MP HR = 0.158   P(HR>1) = 1.000
  shape k = 1.345 [0.97, 1.74]  (k>1 => hazard rises with time; CI vs 1 tells Weibull-vs-exponential)

# survival curves: Kaplan-Meier (nonparametric) + fitted Weibull, by group
fig, ax = plt.subplots(1, 2, figsize=(12, 4.6))
bm = mle; kk = np.exp(bm[2]); tt = np.linspace(0, 36, 250)
for g, lab, col in [(0,'6-MP (treated)','steelblue'), (1,'control (placebo)','firebrick')]:
    m = trt==g
    tk, Sk = km_estimator(t[m], delta[m])
    ax[0].step(tk, Sk, where='post', color=col, lw=1.8, label='KM '+lab)
    lam = np.exp(bm[0] + bm[1]*g); ax[0].plot(tt, np.exp(-lam*tt**kk), '--', color=col, lw=1.2)
    ct = t[m & (delta==0)]                                   # censoring marks at KM height
    if len(ct):
        idx = np.searchsorted(tk, ct, side='right')-1
        ax[0].plot(ct, Sk[np.clip(idx,0,len(Sk)-1)], '+', color=col, ms=9, mew=1.5)
ax[0].set_xlabel('weeks'); ax[0].set_ylabel('S(t) = P(still in remission)'); ax[0].set_ylim(0,1.02)
ax[0].set_title('Kaplan-Meier (step) vs fitted Weibull (dashed); + = censored'); ax[0].legend(fontsize=8)
ax[1].hist(hr, bins=50, density=True, color='seagreen', alpha=.8)
ax[1].axvline(hr.mean(), color='k', ls='--', lw=1, label=f'mean {hr.mean():.1f}')
ax[1].axvline(1, color='firebrick', ls=':', lw=1, label='no effect (HR=1)')
ax[1].set_xlabel('hazard ratio  control vs 6-MP'); ax[1].set_yticks([])
ax[1].set_title(f'Treatment effect  (P(HR>1) = {(hr>1).mean():.3f})'); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.savefig('survival_leuk.png', dpi=120, bbox_inches='tight'); plt.show()

The hazard rate

The survival curve answers "what fraction are still in remission?"; the hazard $h(t)=f(t)/S(t)$ answers a sharper question — "for someone still in remission at week $t$, how fast are they relapsing right now?" — the instantaneous risk (a rate, not a probability). For the Weibull, $h(t\mid x)=\lambda\,k\,t^{k-1}$.

• Left — fitted hazard by group. Both curves rise with time (since $k\approx1.35>1$: the longer in remission, the higher the instantaneous relapse risk), and the control hazard sits a constant ~5.7× above 6-MP at every $t$. That fixed vertical gap is the proportional-hazards assumption — and the hazard ratio.
• Right — cumulative hazard $H(t)=-\log S(t)=\lambda t^k$. The nonparametric Nelson-Aalen estimate (step) tracks the fitted Weibull (dashed): a data-vs-fit check on the hazard scale. (We plot the cumulative hazard rather than the raw instantaneous empirical hazard, which is too noisy to estimate from 42 patients without smoothing.)

# the hazard view: fitted hazard h(t) by group + cumulative hazard (Nelson-Aalen vs fitted)
def nelson_aalen(tt_, dd_):
    o = np.argsort(tt_); tt_, dd_ = tt_[o], dd_[o]
    H = 0.0; ts = [0.0]; Hs = [0.0]
    for ti in np.unique(tt_[dd_ == 1]):
        H += np.sum((tt_ == ti) & (dd_ == 1)) / np.sum(tt_ >= ti); ts.append(ti); Hs.append(H)
    return np.array(ts), np.array(Hs)

fig, ax = plt.subplots(1, 2, figsize=(12, 4.6))
grid = np.linspace(0.1, 36, 250); kk = np.exp(bm[2])
for g, lab, col in [(0,'6-MP (treated)','steelblue'), (1,'control (placebo)','firebrick')]:
    lam = np.exp(bm[0] + bm[1]*g)
    ax[0].plot(grid, lam*kk*grid**(kk-1), color=col, lw=1.9, label=lab)         # fitted hazard
    m = trt==g; ta, Ha = nelson_aalen(t[m], delta[m])
    ax[1].step(ta, Ha, where='post', color=col, lw=1.6, label='Nelson-Aalen '+lab)
    ax[1].plot(grid, lam*grid**kk, '--', color=col, lw=1.2)                      # fitted cum. hazard
ax[0].set_xlabel('weeks'); ax[0].set_ylabel('hazard  h(t)  (relapse rate per week)')
ax[0].set_title(f'Fitted hazard by group  (rises since k={kk:.2f}>1; gap = HR)'); ax[0].legend(fontsize=8)
ax[1].set_xlabel('weeks'); ax[1].set_ylabel('cumulative hazard  H(t) = $-\\log S(t)$')
ax[1].set_title('Cumulative hazard: Nelson-Aalen (step) vs Weibull (dashed)'); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.savefig('survival_hazard.png', dpi=120, bbox_inches='tight'); plt.show()

Results & interpretation

• 6-MP works, strongly. The control group's relapse hazard is ~5–6× higher than 6-MP's (HR control-vs-6-MP ≈ 5.7, equivalently 6-MP HR ≈ 0.18), with P(HR>1) = 1.00 — an unambiguous treatment benefit. The Kaplan-Meier curves separate immediately and the fitted Weibull tracks them closely.
• Censoring matters. Ignoring the 12 censored 6-MP patients (treating their last visit as a relapse) would understate the drug's benefit; the censored likelihood credits them with surviving at least that long — visible as the + marks sitting high on the 6-MP curve.
• Shape: $k \approx 1.35$ (CI roughly 0.97–1.74) — a mildly increasing hazard, though $k=1$ (exponential, constant hazard) is borderline plausible. The exponential fit (fix_k=True) gives a similar HR, so the treatment conclusion is robust to the baseline shape.

Takeaways

• Censoring is the whole new idea. Events contribute $h(t)S(t)$, censored cases only $S(t)$ — one indicator $\delta$ switches between them. Everything else (MLE, RW-Metropolis, HR = $e^\beta$) is familiar GLM machinery.
• Survival is a count model underneath. The Weibull log-likelihood is a Poisson kernel with offset $\log t^k$ — the piecewise-exponential trick — so this connects directly to count_reg/Pumps.
• Cross-checks: PyMC (Weibull via a censored log-likelihood) and R (survival::survreg Weibull + coxph) follow; the semiparametric Cox model estimates the same HR without assuming a Weibull shape.