Kaplan–Meier estimator

Problem definition

Let ${\displaystyle \tau \ge 0}$ be a random variable as the time that passes between the start of the possible exposure period, ${\displaystyle t_0}$, and the time that the event of interest takes place, ${\displaystyle t_1}$. As indicated above, the goal is to estimate the survival function ${\displaystyle S}$ underlying ${\displaystyle \tau }$. Recall that this function is defined as

${\displaystyle S(t)=\mathrm{Prob}(\tau >t)}$, where ${\displaystyle t=0,1,\dots }$ is the time.

Let ${\displaystyle {\tau }_1,\dots ,{\tau }_n\ge 0}$ be independent, identically distributed random variables, whose common distribution is that of ${\displaystyle \tau }$: ${\displaystyle {\tau }_j}$ is the random time when some event ${\displaystyle j}$ happened. The data available for estimating ${\displaystyle S}$ is not ${\displaystyle ({\tau }_j{)}_{j=1,\dots ,n}}$, but the list of pairs ${\displaystyle (({\tilde{\tau }}_j,c_j){)}_{j=1,\dots ,n}}$ where for ${\displaystyle j\in [n]:=\{1,2,\dots ,n\}}$, ${\displaystyle c_j\ge 0}$ is a fixed, deterministic integer, the censoring time of event ${\displaystyle j}$ and ${\displaystyle {\tilde{\tau }}_j=min({\tau }_j,c_j)}$. In particular, the information available about the timing of event ${\displaystyle j}$ is whether the event happened before the fixed time ${\displaystyle c_j}$ and if so, then the actual time of the event is also available. The challenge is to estimate ${\displaystyle S(t)}$ given this data.

Derivation of the Kaplan–Meier estimator

Two derivations of the Kaplan–Meier estimator are shown. Both are based on rewriting the survival function in terms of what is sometimes called hazard, or mortality rates. However, before doing this it is worthwhile to consider a naive estimator.

A naive estimator

To understand the power of the Kaplan–Meier estimator, it is worthwhile to first describe a naive estimator of the survival function.

Fix ${\displaystyle k\in [n]:=\{1,\dots ,n\}}$ and let ${\displaystyle t>0}$. A basic argument shows that the following proposition holds:

Proposition 1: If the censoring time ${\displaystyle c_k}$ of event ${\displaystyle k}$ exceeds ${\displaystyle t}$ (${\displaystyle c_k\ge t}$), then ${\displaystyle {\tilde{\tau }}_k\ge t}$ if and only if ${\displaystyle {\tau }_k\ge t}$.

Let ${\displaystyle k}$ be such that ${\displaystyle c_k\ge t}$. It follows from the above proposition that

${\displaystyle \mathrm{Prob}({\tau }_k\ge t)=\mathrm{Prob}({\tilde{\tau }}_k\ge t).}$

Let ${\displaystyle X_k=\mathbb{I}({\tilde{\tau }}_k\ge t)}$ and consider only those ${\displaystyle k\in C(t):=\{k:c_k\ge t\}}$, i.e. the events for which the outcome was not censored before time ${\displaystyle t}$. Let ${\displaystyle m(t)=|C(t)|}$ be the number of elements in ${\displaystyle C(t)}$. Note that the set ${\displaystyle C(t)}$ is not random and so neither is ${\displaystyle m(t)}$. Furthermore, ${\displaystyle (X_k{)}_{k\in C(t)}}$ is a sequence of independent, identically distributed Bernoulli random variables with common parameter ${\displaystyle S(t)=\mathrm{Prob}(\tau \ge t)}$. Assuming that ${\displaystyle m(t)>0}$, this suggests to estimate ${\displaystyle S(t)}$ using

${\displaystyle {\hat{S}}_{\text{naive}}(t)=\frac{1}{m(t)}\sum _{k:c_k\ge t}{X}_k=\frac{|\{1\le k\le n:{\tilde{\tau }}_k\ge t\}|}{|\{1\le k\le n:c_k\ge t\}|}=\frac{|\{1\le k\le n:{\tilde{\tau }}_k\ge t\}|}{m(t)},}$

where the second equality follows because ${\displaystyle {\tilde{\tau }}_k\ge t}$ implies ${\displaystyle c_k\ge t}$, while the last equality is simply a change of notation.

