Turnbull

Turnbull’s self-consistent NPMLE for interval-censored survival data.

Usage

Turnbull(
    *,
    tol=1e-09,
    max_iter=10000,
)

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 an interval-censored Surv response, built with gw.Surv(time=lower, time2=upper, type="interval2") or from a gw.event_time() vector via gw.as_surv(). With a data frame, name its columns in a formula such as "Surv(lower, upper, type='interval2')" and pass data=. Right- and left-censored responses are degenerate cases. Fitting a right-censored response through Turnbull reproduces KaplanMeier exactly, since there is then no genuine interval ambiguity. Left-truncated (counting-process) and multi-state responses are not supported.

Parameters

tol: float = 1e-09: Convergence tolerance on the largest change in any atom’s probability mass between EM iterations (default 1e-9).
max_iter: int = 10000: Maximum number of EM iterations (default 10000).

Returns

Fitted estimator: Call fit() to produce a fitted estimator with cached results (interval_low_, interval_high_, prob_mass_, survival_), accessible as aligned arrays or exported to DataFrames.

Details

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.

Examples

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(
    time=[0, 4, 7, 0, 3, 5],
    time2=[4, float("inf"), 7, 2.5, 6, 5],
    type="interval2",
)

# Fit the Turnbull NPMLE
tb = gw.Turnbull().fit(y)
tb

Turnbull (self-consistent NPMLE for interval-censored data)

     n  atoms  ambiguous  iters  converged
all  6      8          6    120       True

The fitted curve, one row per maximal intersection interval, is available via to_frame:

tb.to_frame(format="polars")

shape: (8, 5)

strata	interval_low	interval_high	prob_mass	estimate
str	f64	f64	f64	f64
"all"	0.0	2.5	0.333333	0.666667
"all"	2.5	3.0	4.0750e-38	0.666667
"all"	3.0	4.0	6.3494e-9	0.666667
"all"	4.0	5.0	2.9652e-26	0.666667
"all"	5.0	5.0	0.444444	0.222222
"all"	5.0	6.0	2.9652e-26	0.222222
"all"	6.0	7.0	9.1637e-74	0.222222
"all"	7.0	7.0	0.222222	1.1102e-16

Attributes

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_	The curve label of each row: `"all"` for one curve, `"sex=1"` and so on for groups.
survival_	Survival estimate just after each interval (`1 - cumsum(prob_mass_)`).

interval_high_

Upper bound of each maximal intersection interval.

interval_high_: Array

interval_low_

Lower bound of each maximal intersection interval.

interval_low_: Array

prob_mass_

Probability mass assigned to each maximal intersection interval.

prob_mass_: Array

strata_

The curve label of each row: "all" for one curve, "sex=1" and so on for groups.

strata_: Array

survival_

Survival estimate just after each interval (1 - cumsum(prob_mass_)).

survival_: Array

Methods

Name	Description
fit()	Fit Turnbull’s NPMLE to interval-censored survival data.
median()	Median survival time, with its identifiability bounds, for each curve.
predict()	Read the survival curve at given times, for each curve.
quantile()	Return survival-time quantiles, with their identifiability bounds, for each curve.
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()

Fit Turnbull’s NPMLE to interval-censored survival data.

Usage

Source

fit(
    surv,
    *,
    data=None,
    by=None,
    weights=None,
)

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).

Parameters

surv: Surv | Outcome | str: An interval-, right-, or left-censored Surv response built with Surv(). Left-truncated (counting-process) and multi-state responses raise NotImplementedError. An Outcome or a formula string such as "Surv(lower, upper, type='interval2') ~ sex" is also accepted, with its columns read from data. The right-hand side names the by column(s).
by: Any = None: Optional grouping variable (e.g., a column or array). Produces one fit per unique value of by. Default (None): a single, unstratified fit.
weights: Any = None: Optional case weights. Must have the same length as surv. Default (None): unit weights.
data: Any = None: A data frame (pandas, Polars, PyArrow, DuckDB, a lazy frame, …) holding the columns named by the response, by, and weights. When surv is an Outcome or a formula, rows with a missing value in any column used are dropped before fitting.

