Cox via piecewise-exponential = Poisson — the count↔survival bridge¶

Fitting a semiparametric survival model with the count machinery¶

Module: coxpois.py. The Cox proportional-hazards model can be fit as a Poisson GLM on expanded person-time data — the form BUGS Vol I "Leuk" actually uses, and the PyMC "Bayesian survival analysis" (mastectomy) approach. We apply it to the same Gehan leukemia data as survival_python.ipynb, but now with no parametric baseline — the baseline hazard is estimated nonparametrically.

The identity¶

Split time at the distinct event times into intervals $g$ with a piecewise-constant baseline hazard $\lambda_g$. Expand each subject into one row per interval they're at risk in. Then $$y_{ig}\sim\text{Poisson}(\mu_{ig}),\qquad \log\mu_{ig}=\underbrace{\log(\text{exposure}_{ig})}_{\text{offset}}+\lambda_g+x_i'\beta,$$ where $y_{ig}=1$ if subject $i$ failed in interval $g$ (else 0), and $\text{exposure}_{ig}$ = time at risk in $g$. This is exactly a Poisson regression with an offset and interval intercepts — the count_reg model. As the cuts approach every event time, $\hat\beta\to$ the Cox partial-likelihood estimate (Holford 1980; Laird & Olivier 1981).

  • $\lambda_g$ = log baseline hazard in interval $g$ → a nonparametric (step) baseline, no Weibull assumed.
  • $e^{\beta}$ = hazard ratio, the same quantity as coxph / the Weibull fit — but obtained from a count model.

The Poisson fitters (pois_mle, pois_rwm) live in coxpois.py, self-contained with a log-exposure offset (the same Poisson model as count_reg_mcmc, kept separate so this project doesn't alter the older one).

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('leuk.csv')
control = (d['treat'] == 'control').astype(float).to_numpy()
t = d['time'].to_numpy(float); status = d['cens'].to_numpy(float)

cut = cutpoints(t, status); G = len(cut) + 1
subj, itv, expo, ev = survsplit(t, status, cut)          # person-time expansion
X = build_design(itv, G, [control[subj]]); off = np.log(expo)
print('%d subjects -> %d person-time rows over G=%d intervals; events %d' % (len(d), len(ev), G, int(ev.sum())))

b_mle, se = pois_mle(X, ev, off)
o = pois_rwm(ev, X, off, R=12000, burn=3000, seed=1); bt = o['beta'][:, -1]
print('\nTREATMENT (control vs 6-MP), 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']))
print('  --- same effect from coxph (1.57, HR 4.8) and the Weibull fit (1.73, HR 5.6) ---')

42 subjects -> 421 person-time rows over G=17 intervals; events 30

TREATMENT (control vs 6-MP), via the Poisson bridge:
  Poisson MLE   coef 1.648 (se 0.433)   HR 5.20

  Poisson Bayes coef 1.631 [0.77, 2.49]  HR 5.11  P(HR>1) 1.000  (acc 0.26)
  --- same effect from coxph (1.57, HR 4.8) and the Weibull fit (1.73, HR 5.6) ---
In [2]:
# baseline cumulative hazard: nonparametric (piecewise from the Poisson) vs Weibull
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
# Weibull baseline (6-MP reference) for comparison
Xw = np.column_stack([np.ones(len(d)), control]); mw, _, _ = weibull_mle(Xw, t, status)
kw = np.exp(mw[2]); lam0_w = np.exp(mw[0]); bt_mean = bt.mean()
grid = np.linspace(1, 35, 200); 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('weeks'); ax[0].set_ylabel('baseline cumulative hazard $H_0(t)$'); ax[0].set_title('Baseline hazard: semiparametric vs Weibull'); ax[0].legend(fontsize=8)
# fitted survival by group (from the Cox-Poisson baseline) vs KM
for g, lab, col in [(0,'6-MP','steelblue'), (1,'control','firebrick')]:
    tk, Sk = km_estimator(t[control==g], status[control==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('weeks'); 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('coxpois_leuk.png', dpi=120, bbox_inches='tight'); plt.show()
No description has been provided for this image

Results¶

  • The bridge holds: fitting a Poisson regression on person-time data recovers the treatment hazard ratio — control vs 6-MP HR ≈ 5 (coef ≈ 1.6, P(HR>1) ≈ 1) — the same effect as coxph (HR 4.8) and the Weibull fit (HR 5.6). Survival was estimated with the count sampler.
  • Nonparametric baseline: the left panel shows the piecewise (step) baseline cumulative hazard the Poisson estimates — no Weibull shape imposed. It tracks the Weibull baseline closely here, which is why the Weibull was a reasonable parametric choice for this data; but the Cox-Poisson didn't have to assume it.
  • The fitted survival curves (right) sit on the Kaplan-Meier steps for both arms.

Takeaways¶

  • Survival = count model. A semiparametric Cox fit is a Poisson GLM with a log-exposure offset and interval intercepts — the count_reg machinery applied to expanded follow-up time. The hazard ratio is just a Poisson regression coefficient.
  • Why it matters: this is the formulation behind BUGS Leuk; it also means anything you can do for Poisson GLMs (random effects, splines on $\lambda_g$, extra covariates) immediately gives richer survival models — e.g. adding a patient random effect turns this into the frailty model as a hierarchical Poisson.
  • Cross-checks: PyMC (pm.Poisson on the same expansion) and R (survival::survSplit + glm(poisson), with coxph as the reference) follow.