The quality of this estimate is governed by the size of ${\displaystyle m(t)}$. This can be problematic when ${\displaystyle m(t)}$ is small, which happens, by definition, when a lot of the events are censored. A particularly unpleasant property of this estimator, that suggests that perhaps it is not the "best" estimator, is that it ignores all the observations whose censoring time precedes ${\displaystyle t}$. Intuitively, these observations still contain information about ${\displaystyle S(t)}$: For example, when for many events with ${\displaystyle c_k<t}$, ${\displaystyle {\tau }_k<c_k}$ also holds, we can infer that events often happen early, which implies that ${\displaystyle \mathrm{Prob}(\tau \le t)}$ is large, which, through ${\displaystyle S(t)=1-\mathrm{Prob}(\tau \le t)}$ means that ${\displaystyle S(t)}$ must be small. However, this information is ignored by this naive estimator. The question is then whether there exists an estimator that makes a better use of all the data. This is what the Kaplan–Meier estimator accomplishes. Note that the naive estimator cannot be improved when censoring does not take place; so whether an improvement is possible critically hinges upon whether censoring is in place.

The plug-in approach

By elementary calculations,

${\displaystyle \begin{array}{cc}S(t)& =\mathrm{Prob}(\tau >t\mid \tau >t-1)\mathrm{Prob}(\tau >t-1)\\ & =(1-\mathrm{Prob}(\tau \le t\mid \tau >t-1))\mathrm{Prob}(\tau >t-1)\\ & =(1-\mathrm{Prob}(\tau =t\mid \tau \ge t))\mathrm{Prob}(\tau >t-1)\\ & =q(t)S(t-1),\end{array}}$

where the second to last equality used that ${\displaystyle \tau }$ is integer valued and for the last line we introduced

${\displaystyle q(t)=1-\mathrm{Prob}(\tau =t\mid \tau \ge t).}$

By a recursive expansion of the equality ${\displaystyle S(t)=q(t)S(t-1)}$, we get

${\displaystyle S(t)=q(t)q(t-1)\cdots q(0).}$

Note that here ${\displaystyle q(0)=1-\mathrm{Prob}(\tau =0\mid \tau >-1)=1-\mathrm{Prob}(\tau =0)}$.

The Kaplan–Meier estimator can be seen as a "plug-in estimator" where each ${\displaystyle q(s)}$ is estimated based on the data and the estimator of ${\displaystyle S(t)}$ is obtained as a product of these estimates.

It remains to specify how ${\displaystyle q(s)=1-\mathrm{Prob}(\tau =s\mid \tau \ge s)}$ is to be estimated. By a similar reasoning that lead to the construction of the naive estimator above, for any ${\displaystyle k\in [n]}$ such that ${\displaystyle c_k>s}$, ${\displaystyle \mathrm{Prob}(\tau =s)=\mathrm{Prob}({\tilde{\tau }}_k=s)}$ and ${\displaystyle \mathrm{Prob}(\tau \ge s)=\mathrm{Prob}({\tilde{\tau }}_k\ge s)}$ both hold. Hence, for any ${\displaystyle k\in [n]}$ such that ${\displaystyle c_k>s}$,

${\displaystyle \mathrm{Prob}(\tau =s|\tau \ge s)=\mathrm{Prob}({\tilde{\tau }}_k=s)/\mathrm{Prob}({\tilde{\tau }}_k\ge s).}$

Hence we arrive at the estimator

${\displaystyle \hat{q}(s)=1-\frac{|\{1\le k\le n:c_k>s,{\tilde{\tau }}_k=s\}|}{|\{1\le k\le n:c_k>s,{\tilde{\tau }}_k\ge s\}|}}$

(think of estimating the numerator and denominator separately in the definition of the "hazard rate" ${\displaystyle \mathrm{Prob}(\tau =s|\tau \ge s)}$). The Kaplan–Meier estimator is then given by

${\displaystyle \hat{S}(t)=\prod _{s=0}^t\hat{q}(s).}$

The form of the estimator stated at the beginning of the article can be obtained by some further algebra. For this, write ${\displaystyle \hat{q}(s)=1-d(s)/n(s)}$ where, using the actuarial science terminology, ${\displaystyle d(s)=|\{1\le k\le n:c_k>s,{\tilde{\tau }}_k=s\}|}$ is the number of known deaths at time ${\displaystyle s}$, while ${\displaystyle n(s)=|\{1\le k\le n:c_k>s,{\tilde{\tau }}_k\ge s\}|}$ is the number of those persons who are alive (and not being censored) at time ${\displaystyle s}$.

Note that if ${\displaystyle d(s)=0}$, ${\displaystyle \hat{q}(s)=1}$. This implies that we can leave out from the product defining ${\displaystyle \hat{S}(t)}$ all those terms where ${\displaystyle d(s)=0}$. Then, letting ${\displaystyle 0\le t_1<t_2<\cdots <t_m}$ be the times ${\displaystyle s}$ when ${\displaystyle d(s)>0}$, ${\displaystyle d_i=d(t_i)}$ and ${\displaystyle n_i=n(t_i)}$, we arrive at the form of the Kaplan–Meier estimator given at the beginning of the article:

