Estimating survival with Kaplan-Meier

The Kaplan-Meier estimator is the workhorse of survival analysis. It answers the most basic question you can ask of time-to-event data: what fraction of subjects are still event-free at each point in time? The result is the survival curve, a step function that starts at 1.0 and drops at each observed event. This page shows how to estimate it, quantify its uncertainty, summarize it with medians and restricted means, and compare subgroups.

Every example here uses the bundled lung dataset and the Surv response introduced in Survival data and the Surv object.

import greenwood as gw

# Load the bundled lung dataset as a Polars DataFrame
lung = gw.load_dataset("lung", backend="polars")

# Build a right-censored response (status == 2 marks a death)
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

# Display the response summary
y

Surv(type=right, n=228, events=165)

The response y is a Surv object holding one entry per subject: the follow-up time and whether that time ended in an event or a censoring. Displaying it shows a compact summary where / marks a censored observation. Everything on this page is estimated from this single object.

Estimating the survival function

You fit the estimator by calling fit() with a response. The KaplanMeier object follows the familiar fit-then-inspect pattern, so we fit it first and set it aside.

# Load data and build the response
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

# Fit the Kaplan-Meier estimator
km = gw.KaplanMeier().fit(y)

# Print the median survival and confidence interval
km

KaplanMeier (Kaplan-Meier survival estimate)

    n  events  median  0.95LCL  0.95UCL
  228     165     310      285      363

The fitted km object holds the whole curve internally. The most direct way to read it is to_frame(), which lays the estimate out one row per observed time. We show the first few rows.

# Export the full survival curve as a tidy one-row-per-time DataFrame
km.to_frame(format="polars")

	time f64	n_risk f64	n_event f64	n_censor f64	estimate f64	std_error f64	conf_low f64	conf_high f64
PolarsRows186Columns8
0	5	228	1	0	0.995614035088	0.00437633599855	0.987073416741	1
1	11	227	3	0	0.982456140351	0.00869464259269	0.965561897075	0.99964597882
2	12	224	1	0	0.978070175439	0.00969918321668	0.95924367691	0.997266170327
3	13	223	2	0	0.969298245614	0.0114246495327	0.947163003131	0.991950789721
4	15	221	1	0	0.964912280702	0.0121858004893	0.941321714601	0.989094052551
5	26	220	1	0	0.960526315789	0.0128955847988	0.935581072627	0.986136669838
6	30	219	1	0	0.956140350877	0.0135620698334	0.929925266889	0.983094451916
7	31	218	1	0	0.951754385965	0.0141913574455	0.92434233917	0.979979357017
183	965	3	0	1	0.0503455680708	0.0228480489161	0.020685460199	0.122534195517
184	1010	2	0	1	0.0503455680708	0.0228480489161	0.020685460199	0.122534195517
185	1022	1	0	1	0.0503455680708	0.0228480489161	0.020685460199	0.122534195517

Each row corresponds to a distinct observed time. The columns tell a complete story: n_risk is the number of subjects still under observation just before that time, n_event and n_censor count the events and censorings that occurred at it, estimate is the Kaplan-Meier survival probability, and conf_low and conf_high bound it. The estimate only steps down at times where events occur. Censorings reduce the risk set but do not move the curve.

The survival probability at any set of times is available through predict, which evaluates the step function. It returns a plain array with one probability per requested time, in the same order you asked for them.

# Evaluate the step-function survival curve at specific time points (days)
km.predict([180, 365, 730])

array([0.72167065, 0.40924162, 0.1156931 ])

Reading the survival probability

A survival probability of 0.53 at 365 days means the estimator predicts that 53 percent of subjects remain event-free one year in. Because the curve is a step function, the value holds constant between event times.

Confidence intervals

The confidence band around the curve comes from Greenwood’s variance formula. Greenwood supports three transforms for the interval, chosen with conf_type=. The default is "log", which matches the default in R’s survival package. The "log-log" transform keeps the limits inside the valid range from 0 to 1 more reliably at the tails, and "plain" gives a symmetric interval on the probability scale.

The default transform is not the same for every estimator

