Interval-censored survival — PyMC cross-check

Companion to icsurv_python.ipynb. The interval likelihood $S(l)-S(u)$ is encoded directly with pm.Potential; NUTS samples $(\beta,\log k)$.

import numpy as np, pandas as pd, pymc as pm, pytensor.tensor as pt, arviz as az
from icsurv import icweib_rwm
d = pd.read_csv('mice_tumor_ic.csv'); ge=(d['grp']=='ge').astype(float).to_numpy()
X = np.column_stack([np.ones(len(d)), ge]); l=d['l'].to_numpy(float); u=d['u'].to_numpy(float)
uinf = np.isinf(u) | (u > 1e8)
SC = 100.0                                                  # rescale time (days/100) so the intercept is O(1) -> better NUTS geometry
ls = l / SC; us_safe = np.where(uinf, 1.0, u / SC)
with pm.Model() as m:
    beta = pm.Normal('beta', 0, 20, shape=2); logk = pm.Normal('logk', 0, 5); k = pm.math.exp(logk)
    lam = pm.math.exp(pt.dot(X, beta))
    Sl = pm.math.exp(-lam * ls**k)
    Su = pt.switch(uinf, 0.0, pm.math.exp(-lam * us_safe**k))
    pm.Potential('lik', pt.sum(pt.log(Sl - Su)))
    pm.Deterministic('HR', pm.math.exp(beta[1])); pm.Deterministic('k', k)
    idata = pm.sample(2000, tune=2000, chains=4, target_accept=0.99, random_seed=1, progressbar=False)
po=idata.posterior; hr=po['HR'].values.ravel()
print('PyMC: germ-free HR %.2f [%.2f,%.2f]  P(HR>1) %.3f  k %.3f  max r_hat %.3f'
      %(hr.mean(),np.percentile(hr,2.5),np.percentile(hr,97.5),(np.log(hr)>0).mean(),float(po['k'].mean()),
        float(az.summary(idata,var_names=['beta','logk'])['r_hat'].max())))
g=icweib_rwm(X,l,np.where(uinf,np.inf,u),np.zeros(len(d),bool),seed=1)['draws']
print('from-scratch: HR %.2f  k %.3f  (survreg/icenReg: HR ~2.2, k ~2.03)'%(np.exp(g[:,1].mean()),np.exp(g[:,2].mean())))

Initializing NUTS using jitter+adapt_diag...

Multiprocess sampling (4 chains in 4 jobs)

NUTS: [beta, logk]

Sampling 4 chains for 2_000 tune and 2_000 draw iterations (8_000 + 8_000 draws total) took 15 seconds.

PyMC: germ-free HR 2.55 [1.17,4.82]  P(HR>1) 0.994  k 1.736  max r_hat 1.010

from-scratch: HR 2.27  k 1.787  (survreg/icenReg: HR ~2.2, k ~2.03)

Results

PyMC (interval likelihood via pm.Potential) reproduces the from-scratch fit. It reports the germ-free tumour HR ≈ 2.5 as the posterior mean of $e^\beta$ — the from-scratch and survreg quote the plug-in $e^{\bar\beta}\approx 2.2$, a touch lower because the HR posterior is right-skewed. The Bayesian shape k ≈ 1.74 (PyMC and from-scratch agree) sits a little below the MLE/survreg value of 2.03; r̂ ≈ 1.0. The $S(l)-S(u)$ contribution is all that's needed; pt.switch sets $S(u)=0$ for the right-censored mice.