Cox via Poisson — R cross-check (survSplit + glm poisson, vs coxph)

Companion to coxpois_python.ipynb / coxpois_pymc.ipynb. R's survival::survSplit does the person-time expansion; a plain glm(..., family=poisson) with interval factors and a log(exposure) offset then is the Cox fit. coxph is the reference.

suppressMessages(library(survival))
d <- read.csv('leuk.csv'); d$control <- as.integer(d$treat=='control'); d$id <- seq_len(nrow(d))
cut <- sort(unique(d$time[d$cens==1])); cut <- cut[-length(cut)]            # interior event-time cutpoints
sp <- survSplit(Surv(time, cens) ~ ., data=d, cut=cut, episode='interval', start='tstart')
sp$expo <- sp$time - sp$tstart; sp <- sp[sp$expo > 0, ]
gl <- glm(cens ~ factor(interval) + control + offset(log(expo)), family=poisson, data=sp)
bt <- coef(gl)[['control']]; se <- sqrt(diag(vcov(gl)))[['control']]
cat(sprintf('Poisson glm (Cox bridge): treat coef = %.3f (se %.3f)  HR = %.2f\n', bt, se, exp(bt)))
cx <- coxph(Surv(time, cens) ~ control, data=d)
cat(sprintf('coxph (reference)       : treat coef = %.3f          HR = %.2f\n', coef(cx)[['control']], exp(coef(cx)[['control']])))
cat(sprintf('rows after split: %d  (from %d subjects, %d intervals)\n', nrow(sp), nrow(d), length(cut)+1))

Poisson glm (Cox bridge): treat coef = 1.648 (se 0.433)  HR = 5.20

coxph (reference)       : treat coef = 1.572          HR = 4.82

rows after split: 421  (from 42 subjects, 17 intervals)

Results

• The Poisson glm on the survSplit person-time data returns treat HR ≈ 5 — matching the from-scratch and PyMC Poisson fits exactly (same model).
• coxph gives HR ≈ 4.8 — the partial-likelihood estimate the piecewise-exponential Poisson approximates (the small gap is tie-handling: Breslow vs the interval baseline). As the cuts → every event time, the two coincide.
• Three Poisson engines + coxph agree: the Cox model really is a Poisson GLM on expanded follow-up time. This is the literal form of BUGS "Leuk."