KaplanMeier defaults to conf_type="log", matching R. AalenJohansen (the cumulative incidence estimator for competing risks, see Competing risks) defaults to conf_type="plain" instead, matching R’s default for that estimator. If you fit both in the same analysis, set conf_type= explicitly on each rather than relying on the defaults to agree.

You set the transform when you construct the estimator, so we fit a second KaplanMeier with conf_type="log-log" and keep it separate from the default fit above.

# Load data and build the response
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

# Fit with the log-log confidence interval transform
km_loglog = gw.KaplanMeier(conf_type="log-log").fit(y)

# Print the median survival with log-log bounds
km_loglog

KaplanMeier (Kaplan-Meier survival estimate)

    n  events  median  0.95LCL  0.95UCL
  228     165     310      284      361

Pulling out just the estimate and its interval limits lets you compare this fit against the default one. We select those columns to keep the table narrow.

# Compare confidence limits from the log-log transform against the plain estimate
km_loglog.to_frame(format="polars")[["time", "estimate", "conf_low", "conf_high"]]

	time f64	estimate f64	conf_low f64	conf_high f64
PolarsRows186Columns4
0	5	0.995614035088	0.969277036219	0.999381011439
1	11	0.982456140351	0.953935230206	0.99337913279
2	12	0.978070175439	0.948119884093	0.990813247959
3	13	0.969298245614	0.936681974228	0.985244433081
4	15	0.964912280702	0.931066246164	0.982296706677
5	26	0.960526315789	0.925513676153	0.979263834158
6	30	0.956140350877	0.920018759486	0.976158015297
7	31	0.951754385965	0.914576283437	0.972988696334
183	965	0.0503455680708	0.0178661710929	0.108662176031
184	1010	0.0503455680708	0.0178661710929	0.108662176031
185	1022	0.0503455680708	0.0178661710929	0.108662176031

The point estimates are identical across transforms and only the interval limits differ. Choose the transform before fitting, and report which one you used.

Customizing the confidence level

By default, Greenwood computes 95% confidence intervals. You can adjust this with the conf_level= parameter, useful for regulatory submissions or sensitivity analyses that require different coverage levels (e.g., 90% or 99%).

# Fit with 90% and 99% confidence levels for sensitivity comparison
km_90 = gw.KaplanMeier(conf_level=0.90).fit(y)
km_99 = gw.KaplanMeier(conf_level=0.99).fit(y)

# Inspect the narrower 90% interval on the first three rows
km_90.to_frame(format="polars")[["time", "estimate", "conf_low", "conf_high"]].head(3)

	time f64	estimate f64	conf_low f64	conf_high f64
PolarsRows3Columns4
0	5	0.995614035088	0.988441563193	1
1	11	0.982456140351	0.968258314093	0.996862153069
2	12	0.978070175439	0.962245848411	0.994154736715

Note that narrower confidence levels (0.90) produce tighter bands, while wider levels (0.99) give more conservative coverage. The confidence level affects the median and quantile intervals as well.

Comparing confidence interval transforms

The three confidence interval transforms give the same point estimates. They differ only in how they construct the interval. Here is a comparison to help you choose:

import polars as pl

# Fit the same data with all three CI transforms
km_plain = gw.KaplanMeier(conf_type="plain").fit(y)
km_log = gw.KaplanMeier(conf_type="log").fit(y)
km_loglog = gw.KaplanMeier(conf_type="log-log").fit(y)

# Build a combined table showing every 30th row from each transform
results = []
for km, name in [(km_plain, "plain"), (km_log, "log"), (km_loglog, "log-log")]:
    df = km.to_frame(format="polars")
    # Select every 30th row to show interval differences across follow-up
    df_subset = df[::30].select(["time", "estimate", "conf_low", "conf_high"])
    df_subset = df_subset.with_columns(pl.lit(name).alias("transform"))
    results.append(df_subset)

# Stack all three into one frame for side-by-side comparison
pl.concat(results)

	time f64	estimate f64	conf_low f64	conf_high f64	transform str
