Bayesian Cox via Poisson

Model

The Cox proportional-hazards model — the semiparametric workhorse of survival analysis — can be fit as an ordinary Poisson GLM on expanded person-time data. This is the literal form used by BUGS Vol I Leuk and PyMC's "Bayesian survival analysis" case study, and it closes the loop on the recurring theme of this site: survival is a count model. The earlier Weibull example noted its likelihood was a Poisson kernel; here that observation becomes the estimation method itself, with no parametric baseline assumed.

y_{ig} \sim \text{Poisson}(\mu_{ig}), \qquad \log\mu_{ig} = \underbrace{\log(\text{exposure}_{ig})}_{\text{offset}} + \lambda_g + x_i'\beta

$y_{ig}$ — 1 if subject $i$ failed in interval $g$ , else 0
$\lambda_g$ — log baseline hazard in interval $g$ (nonparametric step baseline)
$e^{\beta}$ — hazard ratio; $\text{exposure}_{ig}$ — time at risk in interval $g$

The piecewise-exponential = Poisson identity

Split the time axis at the distinct event times into intervals $g$ with a piecewise-constant baseline hazard $\lambda_g$ , and expand each subject into one row per interval they are at risk in. Then $y_{ig}$ (1 if subject $i$ failed in interval $g$ ) is Poisson with a log-exposure offset — exactly the count regression model with interval intercepts. As the cuts approach every event time, $\hat\beta$ converges to the Cox partial-likelihood estimate (Holford 1980; Laird & Olivier 1981): $\lambda_g$ is a nonparametric step baseline (no Weibull shape imposed), and $e^{\beta}$ is the hazard ratio — the same quantity as coxph, obtained from a count model.

Why the bridge matters

Why it matters: anything available for Poisson GLMs immediately enriches the survival model — splines on $\lambda_g$ , extra covariates, or random effects. Adding a patient random intercept to this Poisson turns it into the frailty model as a hierarchical Poisson — the same modular pieces, recombined.

Notebooks

Two datasets, identical machinery (the same coxpois.py, only the covariate changes). Leuk (Gehan leukemia, the Weibull project's data): the Poisson-on-person-time fit recovers control-vs-6-MP HR $\approx 5$ — matching coxph (4.8) and the Weibull fit (5.6), but with the baseline estimated nonparametrically (it tracks the Weibull baseline closely, which is why the parametric choice was reasonable — without having to assume it). Mastectomy (PyMC's canonical example; 44 breast-cancer patients): metastasis raises the death hazard ~2× (MLE/coxph HR $\approx 2.3$ , the Bayesian posterior a milder HR $\approx 1.9\text{–}2.0$ with $P(\text{HR}>1)\approx 0.92$ ) — its credible interval still includes 1, an uncertain effect with only 26 deaths. Each is cross-checked three ways — from-scratch Poisson (MLE + RW-Metropolis), PyMC pm.Poisson on the same expansion, and R survSplit + glm(poisson) — with coxph as the partial-likelihood reference the piecewise-exponential Poisson approximates.

Python Leuk — Cox via Poisson from scratch: person-time expansion, MLE + RW-Metropolis (Python) Open → Python Leuk — PyMC cross-check: pm.Poisson on the person-time data (Python) Open → R Leuk — R cross-check: survSplit + glm(poisson) vs coxph reference (R) Open → Python Mastectomy — Cox via Poisson: metastasis effect on breast-cancer survival (Python) Open → Python Mastectomy — PyMC cross-check: pm.Poisson (Python) Open → R Mastectomy — R cross-check: survSplit + glm(poisson) vs coxph (R) Open →

Downloads

coxpois.py Person-time expansion + offset-Poisson MLE & RW-Metropolis (NumPy / scipy) leuk.csv Gehan leukemia trial — 42 patients, 6-MP vs placebo, weeks to relapse mastectomy.csv Mastectomy breast-cancer survival — 44 women, metastasis covariate

References

Cox, D. R. (1972). Regression models and life-tables. Journal of the Royal Statistical Society: Series B 34(2), 187–220. — the Cox proportional-hazards model and partial likelihood this method reproduces
Holford, T. R. (1980). The analysis of rates and of survivorship using log-linear models. Biometrics 36(2), 299–305. — the piecewise-exponential = Poisson-GLM equivalence
Laird, N. & Olivier, D. (1981). Covariance analysis of censored survival data using log-linear analysis techniques. Journal of the American Statistical Association 76(374), 231–240. — the person-time Poisson formulation of Cox regression
Aitkin, M. & Clayton, D. (1980). The fitting of exponential, Weibull and extreme value distributions to complex censored survival data using GLIM. Journal of the Royal Statistical Society: Series C 29(2), 156–163. — fitting survival models as generalized linear models
Freireich, E. J. et al. (1963). The effect of 6-mercaptopurine on the duration of steroid-induced remissions in acute leukemia. Blood 21(6), 699–716. — the leukemia remission trial (the Leuk data)