Returns

Turnbull: The fitted estimator itself (for method chaining), with cached results (interval_low_, interval_high_, prob_mass_, survival_, …) as attributes.

Examples

import greenwood as gw

y = gw.Surv(
    time=[0, 4, 7, 0, 3, 5],
    time2=[4, float("inf"), 7, 2.5, 6, 5],
    type="interval2",
)
tb = gw.Turnbull().fit(y)
tb.survival_

array([6.66666670e-01, 6.66666670e-01, 6.66666664e-01, 6.66666664e-01,
       2.22222223e-01, 2.22222223e-01, 2.22222223e-01, 1.11022302e-16])

median()

Median survival time, with its identifiability bounds, for each curve.

Usage

Source

median(
    *,
    format=None,
)

Parameters

format: str | None = None: Output format: None (default), "pandas", "polars", or "pyarrow". When None, a backend is auto-detected (Polars, then Pandas, then PyArrow).

Returns

pandas.DataFrame, polars.DataFrame, or pyarrow.Table: One row per curve, with columns strata, time, time_low, and time_high. See quantile() for what the bounds mean.

Examples

import greenwood as gw

y = gw.Surv(
    time=[0, 4, 7, 0, 3, 5],
    time2=[4, float("inf"), 7, 2.5, 6, 5],
    type="interval2",
)
tb = gw.Turnbull().fit(y)
tb.median(format="polars")

shape: (1, 4)

strata	time	time_low	time_high
str	f64	f64	f64
"all"	5.0	5.0	5.0

predict()

Read the survival curve at given times, for each curve.

Usage

Source

predict(
    times,
    *,
    format=None,
)

Parameters

times: Any: One or more times at which to read the curve.
format: str | None = None: Output format: None (default), "pandas", "polars", or "pyarrow". When None, a backend is auto-detected (Polars, then Pandas, then PyArrow).

Returns

pandas.DataFrame, polars.DataFrame, or pyarrow.Table: One row per curve and time, with columns strata, time, and estimate. The estimate is nan at any time strictly inside a non-degenerate (ambiguous) interval. The layout is the same for one curve (strata is "all") as for several.

Details

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.

Examples

import greenwood as gw

y = gw.Surv(
    time=[0, 4, 7, 0, 3, 5],
    time2=[4, float("inf"), 7, 2.5, 6, 5],
    type="interval2",
)
tb = gw.Turnbull().fit(y)

# nan at t=1 (strictly inside the ambiguous (0, 2.5) region)
tb.predict(times=[1, 2.5, 5, 7], format="polars")

shape: (4, 3)

strata	time	estimate
str	f64	f64
"all"	1.0	NaN
"all"	2.5	0.666667
"all"	5.0	0.222222
"all"	7.0	1.1102e-16

quantile()

Return survival-time quantiles, with their identifiability bounds, for each curve.

Usage

Source

quantile(
    p,
    *,
    format=None,
)

Parameters

p: Any: One or more proportions between 0 and 1 (for example, 0.5 for the median).
format: str | None = None: Output format: None (default), "pandas", "polars", or "pyarrow". When None, a backend is auto-detected (Polars, then Pandas, then PyArrow).

Returns

pandas.DataFrame, polars.DataFrame, or pyarrow.Table: One row per curve and proportion, with columns strata, prob, time, time_low, and time_high. When the crossing falls inside a non-degenerate (ambiguous) interval, time is nan and time_low/time_high are that interval’s bounds. They bound what the data can identify and are not confidence limits. Otherwise all three are equal. The layout is the same for one curve (strata is "all") as for several.

Details

The quantile is found by inverting the step-function survival curve, exactly as for KaplanMeier, but reporting the identifiability bracket rather than guessing a point inside it.

Examples

import greenwood as gw

y = gw.Surv(
    time=[0, 4, 7, 0, 3, 5],
    time2=[4, float("inf"), 7, 2.5, 6, 5],
    type="interval2",
)
tb = gw.Turnbull().fit(y)

