# 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](../reference/Surv.md#greenwood.Surv) response introduced in [Survival data and the Surv object](survival-data.md).
``` python
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](../reference/Surv.md#greenwood.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](../reference/AFT.md#greenwood.AFT.fit) with a response. The [KaplanMeier](../reference/KaplanMeier.md#greenwood.KaplanMeier) object follows the familiar fit-then-inspect pattern, so we fit it first and set it aside.
``` python
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()](../reference/AFT.md#greenwood.AFT.to_frame), which lays the estimate out one row per observed time. We show the first few rows.
``` python
km.to_frame(format="polars")
```
PolarsRows186Columns8
|
|
time
f64
|
n_risk
f64
|
n_event
f64
|
n_censor
f64
|
estimate
f64
|
std_error
f64
|
conf_low
f64
|
conf_high
f64
|
| 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](../reference/AFT.md#greenwood.AFT.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.
``` python
km.predict([180, 365, 730])
```
array([0.72167065, 0.40924162, 0.1156931 ])
> **Note: 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](../reference/KaplanMeier.md#greenwood.KaplanMeier) with `conf_type="log-log"` and keep it separate from the default fit above.
``` python
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.
``` python
km_loglog.to_frame(format="polars")[["time", "estimate", "conf_low", "conf_high"]]
```
PolarsRows186Columns4
|
|
time
f64
|
estimate
f64
|
conf_low
f64
|
conf_high
f64
|
| 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](../reference/RMSTResult.md#greenwood.RMSTResult.conf_level) parameter, useful for regulatory submissions or sensitivity analyses that require different coverage levels (e.g., 90% or 99%).
``` python
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)
```
PolarsRows3Columns4
|
|
time
f64
|
estimate
f64
|
conf_low
f64
|
conf_high
f64
|
| 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:
``` python
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)
```
PolarsRows21Columns5
|
|
time
f64
|
estimate
f64
|
conf_low
f64
|
conf_high
f64
|
transform
str
|
| 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
``` python
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](../reference/KaplanMeier.md#greenwood.KaplanMeier.quantile).
``` python
km.quantile(0.25)
```
170.0
> **Tip: 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.
``` python
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.
``` python
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](cox-diagnostics.md).
# 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.
``` python
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.
``` python
km_sex.to_frame(format="polars").group_by("strata").head(2)
```
PolarsRows4Columns9
|
|
strata
obj
|
time
f64
|
n_risk
f64
|
n_event
f64
|
n_censor
f64
|
estimate
f64
|
std_error
f64
|
conf_low
f64
|
conf_high
f64
|
| 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](../reference/KaplanMeier.md#greenwood.KaplanMeier.median) and [rmst](../reference/KaplanMeier.md#greenwood.KaplanMeier.rmst) return a dictionary keyed by stratum when the fit is grouped.
``` python
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](comparing-groups.md).
## 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.
``` python
# 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](../reference/KaplanMeier.md#greenwood.KaplanMeier), so we fit it in its own step.
``` python
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()](../reference/AFT.md#greenwood.AFT.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.
``` python
na.to_frame(format="polars")
```
PolarsRows186Columns7
|
|
time
f64
|
n_risk
f64
|
n_event
f64
|
estimate
f64
|
std_error
f64
|
conf_low
f64
|
conf_high
f64
|
| 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](comparing-groups.md) tests whether survival differs between strata with the log-rank family.
- [Visualizing survival](visualization.md) turns these curves into publication-quality figures with plotnine, including numbers-at-risk tables.
- [Cox regression](cox-regression.md) moves from describing groups to modeling the effect of continuous and multiple covariates.