PolarsRows21Columns5
0	5	0.995614035088	0.987036574146	1	plain
1	135	0.815253696589	0.76479399977	0.865713393409	plain
2	199	0.680272862223	0.619250223099	0.741295501347	plain
3	269	0.580686970067	0.514745920937	0.646628019198	plain
4	353	0.434044147155	0.363908912208	0.504179382101	plain
5	524	0.263190302049	0.195465597373	0.330915006725	plain
6	821	0.0671274240944	0.0210556916661	0.113199156523	plain
7	5	0.995614035088	0.987073416741	1	log
18	353	0.434044147155	0.36316544851	0.502729471467	log-log
19	524	0.263190302049	0.198161004017	0.332592324663	log-log
20	821	0.0671274240944	0.0307283845558	0.123059694935	log-log

Transform summary:

log (default): symmetric on the log scale and matches R’s survival package (this is a good all-purpose choice).
log-log: bounds intervals more reliably at the tails (best when survival is near 0 or 1).
plain: symmetric on probability scale (intuitive but can exceed [0,1] at extreme tails).

Choose the transform before fitting and report which you used.

Median survival and quantiles

# Load data and build the response
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
km = gw.KaplanMeier().fit(y)

# Return the median survival time with its confidence interval
km.median(ci=True)

(310.0, 285.0, 363.0)

The result is the median together with its confidence limits, obtained by inverting the confidence band. Other quantiles work the same way through quantile().

# Return the 25th-percentile survival time
km.quantile(0.25)

170.0

When the median is not defined

If the survival curve never falls to 0.5 within the observed follow-up, the median is not estimable and Greenwood returns nan. This is not an error, rather it is an honest statement that the data do not reach the median.

Restricted mean survival time

The mean survival time is usually not estimable from censored data, because the tail of the curve is unknown. The restricted mean survival time, or RMST, solves this by measuring the area under the survival curve up to a chosen horizon tau. It has a direct interpretation as the average event-free time over the first tau units, and it is a good summary when survival curves cross, which makes hazard ratios hard to interpret.

# Load data and build the response
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
km = gw.KaplanMeier().fit(y)

# Average event-free days in the first year, with a 95% confidence interval
km.rmst(365, ci=True)

(263.22186648200665, 247.93638355232736, 278.50734941168594)

This reports the average number of event-free days during the first year, with a confidence interval. Choose the tau value based on a clinically or practically meaningful horizon, and use the same value when comparing groups.

Restricted mean residual life

A related question is how much additional time a subject can expect, given that they have already survived to some landmark time s. The restricted mean residual life answers this: it is the area under the conditional survival curve from s to tau, divided by the survival at s. In other words, it is the restricted mean survival measured from s onward, among subjects still event-free at s. Setting s=0 recovers the restricted mean survival time.

# Load data and build the response
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
km = gw.KaplanMeier().fit(y)

# Expected additional event-free days between day 180 and 730 (s and tau values), for
# survivors at day 180
km.rmrl(180, 730, ci=True)

(275.7027711565545, 242.6533843723124, 308.75215794079656)

Here the value is the expected number of additional event-free days between 180 and 730 days, for a subject known to be alive at 180 days. This is useful for updating a prognosis partway through follow-up, and it complements the conditional survival curves shown in Cox model diagnostics.

Comparing subgroups

Passing a grouping variable to by fits a separate curve within each stratum, which is the starting point for any subgroup comparison. The tidy output gains a strata column.

# Load data and build the response
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

# Fit a separate Kaplan-Meier curve for each level of sex
km_sex = gw.KaplanMeier().fit(y, by=lung["sex"])

# Print per-stratum median survival
km_sex

KaplanMeier (Kaplan-Meier survival estimate)

     n  events  median  0.95LCL  0.95UCL
1  138     112     270      212      310
2   90      53     426      348      550

The grouped fit km_sex holds one curve per level of sex. Its tidy output gains a strata column identifying which curve each row belongs to. We show the first couple of rows within each stratum so you can see both curves side by side without printing every time point.

# Show the first two rows from each stratum to see both curves side by side
km_sex.to_frame(format="polars").group_by("strata").head(2)

	strata obj	time f64	n_risk f64	n_event f64	n_censor f64	estimate f64	std_error f64	conf_low f64	conf_high f64
PolarsRows4Columns9
0	2	5	90	1	0	0.988888888889	0.011049210289	0.967468240209	1
1	2	60	89	1	0	0.977777777778	0.0155379088618	0.947793404603	1
2	1	11	138	3	0	0.978260869565	0.0124139182765	0.954230116249	1
3	1	12	135	1	0	0.971014492754	0.014281169221	0.943423496588	0.999412404448

