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

lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
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 a trailing + 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.

lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

km = gw.KaplanMeier().fit(y)
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.

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.

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.

lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

km_loglog = gw.KaplanMeier(conf_type="log-log").fit(y)
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.

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

km_90 = gw.KaplanMeier(conf_level=0.90).fit(y)
km_99 = gw.KaplanMeier(conf_level=0.99).fit(y)

# Compare intervals at a few time points
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 KM with each transform
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)

# Compare interval widths at selected quantile times
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)

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; matches R’s survival package. 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

lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
km = gw.KaplanMeier().fit(y)

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.

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

lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
km = gw.KaplanMeier().fit(y)

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 tau based on a clinically or practically meaningful horizon, and use the same tau 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.

lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
km = gw.KaplanMeier().fit(y)

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.

lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

km_sex = gw.KaplanMeier().fit(y, by=lung["sex"])
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.

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.

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.

# Example with case weights: create sample weights for demonstration
import numpy as np

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

km_weighted = gw.KaplanMeier().fit(y, weights=weights)
km_weighted

KaplanMeier (Kaplan-Meier survival estimate)

      n  events  median  0.95LCL  0.95UCL
  281.3   204.5     320      286      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.

lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

na = gw.NelsonAalen().fit(y)
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.

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.