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.

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; 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	1	11	138	3	0	0.978260869565	0.0124139182765	0.954230116249	1
1	1	12	135	1	0	0.971014492754	0.014281169221	0.943423496588	0.999412404448
2	2	5	90	1	0	0.988888888889	0.011049210289	0.967468240209	1
3	2	60	89	1	0	0.977777777778	0.0155379088618	0.947793404603	1

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
  280.8   202.7     310      285      353

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.

Nelson-Aalen cumulative hazard

Closely related to survival is the cumulative hazard, the accumulated risk of the event over time. The Nelson-Aalen estimator provides it directly, and it is often preferred for diagnostics because it is roughly linear when the hazard is constant.

The estimator follows the same fit-then-inspect pattern as KaplanMeier, so we fit it in its own step.

# 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 Nelson-Aalen cumulative hazard estimator
na = gw.NelsonAalen().fit(y)

# Print a summary of the cumulative hazard
na

NelsonAalen (Nelson-Aalen cumulative hazard estimate)

    n  events  max cumhaz
  228     165       2.889

Reading na with to_frame() gives the same table layout as before, but now the estimate column holds the cumulative hazard rather than the survival probability. Unlike survival, it starts at 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

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.