Per-stratum summaries follow naturally. median and rmst return a dictionary keyed by stratum when the fit is grouped.

# Per-stratum median survival times as a dictionary keyed by stratum label
km_sex.median()

{1: 270.0, 2: 426.0}

To test whether the difference between groups is statistically significant, rather than just describing it, use the log-rank family covered in Comparing groups.

Handling case weights

When your data contains frequency weights (e.g., summarized counts rather than individual records), or case weights for design-based analysis, pass them via the weights parameter.

import numpy as np

# Create sample weights for demonstration
weights = np.random.uniform(0.5, 2.0, len(y))

# Fit with frequency/case weights applied to each observation
km_weighted = gw.KaplanMeier().fit(y, weights=weights)

# Display the weighted median survival
km_weighted

KaplanMeier (Kaplan-Meier survival estimate)

      n  events  median  0.95LCL  0.95UCL
  288.1   207.9     329      288      364

The estimates adjust for the weights using the weighted Kaplan-Meier formula. Confidence intervals and other summaries (median, RMST) also account for the weights automatically.

Robust (sandwich) variance for weighted curves

Greenwood’s variance formula treats weighted counts as though they were ordinary integer counts. When weights are non-integer, for example inverse-probability-of-censoring weights or survey design weights, this understates the true uncertainty. The robust (infinitesimal-jackknife) variance corrects for this by computing each subject’s influence on the survival estimate and summing the squared influences.

import numpy as np

np.random.seed(23)
weights = np.random.uniform(0.5, 2.0, len(y))

# Greenwood SE (default): treats weighted counts as ordinary counts
km_green = gw.KaplanMeier().fit(y, weights=weights)

# Robust SE: infinitesimal-jackknife, correct for non-integer weights
km_robust = gw.KaplanMeier(robust=True).fit(y, weights=weights)

The survival estimates are identical. Only the standard errors and confidence intervals change.

# Compare standard errors at the median event time
idx = len(km_green.time_) // 2
t_mid = km_green.time_[idx]
print(f"Time {t_mid}: Greenwood SE = {km_green.std_error_[idx]:.6f}, "
      f"Robust SE = {km_robust.std_error_[idx]:.6f}")

Time 276.0: Greenwood SE = 0.030327, Robust SE = 0.035603

When to use robust=True

Use robust=True whenever weights are non-integer (IPW, survey weights). For unweighted data or integer frequency weights, Greenwood’s formula is efficient and correct.

Clustered robust variance

When subjects are nested within clusters (e.g., patients within hospitals, or recurring events on the same individual), observations within a cluster are correlated. The standard robust variance ignores this correlation. Passing cluster= to fit() accounts for it by summing per-subject influences within each cluster before squaring.

# Load data and build the response
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

# Cluster by institution (inst column) (NaN rows are dropped automatically)
km_cluster = gw.KaplanMeier().fit(y, cluster=lung["inst"])

# Compare clustered SE to unclustered robust SE at the midpoint
km_robust = gw.KaplanMeier(robust=True).fit(y)
idx = len(km_cluster.time_) // 2
t_mid = km_cluster.time_[idx]
print(
    f"Time {t_mid}: Robust SE = {km_robust.std_error_[idx]:.6f}, "
    f"Clustered SE = {km_cluster.std_error_[idx]:.6f}"
)

Time 276.0: Robust SE = 0.033759, Clustered SE = 0.029427

Passing cluster= implies robust=True. The survival estimates are unchanged and only the standard errors and confidence intervals differ.

When to use cluster=

Use cluster= whenever your data has a natural grouping that induces within-group correlation: multi-center trials (cluster by site), recurrent-event data (cluster by subject), or family studies (cluster by family). Without clustering, the robust SE understates uncertainty by treating correlated observations as independent.

Bootstrap confidence intervals

The analytical confidence intervals above (Greenwood, robust, clustered) are fast and well-calibrated for single quantities like the survival curve. However, some summary statistics lack simple closed-form intervals. The median survival difference between two groups is a classic example: each group’s median has an analytical CI, but the difference does not. The bootstrap() function fills this gap by resampling subjects with replacement.

