# Turnbull


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


Usage

``` python
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](Turnbull.md#greenwood.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_](Turnbull.md#greenwood.Turnbull.interval_high_)); see their docstrings for the resulting conservative (upper-bound) bias.

To use this estimator, call `fit()` with a [Surv](Surv.md#greenwood.Surv) response built via [Surv.interval()](Surv.md#greenwood.Surv.interval) (the general case), or [Surv.right()](Surv.md#greenwood.Surv.right)/[Surv.left()](Surv.md#greenwood.Surv.left) (degenerate cases: fitting a right-censored response through [Turnbull](Turnbull.md#greenwood.Turnbull) reproduces [KaplanMeier](KaplanMeier.md#greenwood.KaplanMeier) exactly, since there is then no genuine interval ambiguity). Left-truncated ([Surv.counting()](Surv.md#greenwood.Surv.counting)) 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_](Turnbull.md#greenwood.Turnbull.interval_low_), [interval_high_](Turnbull.md#greenwood.Turnbull.interval_high_), [prob_mass_](Turnbull.md#greenwood.Turnbull.prob_mass_), `survival_`), accessible as aligned arrays or exported to DataFrames.


## Details

Call `fit` with a [Surv](Surv.md#greenwood.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:


``` python
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)
tb
```


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

      n  atoms  ambiguous  iters  converged
      6      8          5     32       True


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


``` python
tb.to_frame(format="polars")
```


shape: (8, 4)

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


## Attributes

| Name | Description |
|----|----|
| [interval_high_](#interval_high_) | Upper bound of each maximal intersection interval. |
| [interval_low_](#interval_low_) | Lower bound of each maximal intersection interval. |
| [prob_mass_](#prob_mass_) | Probability mass assigned to each maximal intersection interval. |
| [strata_](#strata_) | Stratum labels for each row, or `None` for unstratified fits. |
| [survival_](#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\_


Stratum labels for each row, or `None` for unstratified fits.


`strata_: Array | None`


------------------------------------------------------------------------


### survival\_


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


`survival_: Array`


## Methods

| Name | Description |
|----|----|
| [fit()](#fit) | Fit Turnbull's NPMLE to interval-censored survival data. |
| [median()](#median) | Median survival time per stratum (the 0.5-quantile). |
| [predict()](#predict) | Evaluate the survival curve at specified times. |
| [quantile()](#quantile) | Return the `p`-quantile survival time per stratum. |
| [rmrl()](#rmrl) | Restricted mean residual life at `s`, under the right-endpoint convention. |
| [rmst()](#rmst) | Restricted mean survival time up to `tau`, under the right-endpoint convention. |
| [to_frame()](#to_frame) | Return the fitted NPMLE as a DataFrame. |

------------------------------------------------------------------------


### fit()


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


Usage

``` python
fit(
    surv,
    *,
    by=None,
    weights=None,
)
```


Computes the maximal intersection intervals and their probability masses from a [Surv](Surv.md#greenwood.Surv) response, via EM self-consistency. Pass `by=` to fit separate curves per group (stratified analysis).


#### Parameters


`surv: Surv`  
A [Surv](Surv.md#greenwood.Surv) response built with [Surv.interval()](Surv.md#greenwood.Surv.interval), [Surv.right()](Surv.md#greenwood.Surv.right), or [Surv.left()](Surv.md#greenwood.Surv.left). Left-truncated ([Surv.counting()](Surv.md#greenwood.Surv.counting)) and multi-state responses raise `NotImplementedError`.

`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`): uses `surv.weights` if present, otherwise unit weights.


#### Returns


`Turnbull`  
The fitted estimator itself (for method chaining), with cached results ([interval_low_](Turnbull.md#greenwood.Turnbull.interval_low_), [interval_high_](Turnbull.md#greenwood.Turnbull.interval_high_), [prob_mass_](Turnbull.md#greenwood.Turnbull.prob_mass_), `survival_`, …) as attributes.


#### Examples


``` python
import 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.survival_
```


    array([7.69672331e-01, 7.69672331e-01, 7.69672331e-01, 3.96994335e-01,
           3.96994335e-01, 1.66666667e-01, 1.66666667e-01, 1.11022302e-16])


------------------------------------------------------------------------


### median()


Median survival time per stratum (the 0.5-quantile).


Usage

``` python
median()
```


A convenience wrapper around `quantile(0.5)`. See `quantile` for the return shape.


#### Returns


`tuple or dict`  
`(estimate, lower, upper)` for a single stratum, or a `dict` keyed by stratum label for stratified fits.


#### Examples


``` python
import 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.median()
```


    (4.0, 4.0, 4.0)


------------------------------------------------------------------------


### predict()


Evaluate the survival curve at specified times.


Usage

``` python
predict(times)
```


#### Parameters


`times: Any`  
Query times at which to evaluate the curve. Can be a scalar or array-like of floats.


#### Returns


`ndarray or dict`  
For a single stratum: an array (matching `times`' 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.


#### Details

Outside every ambiguous interval, the curve is a well-defined, right-continuous step function, exactly as for [KaplanMeier](KaplanMeier.md#greenwood.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


``` python
import 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)

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


    array([           nan, 7.69672331e-01, 1.66666667e-01, 1.11022302e-16])


------------------------------------------------------------------------


### quantile()


Return the `p`-quantile survival time per stratum.


Usage

``` python
quantile(p)
```


#### Parameters


`p: float`  
Quantile level between 0 and 1 (e.g. `p=0.5` for the median).


#### Returns


`tuple or dict`  
For a single stratum: `(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.


#### Details

The quantile is found by inverting the step-function survival curve, exactly as for [KaplanMeier](KaplanMeier.md#greenwood.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.


#### Examples


``` python
import 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)

# The first-quartile survival time (or its ambiguity bracket)
tb.quantile(0.25)
```


    (4.0, 4.0, 4.0)


------------------------------------------------------------------------


### rmrl()


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


Usage

``` python
rmrl(
    s,
    tau,
)
```


#### Parameters


`s: float`  
The landmark time. Must be non-negative.

`tau: float`  
The upper time limit for the restriction. Must be greater than `s`.


#### Returns


`float or dict`  
The restricted mean residual life for a single stratum, or a `dict` keyed by stratum label for stratified fits. `nan` if everyone has resolved (under the same right-endpoint convention) by time `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_](Turnbull.md#greenwood.Turnbull.interval_high_)); see `rmst` for why, and for the resulting conservative (upper-bound) bias.


#### Examples


``` python
import 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.rmrl(4, 7)
```


    1.839642543409065


------------------------------------------------------------------------


### rmst()


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


Usage

``` python
rmst(tau)
```


#### Parameters


`tau: float`  
The upper time limit for the restriction. Must be positive.


#### Returns


`float or dict`  
The restricted mean survival time for a single stratum, or a `dict` keyed by stratum label for stratified fits.


#### 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_](Turnbull.md#greenwood.Turnbull.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).


#### Examples


``` python
import 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.rmst(7)
```


    4.384836165729157


------------------------------------------------------------------------


### to_frame()


Return the fitted NPMLE as a DataFrame.


Usage

``` python
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 `interval_low`, `interval_high`, `prob_mass`, `estimate`, and optionally `strata`.


#### Raises


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


#### Examples


``` python
import 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")
```


shape: (8, 4)

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