Mastectomy — Cox via piecewise-exponential = Poisson¶

The canonical PyMC survival example; same bridge as coxpois_python.ipynb, different data¶

Methodologically identical to coxpois_python.ipynb (Cox proportional hazards fit as a Poisson GLM on person-time data — see that notebook for the full piecewise-exponential = Poisson derivation). It reuses the same module coxpois.py. The only change is the dataset: this is the mastectomy breast-cancer study that PyMC's "Bayesian Survival Analysis" case study made the standard demonstration of this technique.

The data — mastectomy (breast-cancer survival)¶

44 women followed after a mastectomy for breast cancer. The clinical question: does tumour metastasis (spread beyond the breast, detected by immunohistochemical staining) shorten survival? Survival is right-censored — 18 women were still alive at last follow-up.

field meaning
time months from surgery to death or last follow-up; range 5–225, median 74
event 1 = died (26), 0 = censored / alive at last contact (18)
metastasized yes (32 women) / no (12) — whether the cancer had metastasized

The signal is visible in the raw data: metastasized patients died in 66% of cases (median 64 months) vs 42% (median 124 months) for non-metastasized — roughly half the survival time. The model's job is to quantify that hazard difference after properly crediting the censored survivors.

In [1]:
import numpy as np, pandas as pd, matplotlib.pyplot as plt
from coxpois import survsplit, cutpoints, build_design, pois_mle, pois_rwm
from survival_mcmc import weibull_mle, km_estimator

d = pd.read_csv('mastectomy.csv')
met = (d['metastasized'] == 'yes').astype(float).to_numpy()
t = d['time'].to_numpy(float); ev = d['event'].to_numpy(float)

cut = cutpoints(t, ev); G = len(cut) + 1
subj, itv, expo, evt = survsplit(t, ev, cut)
X = build_design(itv, G, [met[subj]]); off = np.log(expo)
print('%d patients -> %d person-time rows over G=%d intervals; deaths %d' % (len(d), len(evt), G, int(evt.sum())))

b_mle, se = pois_mle(X, evt, off)
o = pois_rwm(evt, X, off, R=12000, burn=3000, seed=1); bt = o['beta'][:, -1]
print('\nMETASTASIS effect, via the Poisson bridge:')
print('  Poisson MLE   coef %.3f (se %.3f)   HR %.2f' % (b_mle[-1], se[-1], np.exp(b_mle[-1])))
print('  Poisson Bayes coef %.3f [%.2f, %.2f]  HR %.2f  P(HR>1) %.3f  (acc %.2f)'
      % (bt.mean(), np.percentile(bt,2.5), np.percentile(bt,97.5), np.exp(bt.mean()), (bt>0).mean(), o['accept']))

44 patients -> 750 person-time rows over G=25 intervals; deaths 26

METASTASIS effect, via the Poisson bridge:
  Poisson MLE   coef 0.835 (se 0.501)   HR 2.31
  Poisson Bayes coef 0.679 [-0.24, 1.82]  HR 1.97  P(HR>1) 0.911  (acc 0.25)
In [2]:
# same two graphs as coxpois_python: baseline cumulative hazard, and KM vs Cox-Poisson survival by group
edges = np.concatenate([[0.0], cut]); lam = b_mle[:G]
def H0_cox(tg):
    H = np.zeros_like(tg, float)
    for j, tt in enumerate(tg):
        h = 0.0
        for g in range(G):
            lo = edges[g]; hi = edges[g+1] if g+1 < G else np.inf
            if tt <= lo: break
            h += np.exp(lam[g]) * (min(tt, hi) - lo)
        H[j] = h
    return H
Xw = np.column_stack([np.ones(len(d)), met]); mw, _, _ = weibull_mle(Xw, t, ev)
kw = np.exp(mw[2]); lam0_w = np.exp(mw[0]); bt_mean = bt.mean()
grid = np.linspace(1, 225, 250); H0c = H0_cox(grid); H0w = lam0_w * grid ** kw

fig, ax = plt.subplots(1, 2, figsize=(12, 4.6))
ax[0].step(grid, H0c, where='post', color='steelblue', lw=1.8, label='piecewise (Poisson) — nonparametric')
ax[0].plot(grid, H0w, '--', color='firebrick', lw=1.4, label='Weibull — parametric')
ax[0].set_xlabel('months'); ax[0].set_ylabel('baseline cumulative hazard $H_0(t)$')
ax[0].set_title('Baseline hazard: semiparametric vs Weibull'); ax[0].legend(fontsize=8)
for g, lab, col in [(0,'not metastasized','steelblue'), (1,'metastasized','firebrick')]:
    tk, Sk = km_estimator(t[met==g], ev[met==g])
    ax[1].step(tk, Sk, where='post', color=col, lw=1.6, label='KM '+lab)
    ax[1].plot(grid, np.exp(-H0c * np.exp(bt_mean*g)), '--', color=col, lw=1.2)
ax[1].set_xlabel('months'); ax[1].set_ylabel('S(t)'); ax[1].set_ylim(0,1.02)
ax[1].set_title('KM (step) vs Cox-Poisson fit (dashed)'); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.savefig('mastectomy_leuk.png'.replace('_leuk',''), dpi=120, bbox_inches='tight'); plt.show()
No description has been provided for this image

Results¶

  • Metastasis raises the death hazard ~2× (MLE HR ≈ 2.3, coef ≈ 0.84; the Bayesian posterior a milder HR ≈ 2.0), with P(HR>1) ≈ 0.91 — evidence of higher risk, though the 95% credible interval still includes 1 (only 44 patients, 26 deaths). This is the well-known result from the PyMC case study.
  • The nonparametric (piecewise) baseline tracks the Weibull baseline closely, and the Cox-Poisson survival curves sit on the Kaplan-Meier steps for both groups — the metastasized arm (red) falls away faster.
  • Same model, different data: this used coxpois.py unchanged — only the dataset and covariate (metastasized) changed versus the leukemia analysis. The hazard ratio is, once again, just a Poisson-regression coefficient on expanded follow-up time.
  • Cross-checks: PyMC (pm.Poisson) and R (survSplit + glm, with coxph reference) follow.