# Load data and build the response
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

# Bootstrap CI for the median survival time
gw.bootstrap(y, "median", n_boot=1000, seed=23)

BootstrapResult (n_boot=1000, ci_type='percentile')

  estimate      se  0.95LCL  0.95UCL
   310.000  23.120  284.000  361.000

For two-group comparisons, pass by= with a _diff statistic. The function resamples within each group and computes the difference on every replicate.

# Bootstrap CI for the difference in median survival between sexes
gw.bootstrap(y, "median_diff", by=lung["sex"], n_boot=1000, seed=23)

BootstrapResult (n_boot=1000, ci_type='percentile')

  estimate      se   0.95LCL  0.95UCL
  -156.000  62.275  -285.000  -59.975

The RMST difference works the same way. This is particularly useful when survival curves cross, making hazard ratios hard to interpret.

# Bootstrap CI for the difference in 1-year RMST between sexes
gw.bootstrap(y, "rmst_diff", by=lung["sex"], tau=365.0, n_boot=1000, seed=23)

BootstrapResult (n_boot=1000, ci_type='percentile')

  estimate      se  0.95LCL  0.95UCL
   -55.970  14.943  -84.704  -28.154

For custom statistics, pass a callable that takes a fitted KaplanMeier and returns a float.

# Bootstrap CI for a custom statistic: area under the curve from 180 to 365 days
gw.bootstrap(y, lambda km: km.rmst(365) - km.rmst(180), n_boot=500, seed=23)

BootstrapResult (n_boot=500, ci_type='percentile')

  estimate     se  0.95LCL  0.95UCL
   105.115  5.553   94.331  116.255

Three interval types are available via ci_type=: "percentile" (default), "normal" (symmetric around the estimate), and "bca" (bias-corrected and accelerated, more accurate for skewed distributions but slower).

Choosing the number of replicates

For percentile intervals, 1000 replicates is usually sufficient. For BCa intervals, 2000 or more may be needed for stability. The seed= parameter ensures reproducibility.

Nelson-Aalen cumulative hazard

Closely related to survival is the cumulative hazard, the accumulated risk of the event over time. The Nelson-Aalen estimator computes it directly as a running sum of events over subjects at risk: H(t) = \sum_{t_i \le t} d_i / n_i. It is often preferred for diagnostics because the cumulative hazard is roughly linear when the underlying hazard is constant, making departures easier to spot.

The estimator follows the same fit-then-inspect pattern as KaplanMeier.

# Fit the Nelson-Aalen cumulative hazard estimator
na = gw.NelsonAalen().fit(y)
na

NelsonAalen (Nelson-Aalen cumulative hazard estimate)

    n  events  max cumhaz
  228     165       2.889

The summary shows the number of subjects, events, and the maximum cumulative hazard reached. You can convert to a survival estimate via S(t) = \exp(-H(t)), though KaplanMeier is typically preferred for direct survival estimation.

Exporting the cumulative hazard table

Reading the fit with to_frame() gives one row per event time. The estimate column holds the cumulative hazard rather than a survival probability. Unlike survival, it starts near 0 and climbs as events accumulate.

# Export the cumulative hazard curve as a tidy DataFrame
na.to_frame(format="polars")

	time f64	n_risk f64	n_event f64	estimate f64	std_error f64	conf_low f64	conf_high f64
PolarsRows186Columns7
0	5	228	1	0.00438596491228	0.00438596491228	0.000617822342514	0.031136278001
1	11	227	3	0.0176018239431	0.00880092787832	0.00660626714774	0.0468985281999
2	12	224	1	0.0220661096574	0.00986844356817	0.00918438227555	0.0530153450504
3	13	223	2	0.0310347195229	0.0117304799526	0.0147948751781	0.0651005029965
4	15	221	1	0.0355596064007	0.012572937651	0.0177825713872	0.0711081417777
5	26	220	1	0.0401050609462	0.0133693649138	0.020866224206	0.0770822692987
6	30	219	1	0.0446712709918	0.0141276393067	0.0240341714416	0.0830285519462
7	31	218	1	0.0492584269551	0.0148536928813	0.0272774597424	0.0889522942756
183	965	3	0	2.88926746252	0.418689933124	2.17489507081	3.83828469797
184	1010	2	0	2.88926746252	0.418689933124	2.17489507081	3.83828469797
185	1022	1	0	2.88926746252	0.418689933124	2.17489507081	3.83828469797