${\displaystyle \hat{S}(t)=\prod _{i:t_i\le t}\left(1-\frac{d_i}{n_i}\right).}$

As opposed to the naive estimator, this estimator can be seen to use the available information more effectively: In the special case mentioned beforehand, when there are many early events recorded, the estimator will multiply many terms with a value below one and will thus take into account that the survival probability cannot be large.

Derivation as a maximum likelihood estimator

Kaplan–Meier estimator can be derived from maximum likelihood estimation of the discrete hazard function. More specifically given ${\displaystyle d_i}$ as the number of events and ${\displaystyle n_i}$ the total individuals at risk at time ${\displaystyle t_i}$, discrete hazard rate ${\displaystyle h_i}$ can be defined as the probability of an individual with an event at time ${\displaystyle t_i}$. Then survival rate can be defined as:

${\displaystyle S(t)=\prod _{i:t_i\le t}(1-h_i)}$

and the likelihood function for the hazard function up to time ${\displaystyle t_i}$ is:

${\displaystyle \mathcal{L}(h_{j:j\le i}\mid d_{j:j\le i},n_{j:j\le i})=\prod _{j=1}^i{h}_j^{d_j}(1-h_j{)}^{n_j-d_j}(\genfrac{}{}{0ex}{}{n_j}{d_j})}$

therefore the log likelihood will be:

${\displaystyle \log (\mathcal{L})=\sum _{j=1}^i\left(d_j\log (h_j)+(n_j-d_j)\log (1-h_j)+\log (\genfrac{}{}{0ex}{}{n_j}{d_j})\right)}$

finding the maximum of log likelihood with respect to ${\displaystyle h_i}$ yields:

${\displaystyle \frac{\mathrm{\partial }\log (\mathcal{L})}{\mathrm{\partial }{h}_i}=\frac{d_i}{{\hat{h}}_i}-\frac{n_i-d_i}{1-{\hat{h}}_i}=0\Rightarrow {\hat{h}}_i=\frac{d_i}{n_i}}$

where hat is used to denote maximum likelihood estimation. Given this result, we can write:

${\displaystyle \hat{S}(t)=\prod _{i:t_i\le t}\left(1-{\hat{h}}_i\right)=\prod _{i:t_i\le t}\left(1-\frac{d_i}{n_i}\right)}$

More generally (for continuous as well as discrete survival distributions), the Kaplan-Meier estimator may be interpreted as a nonparametric maximum likelihood estimator.

Benefits and limitations

The Kaplan–Meier estimator is one of the most frequently used methods of survival analysis. The estimate may be useful to examine recovery rates, the probability of death, and the effectiveness of treatment. It is limited in its ability to estimate survival adjusted for covariates; parametric survival models and the Cox proportional hazards model may be useful to estimate covariate-adjusted survival.

The Kaplan-Meier estimator is directly related to the Nelson-Aalen estimator and both maximize the empirical likelihood.

Statistical considerations

The Kaplan–Meier estimator is a statistic, and several estimators are used to approximate its variance. One of the most common estimators is Greenwood's formula:

${\displaystyle \widehat{\mathrm{Var}}\left(\hat{S}(t)\right)=\hat{S}(t{)}^2\sum _{i:t_i\le t}\frac{d_i}{n_i(n_i-d_i)},}$

where ${\displaystyle d_i}$ is the number of cases and ${\displaystyle n_i}$ is the total number of observations, for ${\displaystyle t_i<t}$.

In some cases, one may wish to compare different Kaplan–Meier curves. This can be done by the log rank test, and the Cox proportional hazards test.

Other statistics that may be of use with this estimator are pointwise confidence intervals, the Hall-Wellner band and the equal-precision band.

Software

• Mathematica: the built-in function SurvivalModelFit creates survival models.
• SAS: The Kaplan–Meier estimator is implemented in the proc lifetest procedure.
• R: the Kaplan–Meier estimator is available as part of the survival package.
• Stata: the command sts returns the Kaplan–Meier estimator.
• Python: the lifelines and scikit-survival packages each include the Kaplan–Meier estimator.
• MATLAB: the ecdf function with the 'function','survivor' arguments can calculate or plot the Kaplan–Meier estimator.
• StatsDirect: The Kaplan–Meier estimator is implemented in the Survival Analysis menu.
• SPSS: The Kaplan–Meier estimator is implemented in the Analyze > Survival > Kaplan-Meier... menu.
• Julia: the Survival.jl package includes the Kaplan–Meier estimator.
• Epi Info: Kaplan–Meier estimator survival curves and results for the log rank test are obtained with the KMSURVIVAL command.