Turnbull’s self-consistent NPMLE for interval-censored survival data.
Usage
Kaplan-Meier requires knowing each subject’s event time exactly (up to right-censoring). When follow-up is periodic instead (clinic visits, equipment inspections), all that is known is that the event fell in a window (lower, upper]. Turnbull’s nonparametric maximum likelihood estimator handles this directly, together with mixtures of exact, left-, right-, and interval-censored observations in the same fit.
Unlike Kaplan-Meier, the NPMLE’s support is not identified everywhere: the data determine a set of maximal intersection intervals, and only the total probability mass on each one is identified, not its placement inside. Turnbull reports both bounds of every such interval. Where an interval degenerates to a single point (an exact death, or a region every subject’s constraint resolves unambiguously), survival is known exactly there. Otherwise it is genuinely unidentified and predict()/quantile() return nan rather than interpolate. rmst()/rmrl() need a single number, so they instead fall back to Turnbull’s own right-endpoint convention (every atom’s mass resolves at interval_high_); see their docstrings for the resulting conservative (upper-bound) bias.
To use this estimator, call fit() with a Surv response built via Surv.interval() (the general case), or Surv.right()/Surv.left() (degenerate cases: fitting a right-censored response through Turnbull reproduces KaplanMeier exactly, since there is then no genuine interval ambiguity). Left-truncated (Surv.counting()) and multi-state responses are not supported.
tol: float = 1e-09Convergence tolerance on the largest change in any atom’s probability mass between EM iterations (default 1e-9).
max_iter: int = 1000010000).
Fitted estimatorfit() to produce a fitted estimator with cached results (interval_low_, interval_high_, prob_mass_, survival_), accessible as aligned arrays or exported to DataFrames.
Call fit with a Surv response. The maximal intersection intervals are found via the Gentleman & Geyer (1994) construction: the finest partition induced by every subject’s lower/upper bounds is grouped into runs of consecutive atoms that every subject’s constraint either fully includes or fully excludes (only their combined mass is identified). Probabilities on that support are then found by the EM self-consistency algorithm, which is monotone in the likelihood and converges to the NPMLE.
Six subjects observed at irregular follow-up visits, so some events are only known to lie in a window rather than at an exact time. Two are exact deaths, one is right-censored (no event observed by the last visit), and three are genuinely interval-censored:
import greenwood as gw
# Build an interval-censored response: (lower, upper] windows, inf marks right-censoring
y = gw.Surv.interval(
lower=[0, 4, 7, 0, 3, 5],
upper=[4, float("inf"), 7, 2.5, 6, 5],
)
# Fit the Turnbull NPMLE
tb = gw.Turnbull().fit(y)
tbTurnbull (self-consistent NPMLE for interval-censored data)
n atoms ambiguous iters converged
6 8 5 32 TrueThe fitted curve, one row per maximal intersection interval, is available via to_frame:
| Name | Description |
|---|---|
| interval_high_ | Upper bound of each maximal intersection interval. |
| interval_low_ | Lower bound of each maximal intersection interval. |
| prob_mass_ | Probability mass assigned to each maximal intersection interval. |
| strata_ |
Stratum labels for each row, or None for unstratified fits.
|
| survival_ |
Survival estimate just after each interval (1 - cumsum(prob_mass_)).
|
Upper bound of each maximal intersection interval.
interval_high_: Array
Lower bound of each maximal intersection interval.
interval_low_: Array
Probability mass assigned to each maximal intersection interval.
prob_mass_: Array
Stratum labels for each row, or None for unstratified fits.
strata_: Array | None
Survival estimate just after each interval (1 - cumsum(prob_mass_)).
survival_: Array
| Name | Description |
|---|---|
| fit() | Fit Turnbull’s NPMLE to interval-censored survival data. |
| median() | Median survival time per stratum (the 0.5-quantile). |
| predict() | Evaluate the survival curve at specified times. |
| quantile() |
Return the p-quantile survival time per stratum.
|
| rmrl() |
Restricted mean residual life at s, under the right-endpoint convention.
|
| rmst() |
Restricted mean survival time up to tau, under the right-endpoint convention.
|
| to_frame() | Return the fitted NPMLE as a DataFrame. |
Fit Turnbull’s NPMLE to interval-censored survival data.
Usage
Computes the maximal intersection intervals and their probability masses from a Surv response, via EM self-consistency. Pass by= to fit separate curves per group (stratified analysis).
surv: SurvA Surv response built with Surv.interval(), Surv.right(), or Surv.left(). Left-truncated (Surv.counting()) and multi-state responses raise NotImplementedError.
by: Any = NoneOptional grouping variable (e.g., a column or array). Produces one fit per unique value of by. Default (None): a single, unstratified fit.
weights: Any = Nonesurv. Default (None): uses surv.weights if present, otherwise unit weights.
Turnbullsurvival_, …) as attributes.
Median survival time per stratum (the 0.5-quantile).
Usage
A convenience wrapper around quantile(0.5). See quantile for the return shape.
tuple or dict(estimate, lower, upper) for a single stratum, or a dict keyed by stratum label for stratified fits.
Evaluate the survival curve at specified times.
Usage
times: Any ndarray or dicttimes’ shape) of survival estimates, with nan at any query time that falls strictly inside a non-degenerate (ambiguous) interval. For stratified fits: a dict keyed by stratum label.
Outside every ambiguous interval, the curve is a well-defined, right-continuous step function, exactly as for KaplanMeier. Strictly inside a non-degenerate interval (lower, upper), the true survival value is not identified by the data, so nan is returned there rather than an arbitrary interpolated guess.
Return the p-quantile survival time per stratum.
Usage
p: floatp=0.5 for the median).
tuple or dict(estimate, lower, upper). estimate is nan when the crossing falls inside a non-degenerate (ambiguous) interval, in which case lower/upper are that interval’s bounds rather than a sampling confidence bound. Otherwise estimate == lower == upper. For stratified fits: a dict keyed by stratum label, with values as above.
The quantile is found by inverting the step-function survival curve, exactly as for KaplanMeier, but returning the identifiability bracket rather than ever guessing a point inside it. This is always a 3-tuple regardless of whether the crossing is ambiguous, so the return type does not depend on the data.
Restricted mean residual life at s, under the right-endpoint convention.
Usage
s: floatThe landmark time. Must be non-negative.
tau: floats.
float or dictdict keyed by stratum label for stratified fits. nan if everyone has resolved (under the same right-endpoint convention) by time s.
Generalizes rmst to a later landmark: \mathrm{RMRL}(s; \tau) = \int_s^\tau S(u)\,du / S(s), under the same right-endpoint convention (every atom’s mass is treated as resolving at interval_high_); see rmst for why, and for the resulting conservative (upper-bound) bias.
Restricted mean survival time up to tau, under the right-endpoint convention.
Usage
tau: float float or dictdict keyed by stratum label for stratified fits.
RMST is the area under the survival curve on [0, \tau]. Unlike predict() and quantile(), which report the genuine identifiability gap as nan, RMST needs a single number, so every atom’s probability mass, ambiguous or not, is treated as resolving exactly at that atom’s right endpoint (interval_high_). This is Turnbull’s own convention for reporting a plottable curve from an otherwise partially-unidentified NPMLE. It is also the most conservative choice for RMST: placing mass as late as possible maximizes the area under the curve, so this systematically reports the largest RMST consistent with the data, not an unbiased point estimate. There is no variance estimator for it (no confidence interval is returned).
Return the fitted NPMLE as a DataFrame.
Usage
Exports one row per maximal intersection interval, with its bounds, probability mass, and the survival estimate just after it.
format: str | None = NoneNone (default), "pandas", "polars", or "pyarrow". When None, a backend is auto-detected (Polars, then Pandas, then PyArrow).
pandas.DataFrame, polars.DataFrame, or pyarrow.Tableinterval_low, interval_high, prob_mass, estimate, and optionally strata.
ImportErrorimport greenwood as gw
y = gw.Surv.interval(lower=[0, 4, 7, 0, 3, 5], upper=[4, float("inf"), 7, 2.5, 6, 5])
tb = gw.Turnbull().fit(y)
tb.to_frame(format="polars")| interval_low | interval_high | prob_mass | estimate |
|---|---|---|---|
| f64 | f64 | f64 | f64 |
| 0.0 | 2.5 | 0.230328 | 0.769672 |
| 2.5 | 3.0 | 2.2313e-19 | 0.769672 |
| 3.0 | 4.0 | 8.4788e-10 | 0.769672 |
| 4.0 | 4.0 | 0.372678 | 0.396994 |
| 4.0 | 5.0 | 1.7247e-19 | 0.396994 |
| 5.0 | 5.0 | 0.230328 | 0.166667 |
| 5.0 | 6.0 | 1.7247e-19 | 0.166667 |
| 7.0 | 7.0 | 0.166667 | 1.1102e-16 |