Each row includes the risk set size (n_risk), the number of events (n_event), the cumulative hazard estimate, its standard error, and confidence limits. The confidence interval uses a log transform by default, which provides better coverage in the tails. Pass conf_type="plain" to NelsonAalen() for Wald-type intervals instead.

Stratified cumulative hazard

Pass by= to fit separate cumulative hazard curves per group, just as with KaplanMeier.

# Stratified cumulative hazard by sex
na_strat = gw.NelsonAalen().fit(y, by=lung["sex"])
na_strat

NelsonAalen (Nelson-Aalen cumulative hazard estimate)

     n  events  max cumhaz
1  138     112       3.163
2   90      53       2.322

The stratified fit reports one curve per group. The to_frame() output includes a strata column so you can filter or plot each group separately.

na_strat.to_frame(format="polars")

	strata obj	time f64	n_risk f64	n_event f64	estimate f64	std_error f64	conf_low f64	conf_high f64
PolarsRows206Columns8
0	1	11	138	3	0.0217391304348	0.0125510928085	0.0070113351632	0.0674036800495
1	1	12	135	1	0.0291465378422	0.0145739361597	0.0109387253179	0.0776617607166
2	1	13	134	2	0.0440719109765	0.0179939711233	0.0197982903792	0.0981061142107
3	1	15	132	1	0.0516476685523	0.0195237060937	0.0246196522029	0.108347658403
4	1	26	131	1	0.0592812563385	0.020962985525	0.0296425977615	0.118554634832
5	1	30	130	1	0.0669735640308	0.0223297639879	0.0348419883768	0.128737149858
6	1	31	129	1	0.0747255020153	0.0236370662789	0.0401992508602	0.13890558983
7	1	53	128	2	0.0903505020153	0.0260917844306	0.051299972923	0.159127047234
203	2	765	3	1	2.32191477538	0.527536959731	1.48748997011	3.62441988347
204	2	821	2	0	2.32191477538	0.527536959731	1.48748997011	3.62441988347
205	2	965	1	0	2.32191477538	0.527536959731	1.48748997011	3.62441988347

Summarizing with tidy and glance

Like KaplanMeier, the Nelson-Aalen estimator works with tidy() and glance(). Calling tidy() returns the same table as to_frame(), while glance() produces a one-row-per-stratum summary with the starting risk set, total events, and the maximum cumulative hazard.

# One-row summary of the unstratified fit
gw.glance(na, format="polars")

	n_start f64	events f64	max_cumhaz f64
PolarsRows1Columns3
0	228	165	2.88926746252

For stratified fits, glance() returns one row per group.

# Per-group summary
gw.glance(na_strat, format="polars")

	strata i64	n_start f64	events f64	max_cumhaz f64
PolarsRows2Columns4
0	1	138	112	3.1628387723
1	2	90	53	2.32191477538

Turnbull’s algorithm for interval-censored data

Both KaplanMeier and NelsonAalen assume you know each subject’s event time exactly, up to right censoring. That assumption breaks down when follow-up is periodic rather than continuous: a clinic visit, an equipment inspection, or a survey wave tells you only that the event happened sometime between two observations, not exactly when. Turnbull handles this directly, and can mix exact, left-, right-, and interval-censored observations in the same fit. The example below has six subjects: two exact deaths, one right-censored subject with no event observed by the last visit, and three whose event is only known to fall somewhere within a window.

# Build an interval-censored response
y_interval = 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_interval)
tb

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

  n  atoms  ambiguous  iters  converged
  6      8          5     32       True

Unlike Kaplan-Meier, the estimator’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 where within it the mass truly sits. to_frame() reports both bounds of every such interval:

# Export one row per maximal intersection interval
tb.to_frame(format="polars")

	interval_low f64	interval_high f64	prob_mass f64	estimate f64