# The first-quartile survival time (or its ambiguity bracket)
tb.quantile(p=0.25, format="polars")

shape: (1, 5)

strata	prob	time	time_low	time_high
str	f64	f64	f64	f64
"all"	0.25	NaN	0.0	2.5

rmrl()

Restricted mean residual life at s, under the right-endpoint convention.

Usage

Source

rmrl(
    s,
    tau,
    *,
    format=None,
)

Parameters

s: float: The landmark time. Must be non-negative.
tau: float: The end of the window. Must be greater than s.
format: str | None = None: Output format: None (default), "pandas", "polars", or "pyarrow". When None, a backend is auto-detected (Polars, then Pandas, then PyArrow).

Returns

pandas.DataFrame, polars.DataFrame, or pyarrow.Table: One row per curve, with columns strata, s, tau, and estimate. The estimate is nan for a curve where everyone has resolved (under the same right-endpoint convention) by s.

Details

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.

Examples

import greenwood as gw

y = gw.Surv(
    time=[0, 4, 7, 0, 3, 5],
    time2=[4, float("inf"), 7, 2.5, 6, 5],
    type="interval2",
)
tb = gw.Turnbull().fit(y)
tb.rmrl(s=4, tau=7, format="polars")

shape: (1, 4)

strata	s	tau	estimate
str	f64	f64	f64
"all"	4.0	7.0	1.666667

rmst()

Restricted mean survival time up to tau, under the right-endpoint convention.

Usage

Source

rmst(
    tau,
    *,
    format=None,
)

Parameters

tau: float: The end of the window. Must be positive.
format: str | None = None: Output format: None (default), "pandas", "polars", or "pyarrow". When None, a backend is auto-detected (Polars, then Pandas, then PyArrow).

Returns

pandas.DataFrame, polars.DataFrame, or pyarrow.Table: One row per curve, with columns strata, tau, and estimate. The layout is the same for one curve (strata is "all") as for several.

Details

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 reports the largest RMST consistent with the data, not an unbiased point estimate. There is no variance estimator for it, so no standard error or confidence limits are reported.

Examples

import greenwood as gw

y = gw.Surv(
    time=[0, 4, 7, 0, 3, 5],
    time2=[4, float("inf"), 7, 2.5, 6, 5],
    type="interval2",
)
tb = gw.Turnbull().fit(y)
tb.rmst(tau=7, format="polars")

shape: (1, 3)

strata	tau	estimate
str	f64	f64
"all"	7.0	4.611111

to_frame()

Return the fitted NPMLE as a DataFrame.

Usage

Source

to_frame(
    *,
    format=None,
)

Exports one row per maximal intersection interval, with its bounds, probability mass, and the survival estimate just after it.

Parameters

format: str | None = None: Output format: None (default), "pandas", "polars", or "pyarrow". When None, a backend is auto-detected (Polars, then Pandas, then PyArrow).

Returns

pandas.DataFrame, polars.DataFrame, or pyarrow.Table: A tidy table with columns strata, interval_low, interval_high, prob_mass, and estimate. The layout is the same for one curve (strata is "all") as for several.

Raises

ImportError: If the requested (or, when auto-detecting, any) DataFrame library is not installed.

Examples

import greenwood as gw

y = gw.Surv(
    time=[0, 4, 7, 0, 3, 5],
    time2=[4, float("inf"), 7, 2.5, 6, 5],
    type="interval2",
)
tb = gw.Turnbull().fit(y)
tb.to_frame(format="polars")

shape: (8, 5)

strata	interval_low	interval_high	prob_mass	estimate
str	f64	f64	f64	f64
"all"	0.0	2.5	0.333333	0.666667
"all"	2.5	3.0	4.0750e-38	0.666667
"all"	3.0	4.0	6.3494e-9	0.666667
"all"	4.0	5.0	2.9652e-26	0.666667
"all"	5.0	5.0	0.444444	0.222222
"all"	5.0	6.0	2.9652e-26	0.222222
"all"	6.0	7.0	9.1637e-74	0.222222
"all"	7.0	7.0	0.222222	1.1102e-16