Kaplan-Meier survival curves and the log-rank test

A renal unit in Kuala Lumpur follows its dialysis patients and asks how long they survive once dialysis starts, and whether diabetic patients fare worse. Patients enter on different dates and the study closes on one fixed day, so when you freeze the data some have died and others are still alive. The outcome is not a plain yes or no. It is a yes-or-no paired with a time: how long until the event, and did the event happen at all. That pairing is time-to-event data, and it needs its own tools.

You met rates and person-time in Part D. Survival analysis uses the same raw material, the time each person is followed, but keeps the whole follow-up curve instead of collapsing it to a single rate. This lesson builds the Kaplan-Meier estimate of the survival curve, reads a median off it, and tests whether two curves differ with the log-rank test.

What makes time-to-event data different

A binary outcome asks one question: did the event happen? Time-to-event asks two: did it happen, and when? A patient who dies at 8 months and one who dies at 80 months both count as a death, yet they tell different stories. Throwing away the timing and reporting only the proportion dead wastes most of what the data hold.

The second feature is harder: you rarely observe everyone's event. The study ends, a patient transfers out, or someone is still alive on the closing date. For those people you know the event has not happened yet, but not when it will. A plain proportion cannot handle that, because you cannot label them dead or fully survived.

Censoring

When the study closes, a patient who is still alive has a true survival time you cannot see. You know it is longer than their follow-up so far, but not by how much. This is right-censoring: the event lies somewhere to the right of the last time you observed the person. A patient followed for 60 months and still alive is censored at 60. A patient who transfers out at 28 months, alive when last seen, is censored at 28.

Censoring is information, not missing data, and how you treat it decides whether your answer is honest. Two tempting shortcuts both bias the result:

1. Dropping censored patients. If you delete everyone still alive and analyse only the deaths, you keep the patients who did worst and discard the survivors. The estimated survival collapses far below the truth.
2. Counting censored patients as having survived the whole study. Treating a transfer at 28 months as if they were event-free at 60 months invents follow-up you never had, and pushes the estimate too high.

The Kaplan-Meier method threads between these errors. It uses every patient for exactly as long as you watched them, and no longer.

The survival function

The object we want to estimate is the survival function , the probability that a patient survives past time .

At the start, before anyone can have the event, : everyone is alive. As time passes the curve can only fall or stay flat, never rise, because once a patient has the event they cannot un-have it. Read as: an estimated 70% of patients are still alive two years after dialysis starts.

The Kaplan-Meier estimator

The Kaplan-Meier estimate turns the follow-up times into a step function. The idea is to chain together the chance of getting through each moment when a death actually occurs. Build it in steps:

1. List the distinct times at which at least one death happens, in order. Nothing happens to the curve between these times.
2. At each death time , count the risk set , the patients still alive and still being followed just before , and the number of deaths at that time.
3. The chance of surviving that instant is .
4. Multiply these conditional survival chances together up to time .

This is the product-limit formula. Censoring enters through the risk set . A censored patient stays in the denominator up to the moment they leave, then drops out of the next risk set. Their departure does not pull the curve down, because no death occurred. The result is a curve that is flat wherever there is no death, with a vertical drop at every death time whose size depends on how many were still at risk.

Reading the curve and the median survival

The curve starts at 1 and steps down. To read the median survival time, find where it first crosses 0.5 and drop straight to the time axis: that is the point by which half the patients have had the event. The median is the standard summary because, unlike a mean, it can be estimated before every patient has died. If the curve never reaches 0.5, the median is not reached within follow-up, and you report that rather than guessing past your data.

Worked example 1: fit and plot a Kaplan-Meier curve

Here is the dialysis cohort: 20 patients, 10 with diabetes and 10 without, with survival time in months and a status flag (1 = died, 0 = censored, still alive or transferred out). We fit the Kaplan-Meier estimate for both groups and plot the two step functions.

webr::install("survival")
library(survival)
time   <- c(8, 12, 14, 16, 18, 20, 24, 26, 30, 36,
            20, 28, 32, 36, 40, 44, 48, 52, 56, 60)
status <- c(1, 1, 0, 1, 1, 1, 0, 1, 1, 1,
            1, 0, 1, 0, 1, 0, 1, 0, 1, 0)
group  <- c(rep("diabetic", 10), rep("non-diabetic", 10))
dat    <- data.frame(time = time, status = status, group = group)

km <- # fit survfit(Surv(time, status) ~ group)

km <- survfit(Surv(time, status) ~ group, data = dat)
plot(km, col = c("firebrick", "steelblue"), lwd = 2,
     xlab = "Months since dialysis start",
     ylab = "Survival probability S(t)")
legend("topright", legend = c("diabetic", "non-diabetic"),
       col = c("firebrick", "steelblue"), lwd = 2)
km
# group=diabetic       median 20
# group=non-diabetic   median 48

import numpy as np
import matplotlib.pyplot as plt
from statsmodels.duration.survfunc import SurvfuncRight
time = np.array([8, 12, 14, 16, 18, 20, 24, 26, 30, 36,
                 20, 28, 32, 36, 40, 44, 48, 52, 56, 60])
status = np.array([1, 1, 0, 1, 1, 1, 0, 1, 1, 1,
                   1, 0, 1, 0, 1, 0, 1, 0, 1, 0])
group = np.array(["diabetic"] * 10 + ["non-diabetic"] * 10)

sf = # SurvfuncRight(time, status) for one group

