Survival Analysis and Kaplan-Meier Curves
Introduction and Background
The year 2020 shall go down in history as a year in which a highly contagious and never-before-seen coronavirus spread across the world with devastating consequences. For the many who were unfortunate enough to contract this potentially deadly virus (now and in the future), likely one of the most important questions on their minds was whether or not they will be OK. However, if an infected patient were to ask their doctor this question, it's highly likely that the answer provided was not a definitive Yes or No. It may have been a "It should be all right," or "I expect things to improve," or possibly even "we're doing the best we can." To the patient, this non-binary sort-of-an-answer might be frustrating to hear, but it really is likely the best possible answer that can be given at the time. This is because biology (and biologic processes, in general) are predominantly stochastic in nature. In other words, they are not deterministic. This kind of non-deterministic behavior is certainly not unique to medicine and biology, but rather arises in many aspects of the world we live in. What's more, depending on the circumstances, different styles of mathematical analysis may be necessary in order to understand or characterize the behavior of the scenario at hand. In the realm of medicine and engineering, arguably the most common method of approach for these kinds of problems is survival analysis, alternatively known as time-to-failure analysis.
What is Survival Analysis and how does it work? At its core, survival analysis is a sub-branch of Statistics that is used for analyzing the amount of time that passes before an event occurs. When applied to biological organisms, as in medical applications, this event is often death (though not always). For non-biological systems (i.e. mechanical systems), this event is often some kind of mechanical failure. The underlying principle behind this field of analysis is the assumption of a fundamental conditional function that relates the probability that a given patient/device/object will survive without experiencing an event past a given amount of time. This fundamental function, or relationship, is referred to as the survival function. While there are a number of statistical arguments justifying the fundamental existence of these survival functions (and their appropriateness in describing these kinds of phenomena), I will refrain from delving into this topic. There are several methods/models that are used to estimate the survival function of time-to-event data; in this document, I will focus exclusively on one of the more common techniques, and explain in detail how it works, why it has such widespread acclaim, and some of the limitations that are inherent in its use.
The Kaplan-Meier Estimator
The Kaplan-Meier Estimator emerged from the groundbreaking paper on Nonparametric Estimation From Incomplete Observations by Edward Kaplan and Paul Meier in 1958 (link). In terms of impact, this paper is one of the most cited scientific papers of all time, and was ranked 11th in 2014 with approximately 38,600 citations, according to an article published in the journal Nature (link). The formula that was proposed in this groundbreaking paper (and later named after the two authors) suggested that the Survival function \(S(t)\) could be non-parametrically estimated in the following way:
\[ \begin{array}{rcl} \hat{S}(t) &=& \displaystyle \prod_{i:\ t_i\le t}\bigg(1-\frac{d_i}{n_i}\bigg) \end{array} \]
Where \(t_i\) is a time in which at least one event occurred, \(d_i\) is the number of events that occurred at \(t_i\), and \(n_i\) is the number of objects of interest that have been known to have not had an event (i.e. survived) up to \(t_i\). As evident from the formula, this is an example of a non-parametric statistic. To best get a feel for how this formula works, it is useful to examine this formula in the context of a simple example.
| ID | Time | Event |
|---|---|---|
| ID1 | 1 | 1 |
| ID2 | 2 | 1 |
| ID3 | 3 | 1 |
| ID4 | 4 | 1 |
| ID5 | 5 | 1 |
In this simple example, there are five participants/objects in the study. The time-to-event duration of ID1 was 1 (unit of time). For ID2, it was 2, and so on. From this time-to-event data, we can calculate \(n_i\) and \(d_i\) is a very simple fashion. From these inputs, the Kaplan-Meier estimator is determined to be the following:
Notice that at each step in above graphic, the survival function dropped by an amount that was equal to the current survival function value multiplied by \(d_i/n_i\). At \(t=1\), \(n_i=5\), \(d_i=1\) and \(\hat{S}(t)=1\), and thus the function dropped by \(1\times 1/5 = 0.2\). At \(t=2\), \(n_i=4\), \(d_i=1\) and \(\hat{S}(t)=0.8\), and thus the function dropped by \(0.8\times 1/4 = 0.2\). At \(t=3\), the function dropped by \(0.6\times 1/3 = 0.2\). And so on. Take a moment to go back to the formula and verify that this behavior is indeed reflected by the provided formulation.
The insight to be gained from this is the fact that the survival function counts all participants who have survived up until a given time/duration \(t\) (as opposed to only those participants at or around time \(t\)), only needs to be calculated for the points in which there was an event (because whenever there are no events, the survival function drops by \(\hat{S}(t)\times 0/n_i = 0\) and thus can be completely ignored), it is cumulative (ensuring that all data points can have an effect on the estimated result, suggesting that the information that is available in the data is used effectively), and it is quite easy to calculate.
Why does this product of ratios provide accurate estimates of the survival probability larger than some given duration of time? The simplest explanation of this can be described from the perspective of conditional probability. Suppose that you are trying to find \(\hat{S}(3) = \mathsf{P}(T > 3)\), where \(T\) is the true/unknown time-to-death of the object of interest. By using conditional probability, we can expand this out in the following way.
\[
\begin{array}{rcll}
\hat{S}(3) &=& \mathsf{P}(T > 3)
\\
&=& \mathsf{P}(T > 3\ \cap \ T>2) & \text{Since }\ \mathsf{P}(T > 3) = \mathsf{P}(T > 3 \ \cap \ T>2)
\\
&=& \mathsf{P}(T > 3\ |\ T>2) \times \mathsf{P}(T>2) & \text{Bayes Theorem}
\\
&=& \mathsf{P}(T > 3\ |\ T>2) \times \mathsf{P}(T>2\ |\ T>1.5) \times \mathsf{P}(T>1.5) & ...
\end{array}
\]
The fact that the survival probability at 3 years of age can be written to be proportional to that at 2 years of age can be understood as reflecting the fact that the object can only survive 3 years of age after first surviving 2 years of age. One can then repeat this argument as much as necessary to obtain the desired extended product as the result. With this in mind, the original formula should make a lot more sense.
With this introductory example out of the way, we can dive a little deeper into the nuances of this formulation, as well as the reasons why it is one of the most popular approaches to survival analysis. Let's start with the primary problem that Kaplan and Meier built this estimator to combat; namely missing and incomplete time-to-event data.
For any kind of time-to-event scenario, there are always a number of options that a given patient/object-of-interest may be characterized by at a given moment in time (as relevant to the analysis):
- The event has happened, and the time-to-event duration is known (complete data)
- The event has not happened (yet), and the duration of time without an event is known (right-censored data)
- The event has not happened and never will happen due to removal or lack of communication/knowledge, and the known/confirmed duration of time without an event (i.e. time-to-last communication/knowledge) is known (right-censored data)
NOTE: Other possible options, such as interval censoring, are not considered here.
The power behind the Kaplan-Meier Estimator is that unlike other more naive survival function estimation approaches, which may only be able to reliable extract information from data sources/records that are 'complete', the Kaplan Meier estimator can be proven to be able to effectively make use of both complete and right-censored data to craft the estimated survival function. Specifically, this ensures that none of the data is thrown out due to the incompleteness of the data observations. This, in combination with its status as a non-parametric estimator, is arguably the reason behind its' profound popularity in practical applications (particularly medicine research), as more likely than not you will never have time-to-event datasets that are complete (in the sense that all records are complete). As a simple, but profound, example, suppose you want to better understand human lifetimes (i.e. time-to-death). If you don't throw out any potential data, you will never have a complete dataset of time-to-death data, simply because you yourself are still alive (and thus your time-to-death duration is not known). Waiting until such a time that one has complete data in this case (and in many others) is thus completely ridiculous. Rather, if the approach being used truly does require nothing but complete data, then all the data records with incomplete information would need to be removed from consideration, reducing the size of the data sample and, as a consequence, the accuracy of the final result. However, with the Kaplan-Meier Estimator, this is not an issue, as it was specifically formulated in such a way that even incomplete data records can be effectively made use of.
In fact, with respect to the estimated survival function using the Kaplan-Meier Estimator, none of these censored data records actually cause the survival curve to drop, as we have previously demonstrated that the only 'drops' in the survival curve occur when there are a non-zero number of 'events' at some time/duration \(t_i\). However, they do still play a role in the overall calculation by affecting the value of \(n_i\), or the number of objects of interest that have been observed to have survived without an event for a duration longer than \(t_i\). This is because, while the true time-to-event duration of these censored objects may not be known, it is known that they must be larger than the current age of the objects. As a demonstration of the effect of these censored data records, if we modify our original simple example such that the ID3 is no longer an 'event', but rather a right-censored data point, then the result is as follows:
| ID | Time | Event |
|---|---|---|
| ID1 | 1 | 1 |
| ID2 | 2 | 1 |
| ID3 | 3 | 0 |
| ID4 | 4 | 1 |
| ID5 | 5 | 1 |
Here, the right-censored data record is specified by ensuring the value in the Event column associated with this record (namely ID3) is zero.
There are a few things to note here. First, at \(t=3\), there was no drop in the Kaplan-Meier survival function estimate. This is because there were no events at \(t=3\) (i.e. when \(t_i=3\), \(d_i=0\), as can be seen in the sub-graph for \(d_i\)). Second, there is no change in the graph of \(n_i\) for this slightly modified sample dataset. Third, the reduction in the survival function was larger at \(t=4\) than it was at \(t=2\). Specifically, when \(t_i=4\), the survival function at this time was \(\hat{S}(t)=0.6\), with \(d_i=1\) and \(n_i=2\). Plugging in, we find that the overall size of the drop should be equal to \(0.6\times1/2=0.3\). It's important to note that the information provided by the censored data record behaved exactly like the previous version (in which all time-to-event measurements were known) up to the point where no further information was known regarding the time-to-event of the given record (i.e. it is know to be larger than 3, but nothing more than that). At the same time, no assumptions regarding the true time-to-event duration are made, which is why the curve differed slightly after \(t=3\) compared to the original version (because of this incomplete/unknown information). This demonstrates how useful information can be extracted even from incomplete data records using the Kaplan-Meier Estimator in such a way that no implicit assumptions regarding unknown information are made.
Practical Considerations
In addition to all that we have discussed up to this point, there are also a number of practical considerations that one should be aware of when using and applying the Kaplan-Meier Estimator to time-to-event data. The first that, like all mathematical models, the sample size and quality of the data matter. To be more mathematically precise, the confidence of the resulting curve fit will depend on the data that is used. It is standard practice to use a 95% confidence level when constructing confidence intervals for the estimated fit. The exact details regarding how these confidence intervals are constructed lie outside the scope of this document; however, at a high-level, larger datasets with few censored data records have the tightest confidence intervals (i.e. low uncertainty), while small datasets with many censored data records have the largest confidence intervals (i.e. high uncertainty). This is demonstrated in the example below, using time-to-event survival data drawn from a Geometric distribution and censoring determined by a Binomial distribution with probability \(p\).
\[ \begin{array}{rcl} T\sim \text{Geom}(0.05) && E_{vent} \sim B(1, p) \end{array} \]
Conclusion
Just as any other mathematical model, there are pros and cons associated with using the Kaplan-Meier Estimator to estimate the survival function of some phenomenon of interest. The fact that it is a non-parametric approach to estimating the survival function is something that makes it particularly desirable in many cases for a number of reasons. As a non-parametric estimator, it makes very few implicit assumptions regarding the phenomenon of interest, and thus may be used for a wide range of time-to-event applications. Performance-wise, it is also quite robust and is unbiased, at least up to the last observed event in the provided data (note that this may not be the maximum observed duration in the data if the records with the largest duration are censored). As an additional perk, the Kaplan-Meier Estimator can be demonstrated to be the Maximum-Likelihood Estimator for these kinds of time-to-event analyses (specifically, it is the MLE of the Hazard Function).
Despite all this, the Kaplan-Meier Estimator may not always be the optimal choice for a given time-to-event application. For example, there are many applications in the world of engineering in which enough is known about the problem at hand such that parametric approaches to this kind of survival analysis are not only possible, but also have better performance compared to their non-parametric counterparts. Of particularly note in this respect is the improvement in statistical 'Power' that many parametric approaches provide even with the same input data. What's more, the Kaplan-Meier Estimator is not particular flexible with respect to adjusting for covariates that may affect the underlying survival function of interest. And lastly, the Kaplan-Meier Estimator is not able to account for censored data records that are not right-censored (e.g. left or interval-censored data).
In conclusion, the Kaplan-Meier estimate is a powerful tool that can be used to develop an estimate of the survival function defining a specific process or phenomenon when dealing with time-to-event data. Being able to estimate the survival function is very useful and powerful, as it can be used for a wide range of applications, including predicting future performance (such as the mechanical failure rate of servers in a datacenter, or the maintenance needs of a number of assets in the next two years), or comparing the survival rates between different data sub-populations (such as the long-term performance of two different vaccines with respect to conferring immunity or resistance to a pathogen).
## R Session Information:
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19041)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United States.1252
## [3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C
## [5] LC_TIME=English_United States.1252
##
## attached base packages:
## [1] grid stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] survminer_0.4.8 ggpubr_0.4.0 survival_3.2-7 reshape2_1.4.4 gridExtra_2.3 dplyr_1.0.2
## [7] plyr_1.8.6 ggplot2_3.3.2
##
## loaded via a namespace (and not attached):
## [1] zoo_1.8-8 tidyselect_1.1.0 xfun_0.19 purrr_0.3.4 splines_4.0.3
## [6] haven_2.3.1 lattice_0.20-41 carData_3.0-4 colorspace_2.0-0 vctrs_0.3.5
## [11] generics_0.1.0 htmltools_0.5.0 yaml_2.2.1 survMisc_0.5.5 rlang_0.4.9
## [16] pillar_1.4.7 foreign_0.8-80 glue_1.4.2 withr_2.3.0 readxl_1.3.1
## [21] lifecycle_0.2.0 stringr_1.4.0 cellranger_1.1.0 munsell_0.5.0 ggsignif_0.6.0
## [26] gtable_0.3.0 zip_2.1.1 evaluate_0.14 labeling_0.4.2 knitr_1.30
## [31] rio_0.5.16 rmdformats_1.0.0 forcats_0.5.0 curl_4.3 highr_0.8
## [36] broom_0.7.2 Rcpp_1.0.5 xtable_1.8-4 scales_1.1.1 backports_1.2.0
## [41] abind_1.4-5 km.ci_0.5-2 farver_2.0.3 hms_0.5.3 digest_0.6.27
## [46] openxlsx_4.2.3 stringi_1.5.3 bookdown_0.21 rstatix_0.6.0 KMsurv_0.1-5
## [51] tools_4.0.3 magrittr_2.0.1 tibble_3.0.4 crayon_1.3.4 tidyr_1.1.2
## [56] car_3.0-10 pkgconfig_2.0.3 ellipsis_0.3.1 Matrix_1.2-18 data.table_1.13.2
## [61] rmarkdown_2.5 R6_2.5.0 compiler_4.0.3