Interval-censored survival — Weibull from scratch

When the event time is only known to lie in an interval

Module: icsurv.py. The sequel to survival_python.ipynb (right-censoring only). In many studies subjects are checked periodically, so an event is known only to have happened between two visits — interval censoring. We fit a Weibull proportional-hazards model to the classic mouse lung-tumour data (Hoel & Walburg 1972; from R's icenReg), comparing two housing environments.

One likelihood for every censoring type

The event time lies in $(l_i, u_i]$. With $S(t)=e^{-\lambda_i t^{k}}$, $\lambda_i=e^{x_i'\beta}$, each observation contributes $$\Pr(l_i<T\le u_i)=S(l_i)-S(u_i),$$ which covers all cases at once:

type	bounds	contribution
interval	$l<u<\infty$	$S(l)-S(u)$
left	$l=0$	$S(0)-S(u)=F(u)$ — event before first visit
right	$u=\infty$	$S(l)-0=S(l)$ — no event by last visit
exact	$l=u$	density $f(l)=\lambda k l^{k-1}S(l)$

Algorithm

• MLE: maximise $\sum_i \log[S(l_i)-S(u_i)]$ over $(\beta,\log k)$ (computed with a numerically stable log(S(l)-S(u))).
• Bayes: RW-Metropolis on $(\beta,\log k)$, proposal from the MLE Hessian; $e^{\beta}$ is the hazard ratio.

(A NEW module — survival_mcmc.py handled right-censoring only — leaving the older project untouched.)

import numpy as np, pandas as pd, matplotlib.pyplot as plt
from icsurv import simulate_icweib, icweib_mle, icweib_rwm, summary

# ---- Section 1: synthetic recovery (interval + left + right all present) ----
l, u, X, exact, T = simulate_icweib(n=2000, beta=(-7.0, 0.8), k=2.0, seed=1)
print('synthetic: %d obs  | interval %d, right-censored %d, left (l=0) %d'
      % (len(l), int(((l>0)&np.isfinite(u)).sum()), int(np.isinf(u).sum()), int((l==0).sum())))
m, se = icweib_mle(X, l, u, exact)
o = icweib_rwm(X, l, u, exact, seed=1)
print('  truth: beta=(-7.0, 0.8), k=2.0')
print('  MLE  : beta=%s  k=%.3f' % (np.round(m[:2],3), np.exp(m[2])))
print('  Bayes: grp coef %.3f [%.2f, %.2f]  HR %.2f  k %.3f  (acc %.2f)'
      % (o['draws'][:,1].mean(), np.percentile(o['draws'][:,1],2.5), np.percentile(o['draws'][:,1],97.5),
         np.exp(o['draws'][:,1].mean()), np.exp(o['draws'][:,2].mean()), o['accept']))

synthetic: 2000 obs  | interval 1469, right-censored 373, left (l=0) 158

  truth: beta=(-7.0, 0.8), k=2.0
  MLE  : beta=[-7.253  0.825]  k=2.061
  Bayes: grp coef 0.826 [0.73, 0.93]  HR 2.28  k 2.062  (acc 0.64)

Section 2 — mouse lung-tumour data (interval-censored)

144 male mice (Hoel & Walburg 1972) in two housing environments — conventional (ce, 96) vs germ-free (ge, 48) — examined for lung tumours. The tumour-onset time is interval-censored: l and u bracket the unknown onset (in days). There are no exact times: every observation is either left-censored (l=0, tumour already present) or right-censored (u=∞, no tumour by death).

field	meaning
l	lower bound of the onset interval (0 ⇒ tumour already present at first look = left-censored)
u	upper bound (∞ ⇒ no tumour by the end = right-censored)
grp	environment: ce (conventional) or ge (germ-free)

(of the 144: 62 left-censored and 82 right-censored — no exactly-observed onset times, and in fact no finite-width intervals in this sample; the same $S(l)-S(u)$ likelihood still handles the left/right mix in one stroke, including the left-censored cases a plain right-censored Weibull can't use.)

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)
u = np.where((u > 1e8) | np.isinf(u), np.inf, u); exact = np.zeros(len(d), bool)
print('interval %d, right-censored %d, left (l=0) %d' % (int(((l>0)&np.isfinite(u)).sum()), int(np.isinf(u).sum()), int((l==0).sum())))
m, se = icweib_mle(X, l, u, exact)
o = icweib_rwm(X, l, u, exact, seed=1); g = o['draws'][:,1]
print('\nWeibull (interval-censored):')
print('  shape k = %.3f' % np.exp(m[2]))
print('  germ-free vs conventional: PH coef %.3f [%.2f, %.2f]  HR %.2f  P(HR>1) %.3f'
      % (g.mean(), np.percentile(g,2.5), np.percentile(g,97.5), np.exp(g.mean()), (g>0).mean()))
print('  (survreg interval2 benchmark: k 2.03, germ-free HR ~2.2)')

interval 0, right-censored 82, left (l=0) 62

Weibull (interval-censored):
  shape k = 2.028

  germ-free vs conventional: PH coef 0.819 [0.23, 1.42]  HR 2.27  P(HR>1) 0.998
  (survreg interval2 benchmark: k 2.03, germ-free HR ~2.2)

# fitted Weibull survival by group + NAIVE midpoint Kaplan-Meier (which mishandles interval censoring)
def km(times, evf):
    o = np.argsort(times); times, evf = times[o], evf[o]; S = 1.0; ts=[0.]; Ss=[1.]
    for ti in np.unique(times[evf==1]):
        S *= 1 - ((times==ti)&(evf==1)).sum()/(times>=ti).sum(); ts.append(ti); Ss.append(S)
    return np.array(ts), np.array(Ss)
bm = m; kk = np.exp(bm[2]); grid = np.linspace(1, 1000, 250)
fig, ax = plt.subplots(1, 2, figsize=(12, 4.6))
for gval, lab, col in [(0,'conventional (ce)','steelblue'), (1,'germ-free (ge)','firebrick')]:
    msk = ge==gval
    mid = np.where(np.isinf(u[msk]), l[msk], 0.5*(l[msk]+u[msk]))      # naive midpoint / last-seen
    evn = (~np.isinf(u[msk])).astype(float)                            # interval -> 'event', right-cens -> censored
    tk, Sk = km(mid, evn); ax[0].step(tk, Sk, where='post', color=col, lw=1.4, alpha=.6, label='midpoint KM '+lab)
    lam = np.exp(np.array([1.0, gval]) @ bm[:2]); ax[0].plot(grid, np.exp(-lam*grid**kk), '--', color=col, lw=1.8)
ax[0].set_xlabel('days'); ax[0].set_ylabel('S(t)'); ax[0].set_ylim(0,1.02)
ax[0].set_title('Fitted Weibull (dashed) vs naive midpoint KM (step)'); ax[0].legend(fontsize=7)
hr = np.exp(o['draws'][:,1]); p_hr = (o['draws'][:,1] > 0).mean()
ax[1].hist(hr, bins=50, density=True, color='seagreen', alpha=.8)
ax[1].axvline(hr.mean(), color='k', ls='--', lw=1, label=f'mean {hr.mean():.2f}')
ax[1].axvline(1, color='firebrick', ls=':', lw=1, label='no effect')
ax[1].set_xlabel('hazard ratio (germ-free vs conventional)'); ax[1].set_yticks([])
ax[1].set_title(f'Tumour-hazard effect  (P(HR>1) = {p_hr:.3f})'); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.savefig('icsurv.png', dpi=120, bbox_inches='tight'); plt.show()

Clarifying views: the data, the uncertainty, and the "germ-free puzzle"

• The interval data itself — each mouse is a bracket $[l,u]$, not a point (right-censored mice = arrows): makes interval censoring concrete and shows the per-group timing.
• Cumulative tumour incidence $F(t)=1-S(t)$ with 95% credible bands — germ-free reach a higher tumour probability, with uncertainty.
• Why is germ-free the higher-hazard group? Not because sterility causes cancer — the standard explanation is competing risks / longevity: conventional mice die earlier of infection, before a slow lung tumour appears, so fewer are ever observed with one; germ-free mice live long enough for the (age-related) tumour to show up. The panels below give the tumour-observed split by group and the posterior median onset (earlier for germ-free, consistent with the higher hazard). Caveat: the model estimates only the tumour-onset hazard — it can't, by itself, disentangle "more tumour biology" from "less competing death."

# --- the data itself (interval brackets) + cumulative incidence with credible bands ---
fig, ax = plt.subplots(1, 2, figsize=(13, 5))
# (A) interval-segment plot: each mouse [l, u]; right-censored (u=inf) as an arrow
cap = float(np.nanmax(np.where(np.isinf(u), np.nan, u))) * 1.05
order = np.lexsort((l, ge))
for rank, i in enumerate(order):
    col = 'firebrick' if ge[i] == 1 else 'steelblue'
    if np.isinf(u[i]):
        ax[0].plot([l[i], cap], [rank, rank], color=col, lw=0.6, alpha=.6)
        ax[0].plot([cap], [rank], '>', color=col, ms=3)
    else:
        ax[0].plot([l[i], u[i]], [rank, rank], color=col, lw=0.8, alpha=.75)
ax[0].plot([], [], color='steelblue', lw=2, label='conventional (ce)')
ax[0].plot([], [], color='firebrick', lw=2, label='germ-free (ge)')
ax[0].set_xlabel('days'); ax[0].set_ylabel('mouse (sorted: group, then l)')
ax[0].set_title('The interval-censored data\n(bar = onset bracket [l,u];  > = right-censored)'); ax[0].legend(fontsize=8)
# (B) cumulative tumour incidence F(t)=1-S(t) by group, 95% credible band
D = o['draws']; sel = D[np.linspace(0, len(D)-1, 300).astype(int)]; grid = np.linspace(1, 1000, 200)
for gv, lab, col in [(0,'conventional (ce)','steelblue'), (1,'germ-free (ge)','firebrick')]:
    lam = np.exp(sel[:,0] + sel[:,1]*gv); kk = np.exp(sel[:,2])
    F = 1 - np.exp(-lam[:,None] * grid[None,:]**kk[:,None])
    ax[1].plot(grid, F.mean(0), color=col, lw=2, label=lab)
    ax[1].fill_between(grid, np.percentile(F,2.5,0), np.percentile(F,97.5,0), color=col, alpha=.2)
ax[1].set_xlabel('days'); ax[1].set_ylabel('P(tumour onset by t)'); ax[1].set_ylim(0, 1)
ax[1].set_title('Cumulative tumour incidence (95% credible band)'); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.savefig('icsurv_data.png', dpi=120, bbox_inches='tight'); plt.show()

# --- why is germ-free HIGHER hazard? competing-risks view + posterior median onset ---
fig, ax = plt.subplots(1, 2, figsize=(12, 4.6)); gn = ['conventional (ce)', 'germ-free (ge)']
tum = np.array([((ge==gv) & np.isfinite(u)).sum() for gv in (0,1)], float)
notum = np.array([((ge==gv) & np.isinf(u)).sum() for gv in (0,1)], float); tot = tum + notum
xx = np.arange(2)
ax[0].bar(xx, tum, 0.6, label='tumour observed', color='firebrick', alpha=.85)
ax[0].bar(xx, notum, 0.6, bottom=tum, label='no tumour by death (right-cens.)', color='lightgray', edgecolor='k')
for i in range(2): ax[0].text(i, tot[i]+1.5, f'{tum[i]/tot[i]:.0%} with tumour', ha='center', fontsize=8)
ax[0].set_xticks(xx); ax[0].set_xticklabels(gn, fontsize=8); ax[0].set_ylabel('number of mice')
ax[0].set_title('Tumour observed vs not, by group'); ax[0].legend(fontsize=7)
D = o['draws']; med = {}
for gv in (0,1):
    lam = np.exp(D[:,0] + D[:,1]*gv); kk = np.exp(D[:,2]); med[gv] = (np.log(2)/lam)**(1/kk)
ax[1].errorbar(xx, [med[0].mean(), med[1].mean()],
    yerr=[[med[g].mean()-np.percentile(med[g],2.5) for g in (0,1)],
          [np.percentile(med[g],97.5)-med[g].mean() for g in (0,1)]], fmt='o', color='black', capsize=5, ms=8)
ax[1].set_xticks(xx); ax[1].set_xticklabels(gn, fontsize=8); ax[1].set_xlim(-.5, 1.5)
ax[1].set_ylabel('median tumour-onset (days)'); ax[1].set_title('Posterior median onset by group (95% CrI)')
plt.tight_layout(); plt.savefig('icsurv_why.png', dpi=120, bbox_inches='tight'); plt.show()

Results

• Germ-free mice have a higher tumour hazard: HR ≈ 2.2 (germ-free vs conventional), P(HR>1) ≈ 0.99 — they develop lung tumours faster. Weibull shape k ≈ 2 (increasing hazard).
• The fit matches the survreg/icenReg interval-censored benchmark (k ≈ 2.03, HR ≈ 2.2).
• The dashed fitted curves sit close to the naive midpoint Kaplan-Meier here, but midpoint imputation is a shortcut — it invents exact event times and is known to bias estimates when intervals are wide; the interval likelihood uses $S(l)-S(u)$ honestly.

Takeaways

• One formula, every censoring type. $S(l)-S(u)$ collapses to right-censoring ($u=\infty$), left-censoring ($l=0$), interval, and (with the density) exact — so the same code spans the whole censoring spectrum.
• This completes the censoring tour started with right-censored survival and the Tobit (left/interval Gaussian): censoring is one idea — a probability mass over the region the value is known to lie in.
• Cross-checks: PyMC (interval likelihood via pm.Potential) and R (survreg(type='interval2'), icenReg::ic_par) follow.