Mastectomy — R cross-check (survSplit + glm poisson, vs coxph)¶

Companion to mastectomy_python.ipynb / mastectomy_pymc.ipynb. survival::survSplit does the person-time expansion; glm(poisson) with interval factors + a log(exposure) offset is the Cox fit; coxph is the reference.

In [1]:
suppressMessages(library(survival))
d <- read.csv('mastectomy.csv'); d$met <- as.integer(d$metastasized=='yes'); d$id <- seq_len(nrow(d))
cut <- sort(unique(d$time[d$event==1])); cut <- cut[-length(cut)]
sp <- survSplit(Surv(time, event) ~ ., data=d, cut=cut, episode='interval', start='tstart')
sp$expo <- sp$time - sp$tstart; sp <- sp[sp$expo > 0, ]
gl <- glm(event ~ factor(interval) + met + offset(log(expo)), family=poisson, data=sp)
bt <- coef(gl)[['met']]; se <- sqrt(diag(vcov(gl)))[['met']]
cat(sprintf('Poisson glm (Cox bridge): metastasis coef = %.3f (se %.3f)  HR = %.2f\n', bt, se, exp(bt)))
cx <- coxph(Surv(time, event) ~ met, data=d)
cat(sprintf('coxph (reference)       : metastasis coef = %.3f          HR = %.2f  p = %.3f\n',
            coef(cx)[['met']], exp(coef(cx)[['met']]), summary(cx)$coefficients['met','Pr(>|z|)']))
cat(sprintf('rows after split: %d  (from %d patients, %d intervals)\n', nrow(sp), nrow(d), length(cut)+1))

Poisson glm (Cox bridge): metastasis coef = 0.835 (se 0.501)  HR = 2.31

coxph (reference)       : metastasis coef = 0.852          HR = 2.34  p = 0.090

rows after split: 750  (from 44 patients, 25 intervals)

Results¶

  • The Poisson glm on the survSplit person-time data returns metastasis HR ≈ 2.3 — matching the from-scratch and PyMC Poisson fits.
  • coxph gives the same effect (HR ≈ 2.3) as the partial-likelihood reference.
  • Same conclusion three ways, on a second dataset: the Cox model is a Poisson GLM on expanded follow-up time. Metastasis roughly doubles the death hazard.