Mice — R cross-check (survival::survreg)

Companion to mice_python.ipynb / mice_pymc.ipynb. Weibull AFT via survreg (converted to the PH shape/HR) plus Kaplan-Meier by group.

suppressMessages(library(survival))
d <- read.csv('mice.csv'); d$group <- factor(d$group)
S <- Surv(d$time, d$cens)
wb <- survreg(S ~ group, data = d, dist = 'weibull')
sigma <- wb$scale; k <- 1/sigma
cat(sprintf('Weibull shape k = %.3f\n', k))
for (g in 2:4) {
  b_aft <- coef(wb)[paste0('group', g)]; hr <- exp(-b_aft/sigma)
  cat(sprintf('  group %d: HR vs g1 = %.2f\n', g, hr))
}
# median survival per group from the AFT predictor
md <- predict(wb, newdata = data.frame(group = factor(1:4)), type = 'quantile', p = 0.5)
cat('median survival (weeks):', paste(sprintf('g%d=%.1f', 1:4, md), collapse='  '), '\n')

Weibull shape k = 3.254

  group 2: HR vs g1 = 0.31
  group 3: HR vs g1 = 0.71
  group 4: HR vs g1 = 1.48

median survival (weeks): g1=24.2  g2=34.7  g3=26.9  g4=21.4 

km <- survfit(S ~ group, data = d); print(km)
plot(km, col = c('steelblue','seagreen','goldenrod','firebrick'), lwd = 2, xlab='weeks', ylab='S(t)',
     main = 'Kaplan-Meier by group (mice)')
legend('topright', legend = paste('group', 1:4), col = c('steelblue','seagreen','goldenrod','firebrick'), lwd=2, bty='n')

Call: survfit(formula = S ~ group, data = d)

         n events median 0.95LCL 0.95UCL
group=1 20     19     24      21      30
group=2 20     13     32      29      NA
group=3 20     16     26      22      35
group=4 20     17     22      18      28

Results

survreg matches the Bayesian fits: k ≈ 3.25 (increasing hazard), group 2 longest-lived (HR ≈ 0.3, median ≈ 35 wk), group 4 shortest. Three engines agree; the conclusion (group 2 best, strong wear-out shape) is identical to the Python/PyMC fits.