Sampler Module — Source Code

"""
Cox proportional hazards via the PIECEWISE-EXPONENTIAL = POISSON identity -- the count<->survival bridge.
BUGS Vol I "Leuk" is written in exactly this form; the PyMC "Bayesian survival analysis" (mastectomy)
case study uses it too.

Idea
----
Split the time axis at the distinct event times into intervals g with a piecewise-CONSTANT baseline
hazard lambda_g. Expand each subject into one "person-time" row per interval they are at risk in. Then
  y_ig ~ Poisson(mu_ig),   log mu_ig = log(exposure_ig) + lambda_g + x_i' beta
where y_ig = 1 if subject i had the event in interval g (else 0), exposure_ig = time at risk in interval
g (the OFFSET), and lambda_g = log baseline hazard. This is just a Poisson GLM with an offset and
interval intercepts -> fit it with count_reg_mcmc.poisson_*. As the cuts -> every event time, beta is the
Cox partial-likelihood estimate (Holford 1980; Laird & Olivier 1981). No parametric baseline assumed.

survsplit() does the person-time expansion; build_design() assembles [interval dummies | covariates].
The Poisson fitters below are SELF-CONTAINED (same model as count_reg_mcmc's Poisson, re-implemented here
with a log-exposure OFFSET so this project doesn't modify the older count-regression module).
"""
import numpy as np
from scipy.special import gammaln


def survsplit(time, status, cut):
    """Expand to person-time rows. `cut` = interior cutpoints (sorted). Intervals: (0,c1],...,(c_{K-1},inf).
    Returns (subj, interval, exposure, event) arrays, one entry per (subject, at-risk interval)."""
    time = np.asarray(time, float); status = np.asarray(status, float)
    edges = np.concatenate([[0.0], np.asarray(cut, float)]); G = len(edges)   # G intervals (last open)
    subj = []; itv = []; expo = []; ev = []
    for i, (ti, di) in enumerate(zip(time, status)):
        for g in range(G):
            lo = edges[g]; hi = edges[g + 1] if g + 1 < G else np.inf
            if ti <= lo:
                break
            exits = ti <= hi
            subj.append(i); itv.append(g); expo.append(min(ti, hi) - lo); ev.append(di if exits else 0.0)
            if exits:
                break
    return np.array(subj, int), np.array(itv, int), np.array(expo, float), np.array(ev, float)


def cutpoints(time, status):
    """Distinct event times, minus the last (interior breakpoints) so the final interval holds the last
    event plus any later censoring -> every interval contains >=1 event."""
    ev_times = np.unique(np.asarray(time)[np.asarray(status) == 1])
    return ev_times[:-1]


def build_design(itv, G, covars):
    """[interval dummies (G one-hot cols) | covariate columns]. No global intercept."""
    D = np.zeros((len(itv), G))
    D[np.arange(len(itv)), itv] = 1.0
    return np.column_stack([D] + [np.asarray(c, float) for c in covars])


# ── self-contained Poisson regression with a log-exposure offset (the piecewise-exp likelihood) ──
def pois_mle(X, y, offset, tol=1e-10, maxit=200):
    """Newton/IRLS Poisson MLE with offset. log mu = X beta + offset. Returns (beta, se)."""
    k = X.shape[1]; beta = np.zeros(k)
    for _ in range(maxit):
        mu = np.exp(X @ beta + offset)
        I = (X * mu[:, None]).T @ X
        step = np.linalg.solve(I, X.T @ (y - mu)); beta += step
        if np.max(np.abs(step)) < tol:
            break
    cov = np.linalg.inv((X * np.exp(X @ beta + offset)[:, None]).T @ X)
    return beta, np.sqrt(np.diag(cov))


def pois_rwm(y, X, offset, R=10000, burn=2000, prior_sd=10.0, seed=0, s=None):
    """RW-Metropolis for the offset-Poisson; proposal scaled by the Fisher info at the MLE."""
    rng = np.random.default_rng(seed); n, k = X.shape
    b0, _ = pois_mle(X, y, offset)
    I = (X * np.exp(X @ b0 + offset)[:, None]).T @ X
    L = np.linalg.cholesky(np.linalg.inv(I + np.eye(k) / prior_sd ** 2))
    s = (2.38 / np.sqrt(k)) if s is None else s

    def lp(b):
        eta = X @ b + offset
        return np.sum(y * eta - np.exp(eta)) - 0.5 * np.sum(b ** 2) / prior_sd ** 2

    beta = b0.copy(); cur = lp(beta); keep = R - burn; B = np.zeros((keep, k)); acc = 0
    for r in range(R):
        prop = beta + s * (L @ rng.standard_normal(k))
        pr = lp(prop)
        if np.log(rng.random()) < pr - cur:
            beta, cur = prop, pr; acc += 1
        if r >= burn:
            B[r - burn] = beta
    return dict(beta=B, accept=acc / R)