PolarsRows8Columns4
0	0	2.5	0.230327668542	0.769672331458
1	2.5	3	2.23126768588e-19	0.769672331458
2	3	4	8.47883641837e-10	0.76967233061
3	4	4	0.372677995402	0.396994335208
4	4	5	1.72469616433e-19	0.396994335208
5	5	5	0.230327668542	0.166666666667
6	5	6	1.72469616433e-19	0.166666666667
7	7	7	0.166666666667	1.11022302463e-16

Where an interval degenerates to a single point (interval_low equals interval_high), that’s an exact death or a region every subject’s constraint resolves unambiguously, and survival is known exactly there. Elsewhere it genuinely is not identified by the data, so predict() and quantile() return nan rather than silently interpolate a point that the data cannot support. Here, t=1 falls strictly inside the ambiguous window between 0 and 2.5, so it comes back nan:

# Evaluate the survival curve at specific times
tb.predict([1, 2.5, 5, 7])

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

The quantile() method follows the same rule, but always returns a 3-tuple (estimate, lower, upper) regardless of whether the crossing is ambiguous: estimate is nan exactly when it falls inside an unresolved interval, in which case lower and upper bracket it rather than a confidence bound.

# First-quartile survival time, or its ambiguity bracket
tb.quantile(0.25)

(4.0, 4.0, 4.0)

When there is no genuine interval ambiguity in the data (only exact events and right censoring), Turnbull reduces exactly to KaplanMeier, since that is the one case where the NPMLE is fully identified everywhere. The comparison below reuses the right-censored response y from earlier in this page and checks agreement at KaplanMeier’s own event times. The raw arrays aren’t directly comparable, since Turnbull’s atom table also carries a few near-zero-mass artifacts from the general construction.

tb_right = gw.Turnbull().fit(y)
km_check = gw.KaplanMeier().fit(y)
import numpy as np

# Compare survival at KM's own event times
np.allclose(tb_right.predict(km_check.time_), km_check.survival_, atol=1e-8)

True

That equivalence is a useful check on your own data: if predict() and quantile() never return nan for a fit, Turnbull and KaplanMeier will agree everywhere, and you could have used KaplanMeier directly. The three methods above cover the whole contract: to_frame() for the fitted atoms, predict() to evaluate the curve, and quantile() to invert it, all while keeping identifiability limits visible rather than papering over them.

Restricted mean survival and residual life

rmst() and rmrl() need a single number, so unlike predict() and quantile() they do not return nan for ambiguous regions. Instead they fall back to Turnbull’s own convention for reporting a plottable curve: every atom’s probability mass, ambiguous or not, is treated as resolving at that atom’s right endpoint. This is the most conservative choice available: placing mass as late as possible maximizes the area under the curve, so rmst()/rmrl() report the largest value consistent with the data, not an unbiased point estimate. There is no variance estimator for either, so no confidence interval is available.

# Restricted mean survival time up to t=7
tb.rmst(7)

4.384836165729157

Use rmst()/rmrl() when you need one reportable number, but check to_frame() first to see how much of the mass over your chosen horizon actually sits in non-degenerate intervals: the more ambiguity there is, the further the right-endpoint value can drift from an unbiased estimate.

Stratified fits

Pass by= to fit separate curves per group, just as with KaplanMeier and NelsonAalen.

# Fit separate curves for two groups
tb_strat = gw.Turnbull().fit(y_interval, by=[0, 1, 0, 1, 0, 1])
tb_strat.median()

{0: (nan, 3.0, 4.0), 1: (4.0, 4.0, 4.0)}

Turnbull’s algorithm is a specialized tool: use it only when your data are genuinely interval-censored, since it costs more to fit than Kaplan-Meier and its answer is honestly partial where the data are. The alternative, rounding each observation to an interval endpoint and calling KaplanMeier, is faster but silently biased. When your follow-up is truly periodic, the nans and ambiguity brackets you see above are the estimator doing its job.

Next steps

You can now estimate survival curves, summarize them, and split them by group.

Comparing groups tests whether survival differs between strata with the log-rank family.
Visualizing survival turns these curves into publication-quality figures with plotnine, including numbers-at-risk tables.
Cox regression moves from describing groups to modeling the effect of continuous and multiple covariates.