for g, col in [("diabetic", "firebrick"), ("non-diabetic", "steelblue")]:
    m = group == g
    sf = SurvfuncRight(time[m], status[m])
    plt.step(sf.surv_times, sf.surv_prob, where="post", color=col, label=g)
    below = sf.surv_times[sf.surv_prob <= 0.5]
    median = below[0] if len(below) else np.nan
    print(g, "median:", median)
plt.xlabel("Months since dialysis start")
plt.ylabel("Survival probability S(t)")
plt.legend()
plt.show()
# diabetic median: 20
# non-diabetic median: 48

Both curves start at 1 and step down. The diabetic curve falls faster and crosses 0.5 at about 20 months; the non-diabetic curve holds higher and reaches its median near 48 months, roughly twice as long. The flat stretches in the non-diabetic curve are where censored patients left without a death, exactly the censoring handling described above.

Comparing two survival curves: the log-rank test

Two curves separating on a plot is suggestive, not conclusive. The separation could be sampling noise, especially with 10 patients a group. The log-rank test asks whether the gap is more than chance would produce. Its null hypothesis is that the two groups share the same survival function at every time, so neither group has a higher risk of the event.

The test works through every death time. At each one it asks: given each group's risk set then, how many deaths would you expect in the diabetic group if both groups had equal risk? It sums observed and expected deaths in one group across all event times and compares them. A large gap between total observed and expected, relative to its variance, gives a large chi-square on 1 degree of freedom and a small p-value.

The log-rank test assumes little. It does not assume any shape for the curves, nor a constant event rate over time. Its one real requirement is proportional hazards: the ratio of the two groups' hazards stays roughly constant over follow-up. It has most power when the curves separate steadily, and loses power when they cross, because early and late differences cancel inside the sum. It returns a p-value, not an effect size; for how much worse one group fares, pair it with a hazard ratio from a Cox model, which you meet next.

Worked example 2: the log-rank test

Use the same 20-patient cohort and test whether survival differs between the diabetic and non-diabetic groups.

webr::install("survival")
library(survival)
time   <- c(8, 12, 14, 16, 18, 20, 24, 26, 30, 36,
            20, 28, 32, 36, 40, 44, 48, 52, 56, 60)
status <- c(1, 1, 0, 1, 1, 1, 0, 1, 1, 1,
            1, 0, 1, 0, 1, 0, 1, 0, 1, 0)
group  <- c(rep("diabetic", 10), rep("non-diabetic", 10))
dat    <- data.frame(time = time, status = status, group = group)

lr <- # survdiff(Surv(time, status) ~ group)

lr <- survdiff(Surv(time, status) ~ group, data = dat)
lr
# N Observed Expected
# diabetic      10    8     3.11
# non-diabetic  10    5     9.89
# Chisq = 13.7  on 1 degrees of freedom, p = 2e-04

import numpy as np
from statsmodels.duration.survfunc import survdiff
time = np.array([8, 12, 14, 16, 18, 20, 24, 26, 30, 36,
                 20, 28, 32, 36, 40, 44, 48, 52, 56, 60])
status = np.array([1, 1, 0, 1, 1, 1, 0, 1, 1, 1,
                   1, 0, 1, 0, 1, 0, 1, 0, 1, 0])
group = np.array(["diabetic"] * 10 + ["non-diabetic"] * 10)

chisq, pval = # survdiff(time, status, group)

chisq, pval = survdiff(time, status, group)
print("chi-square:", round(chisq, 2), " p-value:", round(pval, 4))
# chi-square: 13.7  p-value: 0.0002

The diabetic group had 8 deaths where about 3.1 were expected under equal risk, and the non-diabetic group had 5 where about 9.9 were expected. The chi-square near 13.7 on 1 degree of freedom gives a p-value around 0.0002. That is strong evidence that survival on dialysis differed between the two groups, with diabetic patients dying sooner. To report by how much, you would fit a Cox model and quote the hazard ratio.

Common mistakes

• Dropping censored patients to get a clean dataset. Deleting everyone still alive keeps only the deaths and crashes the survival estimate downward. Censored patients carry real follow-up time and must stay in.
• Coding the status flag backwards. The convention is 1 for the event and 0 for censored. Swap them and Surv(time, status) treats survivors as deaths and deaths as survivors, flipping the whole curve. Check the event count against what you expect before trusting any output.
• Quoting a median from the thin tail. If the curve only reaches 0.5 when two or three patients remain at risk, that median rests on almost nothing. Report it with its wide interval, or say the median was not reliably reached.
• Reading the log-rank p-value as an effect size. A small p-value says the curves differ, not by how much. A near-significant p-value in a tiny study can still hide a large hazard ratio. Pair the test with a hazard ratio.
• Using the log-rank test when the curves cross. Crossing curves break the proportional-hazards idea the test rewards, and the early and late differences cancel. Look at the plot first; if the curves cross, the log-rank test can miss a real difference.

Tips

• Always plot the Kaplan-Meier curves before running any test. The plot tells you whether the curves separate cleanly, plateau, or cross, and that shapes which test and model are honest.
• Confirm your event coding by checking the number of events. In R, sum(status) should match the deaths you counted by hand.
• Add a number-at-risk row beneath the time axis when you present a curve. It tells the reader where the estimate is solid and where the tail is thin.
• The log-rank test compares whole curves, so it uses early and late deaths alike. When you expect the effect to act mostly early or mostly late, say so in advance rather than fishing after the plot.
• For the size of the difference, move to a Cox proportional-hazards model. The log-rank test answers "is there a difference"; the Cox hazard ratio answers "how big".