Bayesian Cox via Poisson

Model

The Cox proportional-hazards model — the semiparametric workhorse of survival analysis, semiparametric meaning it puts a parametric form on the covariate effects while leaving the baseline hazard completely unspecified — can be fit as an ordinary Poisson GLM on expanded person-time data: rather than one row per subject, one row per subject per stretch of time they were under observation, so the exposure each subject contributes is recorded explicitly rather than implied. 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) — partial because Cox's insight was that at each event one need only ask which member of the risk set (everyone still under observation at that instant) the event happened to, a question whose answer depends on the covariates but not at all on the baseline hazard, so the baseline drops out of the likelihood rather than having to be estimated: $\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)