Interval-censored survival — R cross-check (survreg, icenReg)¶

Companion to icsurv_python.ipynb / icsurv_pymc.ipynb. survival::survreg fits a parametric Weibull with Surv(l, u, type='interval2'); icenReg::ic_par is the dedicated interval-censored parametric tool.

In [1]:
suppressMessages(library(survival))
d <- read.csv('mice_tumor_ic.csv')
d$u2 <- ifelse(is.infinite(d$u) | d$u > 1e8, NA, d$u)            # NA upper = right-censored
fit <- survreg(Surv(l, u2, type='interval2') ~ grp, data=d, dist='weibull')
k <- 1/fit$scale; hr <- exp(-coef(fit)[['grpge']]/fit$scale)
cat(sprintf('survreg  : shape k = %.3f   germ-free PH HR = %.2f\n', k, hr))
print(round(summary(fit)$table, 3))

survreg  : shape k = 2.028   germ-free PH HR = 2.19

             Value Std. Error      z     p
(Intercept)  6.948      0.211 32.937 0.000
grpge       -0.388      0.281 -1.379 0.168
Log(scale)  -0.707      0.382 -1.852 0.064
In [2]:
ok <- suppressWarnings(suppressMessages(require(icenReg)))
if (!ok) { install.packages('icenReg', repos='https://cloud.r-project.org'); library(icenReg) }
fit2 <- ic_par(cbind(l, u) ~ grp, data=transform(read.csv('mice_tumor_ic.csv'), u=ifelse(is.infinite(u)|u>1e8, Inf, u)),
               model='ph', dist='weibull')
cat('icenReg ic_par (Weibull PH):\n'); print(summary(fit2))

icenReg ic_par (Weibull PH):

Model:  Cox PH
Dependency structure assumed: Independence
Baseline:  weibull 
Call: ic_par(formula = cbind(l, u) ~ grp, data = transform(read.csv("mice_tumor_ic.csv"), 
    u = ifelse(is.infinite(u) | u > 1e+08, Inf, u)), model = "ph", 
    dist = "weibull")

          Estimate Exp(Est) Std.Error z-value       p
log_shape   0.7072    2.028    0.3820   1.852 0.06410
log_scale   6.8190  915.000    0.1268  53.760 0.00000
grpge       0.7862    2.195    0.3359   2.341 0.01925

final llk =  -80.3202 
Iterations =  16
In [3]:
# Turnbull NPMLE -- nonparametric (KM-analog) survival for interval-censored data (rendered inline)
dd <- transform(read.csv('mice_tumor_ic.csv'), u = ifelse(is.infinite(u) | u > 1e8, Inf, u))
np <- ic_np(cbind(l, u) ~ grp, data = dd)
options(repr.plot.width = 7, repr.plot.height = 5)
plot(np, col = c('steelblue', 'firebrick'), lwd = 2, xlab = 'days', ylab = 'S(t)',
     main = 'Turnbull NPMLE survival by group (interval-censored)')
legend('bottomleft', legend = c('conventional (ce)', 'germ-free (ge)'),
       col = c('steelblue', 'firebrick'), lwd = 2, bty = 'n')
No description has been provided for this image

Results¶

  • survreg (interval2) gives k ≈ 2.03 and germ-free HR ≈ 2.2 — matching the from-scratch fit.
  • icenReg::ic_par (the purpose-built interval-censored Weibull PH) gives the same hazard ratio.
  • Engines agree on the effect: from-scratch, survreg, and icenReg all give germ-free HR ~2.2×, with PyMC's posterior-mean HR a touch higher (~2.5, averaging $e^\beta$ over a right-skewed posterior). The MLE shape is k ≈ 2.03, the Bayesian fits ~1.74. Estimated honestly from the left- and right-censored onset times via the same $S(l)-S(u)$ likelihood.