Nonparametric Bayesian Survival¶
A free-form hazard — the Gamma-process prior¶
The survival arc fitted parametric hazards (Weibull, exponential) and the Cox model as a piecewise-exponential Poisson GLM. Here the baseline hazard itself is nonparametric: cut time into $K$ intervals and let the hazard be a free constant $\lambda_k$ on each. Written on person-time the piecewise-exponential likelihood is exactly Poisson — with exposure $e_{ik}$ (time subject $i$ spent in interval $k$) and event indicator $d_{ik}$, $$d_{ik}\sim\text{Poisson}\big(\lambda_k\,e^{x_i'\beta}\,e_{ik}\big).$$
The classic nonparametric prior is the Gamma process (Kalbfleisch 1978): independent increments of the cumulative hazard, i.e. $\lambda_k\sim\text{Gamma}(a,b)$ independently. It is conjugate — $$\lambda_k\mid\cdot\ \sim\ \text{Gamma}\big(a+D_k,\ b+E_k\big),\quad D_k=\textstyle\sum_i d_{ik},\ \ E_k=\textstyle\sum_i e_{ik}e^{x_i'\beta},$$ so the baseline is a direct Gibbs draw. A random-walk prior on $\log\lambda_k$ instead borrows strength across intervals for a smooth hazard. The proportional-hazards coefficient $\beta$ is a Metropolis update.
We validate on a known increasing hazard, fit the Gehan leukemia trial (6-mercaptopurine vs placebo — the arc's Leuk data), compare with the frequentist nonparametric estimators (Kaplan–Meier, Nelson–Aalen) and Cox regression, and cross-check the hazard and hazard ratio in PyMC.
import numpy as np, pandas as pd, matplotlib.pyplot as plt
import npsurv as S
rng = np.random.default_rng(2)
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; ORANGE="#dd6b20"; PURP="#6b46c1"; GREY="#718096"
print("Nonparametric survival: a free hazard on each time interval, tied together by a prior.")
Nonparametric survival: a free hazard on each time interval, tied together by a prior.
1. Does it work? — recovering a known hazard shape¶
Simulated survival data from an increasing hazard $h(t)=0.02\,t$ (a wear-out process), with censoring. The nonparametric model should recover the rising shape from the interval hazards — no Weibull assumed.
haz=lambda t: 0.02*t
time,status=S.simulate_survival(500,haz,20,18,rng)
cuts=np.concatenate([[0],np.quantile(time[status==1],np.linspace(.1,1,10))])
res=S.npsurv_gibbs(time,status,cuts,rng,draws=2000,burn=1000)
mid=0.5*(cuts[:-1]+cuts[1:]); lam=res["lam"]; lo,hi=np.percentile(lam,[2.5,97.5],axis=0)
print(f"recovered hazard vs truth: correlation {np.corrcoef(haz(mid),lam.mean(0))[0,1]:.3f} ({len(time)} subjects, {status.sum()} events)")
fig,ax=plt.subplots(figsize=(8.5,4.2))
ax.step(mid, lam.mean(0), where="mid", color=BLUE, lw=2, label="nonparametric hazard (posterior mean)")
ax.fill_between(mid, lo, hi, step="mid", color=BLUE, alpha=.2, label="95% band")
ax.plot(mid, haz(mid), color="k", lw=1.6, ls="--", label="true hazard 0.02 t")
ax.set_xlabel("time"); ax.set_ylabel("hazard h(t)"); ax.set_title("Recovering an increasing hazard, no parametric form")
ax.legend(frameon=False); plt.tight_layout(); plt.show()
print("The interval hazards climb with the truth: the Gamma-process baseline picks up the wear-out shape a Weibull")
print("would have to assume. The band is wide where events are few (the far tail).")
recovered hazard vs truth: correlation 0.979 (500 subjects, 369 events)
The interval hazards climb with the truth: the Gamma-process baseline picks up the wear-out shape a Weibull would have to assume. The band is wide where events are few (the far tail).
2. The Gehan leukemia trial — baseline hazard and survival¶
42 leukemia patients, remission time in weeks, randomised to 6-mercaptopurine or placebo (12 censored). This is the same Gehan/Leuk dataset the survival arc fits with a parametric Weibull proportional-hazards model (Bayesian Weibull Proportional-Hazards Model) and as a semiparametric Cox model via Poisson regression (Bayesian Cox via Poisson) — here the baseline hazard is left fully nonparametric, so the three notebooks read as one progression on identical data. We fit the nonparametric baseline two ways — the independent Gamma process (each interval free) and the random-walk smoothed hazard (neighbours tied) — and read off survival with a credible band.
d=pd.read_csv("gehan.csv"); t=d.time.to_numpy().astype(float); st=d.status.to_numpy().astype(int)
cutsg=np.concatenate([[0],np.quantile(t[st==1],np.linspace(.15,1,7))])
r_gp=S.npsurv_gibbs(t,st,cutsg,rng,draws=4000,burn=2000,smooth=False)
r_sm=S.npsurv_gibbs(t,st,cutsg,rng,draws=4000,burn=2000,smooth=True)
mid=0.5*(cutsg[:-1]+cutsg[1:])
fig,ax=plt.subplots(1,2,figsize=(12.5,4.4))
ax[0].step(mid,r_gp["lam"].mean(0),where="mid",color=BLUE,lw=2,label="Gamma process (independent)")
ax[0].step(mid,r_sm["lam"].mean(0),where="mid",color=RED,lw=2,ls="--",label="random-walk smoothed")
ax[0].set_xlabel("weeks"); ax[0].set_ylabel("hazard"); ax[0].set_title("Baseline hazard: two priors"); ax[0].legend(frameon=False,fontsize=8)
tg=np.linspace(0,35,120); Sc=S.survival_curves(r_gp,tg)
Slo,Smd,Shi=np.percentile(Sc,[2.5,50,97.5],axis=0)
km_t,km_s=S.kaplan_meier(t,st)
ax[1].fill_between(tg,Slo,Shi,color=BLUE,alpha=.2,label="Bayesian NP 95%"); ax[1].plot(tg,Smd,color=BLUE,lw=2,label="Bayesian NP survival")
ax[1].step(np.r_[0,km_t],np.r_[1,km_s],where="post",color="k",lw=1.6,label="Kaplan–Meier")
ax[1].set_xlabel("weeks"); ax[1].set_ylabel("S(t)"); ax[1].set_title("Overall survival (pooled)"); ax[1].legend(frameon=False,fontsize=8); ax[1].set_ylim(0,1)
plt.tight_layout(); plt.show()
print("The independent Gamma-process hazard is jumpy; the random-walk prior smooths it by borrowing across")
print("intervals. Both give a roughly flat-to-gently-rising hazard -- consistent with the arc's finding that an")
print("exponential fits Leuk well -- and the Bayesian survival curve sits right on the Kaplan–Meier step function.")
The independent Gamma-process hazard is jumpy; the random-walk prior smooths it by borrowing across intervals. Both give a roughly flat-to-gently-rising hazard -- consistent with the arc's finding that an exponential fits Leuk well -- and the Bayesian survival curve sits right on the Kaplan–Meier step function.
3. The frequentist nonparametric estimators — Kaplan–Meier & Nelson–Aalen¶
The classical nonparametric survival tools need no prior: Kaplan–Meier estimates $S(t)$ as a product over event times, Nelson–Aalen the cumulative hazard as a sum of $d/n$ increments. They are the frequentist limit of the Gamma-process posterior (a Gamma process with vanishing prior weight is the Nelson–Aalen estimator), so they should coincide with our posterior mean.
na_t,na_H=S.nelson_aalen(t,st)
Hpost=S.cum_hazard(r_gp,tg); Hlo,Hmd,Hhi=np.percentile(Hpost,[2.5,50,97.5],axis=0)
fig,ax=plt.subplots(1,2,figsize=(12.5,4.3))
ax[0].fill_between(tg,Hlo,Hhi,color=GREEN,alpha=.2,label="Bayesian NP 95%"); ax[0].plot(tg,Hmd,color=GREEN,lw=2,label="Bayesian cumulative hazard")
ax[0].step(np.r_[0,na_t],np.r_[0,na_H],where="post",color="k",lw=1.6,label="Nelson–Aalen")
ax[0].set_xlabel("weeks"); ax[0].set_ylabel("H(t)"); ax[0].set_title("Cumulative hazard: Bayesian vs Nelson–Aalen"); ax[0].legend(frameon=False,fontsize=8)
Sc=S.survival_curves(r_gp,tg); Smd=np.percentile(Sc,50,axis=0)
ax[1].plot(tg,Smd,color=BLUE,lw=2,label="Bayesian NP"); ax[1].step(np.r_[0,km_t],np.r_[1,km_s],where="post",color="k",lw=1.6,ls="--",label="Kaplan–Meier")
ax[1].set_xlabel("weeks"); ax[1].set_ylabel("S(t)"); ax[1].set_title("Survival: Bayesian vs Kaplan–Meier"); ax[1].legend(frameon=False,fontsize=8); ax[1].set_ylim(0,1)
plt.tight_layout(); plt.show()
print("The Bayesian nonparametric posterior mean traces Nelson–Aalen and Kaplan–Meier almost exactly -- the prior")
print("adds a credible band and a little smoothing, and would matter more with sparser data, but here the estimators agree.")
The Bayesian nonparametric posterior mean traces Nelson–Aalen and Kaplan–Meier almost exactly -- the prior adds a credible band and a little smoothing, and would matter more with sparser data, but here the estimators agree.
4. Proportional hazards — the treatment effect¶
Now the covariate: does 6-mercaptopurine lower the hazard? With the nonparametric baseline and a proportional-hazards term $e^{x'\beta}$, $\beta$ is sampled by Metropolis and $e^{\beta}$ is the hazard ratio. We compare with the semiparametric Cox model (frequentist), which conditions the baseline away entirely.
X=d[["treat"]].to_numpy().astype(float) # 1 = 6-MP, 0 = control
r_ph=S.npsurv_gibbs(t,st,cutsg,rng,X=X,draws=6000,burn=3000,step=0.25)
hr=np.exp(r_ph["beta"][:,0])
print(f"Bayesian NP proportional hazards: HR(6-MP) = {np.median(hr):.2f} 95% CrI [{np.percentile(hr,2.5):.2f}, {np.percentile(hr,97.5):.2f}]")
# Tie handling matters here: 18 of the 42 event times are tied, and the two standard
# approximations to the partial likelihood give visibly different answers. statsmodels
# defaults to Breslow; R's survival::coxph defaults to Efron. Report both, so the small
# gap against the R notebook is identified as a CONVENTION and not a disagreement.
try:
from statsmodels.duration.hazard_regression import PHReg
n_tied = len(t) - len(np.unique(t))
for ties in ("breslow", "efron"):
cox = PHReg(t, X, status=st, ties=ties).fit()
chr_ = np.exp(cox.params[0]); cse = cox.bse[0]
print("frequentist Cox model (%-7s): HR(6-MP) = %.2f 95%% CI [%.2f, %.2f]%s"
% (ties, chr_, np.exp(cox.params[0]-1.96*cse), np.exp(cox.params[0]+1.96*cse),
" <- the default in R's coxph" if ties == "efron" else ""))
print("(%d of %d event times are tied, which is why the two conventions separate at all.)" % (n_tied, len(t)))
except Exception as e:
print("Cox (statsmodels) unavailable:", str(e)[:60])
tg=np.linspace(0,35,120)
Sc_ctrl=S.survival_curves(r_ph,tg,x=[0.0]); Sc_trt=S.survival_curves(r_ph,tg,x=[1.0])
fig,ax=plt.subplots(figsize=(9,4.4))
for Sc,c,lab,grp in [(Sc_ctrl,RED,"control (placebo)",0),(Sc_trt,BLUE,"6-MP",1)]:
lo,md_,hi=np.percentile(Sc,[2.5,50,97.5],axis=0); ax.fill_between(tg,lo,hi,color=c,alpha=.15); ax.plot(tg,md_,color=c,lw=2.2,label=lab)
kt,ks=S.kaplan_meier(t[X[:,0]==grp],st[X[:,0]==grp]); ax.step(np.r_[0,kt],np.r_[1,ks],where="post",color=c,lw=1,ls=":")
ax.set_xlabel("weeks in remission"); ax.set_ylabel("S(t)"); ax.set_title(f"Survival by arm (HR={np.median(hr):.2f}); dotted = Kaplan–Meier"); ax.legend(frameon=False); ax.set_ylim(0,1)
plt.tight_layout(); plt.show()
print(f"6-MP cuts the hazard of relapse to about a fifth of placebo (HR~{np.median(hr):.2f}) -- the celebrated result of this")
print("1963 trial. The Bayesian PH estimate matches Cox, but also delivers the full baseline hazard and survival curves.")
Bayesian NP proportional hazards: HR(6-MP) = 0.19 95% CrI [0.08, 0.39]
frequentist Cox model (breslow): HR(6-MP) = 0.22 95% CI [0.10, 0.49] frequentist Cox model (efron ): HR(6-MP) = 0.21 95% CI [0.09, 0.47] <- the default in R's coxph (18 of 42 event times are tied, which is why the two conventions separate at all.)
6-MP cuts the hazard of relapse to about a fifth of placebo (HR~0.19) -- the celebrated result of this 1963 trial. The Bayesian PH estimate matches Cox, but also delivers the full baseline hazard and survival curves.
5. Cross-check in PyMC — the piecewise-exponential Poisson GLM¶
The person-time likelihood is Poisson, so the whole model is a Poisson GLM on the expanded data: one row per (subject, interval) at risk, an offset $\log e_{ik}$, an intercept per interval (the log baseline hazard) and the treatment coefficient. This is the same Cox-via-Poisson bridge used in the survival arc; here the interval intercepts are the nonparametric baseline. PyMC recovers the hazard ratio and the baseline.
import pymc as pm
E,D=S.expand_pwe(t,st,cutsg); N,K=E.shape
rows_i,rows_k=np.nonzero(E) # (subject, interval) at risk
expo=E[rows_i,rows_k]; dev=D[rows_i,rows_k]; kk=rows_k; trt=X[rows_i,0]
with pm.Model() as mod:
logh=pm.Normal("logh",np.log(0.05),1.5,shape=K) # log baseline hazard per interval
beta=pm.Normal("beta",0,5)
mu=pm.math.exp(logh[kk]+beta*trt)*expo
pm.Poisson("d",mu=mu,observed=dev)
idata=pm.sample(1000,tune=1500,chains=4,target_accept=0.9,random_seed=3,progressbar=False)
hr_pm=np.exp(idata.posterior["beta"].values).ravel()
lam_pm=np.exp(idata.posterior["logh"].mean(("chain","draw")).values)
print(f"PyMC HR(6-MP) = {np.median(hr_pm):.2f} [{np.percentile(hr_pm,2.5):.2f}, {np.percentile(hr_pm,97.5):.2f}] from-scratch {np.median(hr):.2f}")
print(f"baseline hazard agreement: max |PyMC - Gibbs| = {np.abs(lam_pm-r_ph['lam'].mean(0)).max():.3f}")
fig,ax=plt.subplots(figsize=(8.5,4))
ax.step(mid,r_ph["lam"].mean(0),where="mid",color=BLUE,lw=2.2,label="from-scratch Gibbs")
ax.step(mid,lam_pm,where="mid",color=RED,lw=1.6,ls="--",label="PyMC Poisson GLM")
ax.set_xlabel("weeks"); ax.set_ylabel("baseline hazard"); ax.set_title("Nonparametric baseline hazard: from-scratch vs PyMC"); ax.legend(frameon=False)
plt.tight_layout(); plt.show()
lo_pm, hi_pm = np.percentile(hr_pm, [2.5, 97.5]); lo_fs, hi_fs = np.percentile(hr, [2.5, 97.5])
ov = (min(hi_pm, hi_fs) - max(lo_pm, lo_fs)) / (max(hi_pm, hi_fs) - min(lo_pm, lo_fs))
base = float(r_ph["lam"].mean(0).mean()); gapb = float(np.abs(lam_pm - r_ph["lam"].mean(0)).max())
print("\nThe two are one model, but read the agreement off the intervals rather than the point estimates:")
print(" HR PyMC %.2f [%.2f, %.2f] vs Gibbs %.2f [%.2f, %.2f]"
% (np.median(hr_pm), lo_pm, hi_pm, np.median(hr), lo_fs, hi_fs))
print(" the medians differ by %.0f%%, which sounds like a lot until you notice the two credible"
% (100*abs(np.median(hr_pm)-np.median(hr))/np.median(hr)))
print(" intervals overlap over %.0f%% of their combined range -- with 30 events, HR is simply not" % (100*ov))
print(" pinned down to two significant figures, and no amount of agreeing samplers will change that.")
print(" baseline hazard: max gap %.3f against a mean baseline of %.3f, about %.0f%%."
% (gapb, base, 100*gapb/base))
print("\nWhat IS exact is the correspondence itself: a Gamma-process baseline with a proportional-hazards")
print("term and a piecewise-exponential Poisson GLM are the same likelihood written two ways. The residual")
print("differences are Monte Carlo error and prior choice, not a difference of model.")
g++ not available, if using conda: `conda install gxx`
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [logh, beta]
Sampling 4 chains for 1_500 tune and 1_000 draw iterations (6_000 + 4_000 draws total) took 5 seconds.
PyMC HR(6-MP) = 0.24 [0.11, 0.53] from-scratch 0.19 baseline hazard agreement: max |PyMC - Gibbs| = 0.083
The two are one model, but read the agreement off the intervals rather than the point estimates: HR PyMC 0.24 [0.11, 0.53] vs Gibbs 0.19 [0.08, 0.39] the medians differ by 25%, which sounds like a lot until you notice the two credible intervals overlap over 64% of their combined range -- with 30 events, HR is simply not pinned down to two significant figures, and no amount of agreeing samplers will change that. baseline hazard: max gap 0.083 against a mean baseline of 0.152, about 55%. What IS exact is the correspondence itself: a Gamma-process baseline with a proportional-hazards term and a piecewise-exponential Poisson GLM are the same likelihood written two ways. The residual differences are Monte Carlo error and prior choice, not a difference of model.
6. Summary¶
A nonparametric hazard: cut time into intervals, let each carry a free rate, and tie them with a prior. The Gamma-process prior (independent increments of the cumulative hazard) is conjugate, so the baseline is a direct Gibbs draw; a random-walk prior smooths it. On the Gehan leukemia trial the free hazard recovered a roughly exponential baseline (validating the arc's parametric fit), the Bayesian survival curve sat exactly on Kaplan–Meier and the cumulative hazard on Nelson–Aalen, and the proportional-hazards term returned the famous HR ≈ 0.2 for 6-mercaptopurine — matching the frequentist Cox model but with the full baseline hazard in hand.
The engine is the piecewise-exponential Poisson GLM, confirmed in PyMC: it is the same Cox-via-Poisson bridge the survival arc used, now with a nonparametric baseline rather than a parametric one. So the arc's progression on this one dataset is complete — parametric Weibull (Bayesian Weibull Proportional-Hazards Model), semiparametric Cox (Bayesian Cox via Poisson), and now a fully nonparametric hazard — and the frequentist Kaplan–Meier / Nelson–Aalen / Cox estimators are the priorless limit of the same picture. Next the arc closes with Polya trees, a nonparametric prior directly on the distribution.