This articleneeds additional citations forverification. Please helpimprove this article byadding citations to reliable sources. Unsourced material may be challenged and removed. Find sources: "Survival analysis" – news ·newspapers ·books ·scholar ·JSTOR(April 2021) (Learn how and when to remove this message) |
Survival analysis is a branch ofstatistics for analyzing the expected duration of time until one event occurs, such as death inbiological organisms and failure in mechanical systems.[1] This topic is calledreliability theory,reliability analysis orreliability engineering inengineering,duration analysis orduration modelling ineconomics, andevent history analysis insociology. Survival analysis attempts to answer certain questions, such as what is the proportion of a population which will survive past a certain time? Of those that survive, at what rate will they die or fail? Can multiple causes of death or failure be taken into account? How do particular circumstances or characteristics increase or decrease the probability ofsurvival?
To answer such questions, it is necessary to define "lifetime". In the case of biological survival,death is unambiguous, but for mechanical reliability,failure may not be well-defined, for there may well be mechanical systems in which failure is partial, a matter of degree, or not otherwise localized intime. Even in biological problems, some events (for example,heart attack or other organ failure) may have the same ambiguity. Thetheory outlined below assumes well-defined events at specific times; other cases may be better treated by models which explicitly account for ambiguous events.
More generally, survival analysis involves the modelling of time to event data; in this context, death or failure is considered an "event" in the survival analysis literature – traditionally only a single event occurs for each subject, after which the organism or mechanism is dead or broken.Recurring event orrepeated event models relax that assumption. The study of recurring events is relevant insystems reliability, and in many areas of social sciences and medical research.
Survival analysis is used in several ways:
The following terms are commonly used in survival analyses:
This example uses theAcute Myelogenous Leukemia survival data set "aml" from the "survival" package in R. The data set is from Miller (1997)[2] and the question is whether the standard course of chemotherapy should be extended ('maintained') for additional cycles.
The aml data set sorted by survival time is shown in the box.
observation | time (weeks) | status | x |
---|---|---|---|
12 | 5 | 1 | Nonmaintained |
13 | 5 | 1 | Nonmaintained |
14 | 8 | 1 | Nonmaintained |
15 | 8 | 1 | Nonmaintained |
1 | 9 | 1 | Maintained |
16 | 12 | 1 | Nonmaintained |
2 | 13 | 1 | Maintained |
3 | 13 | 0 | Maintained |
17 | 16 | 0 | Nonmaintained |
4 | 18 | 1 | Maintained |
5 | 23 | 1 | Maintained |
18 | 23 | 1 | Nonmaintained |
19 | 27 | 1 | Nonmaintained |
6 | 28 | 0 | Maintained |
20 | 30 | 1 | Nonmaintained |
7 | 31 | 1 | Maintained |
21 | 33 | 1 | Nonmaintained |
8 | 34 | 1 | Maintained |
22 | 43 | 1 | Nonmaintained |
9 | 45 | 0 | Maintained |
23 | 45 | 1 | Nonmaintained |
10 | 48 | 1 | Maintained |
11 | 161 | 0 | Maintained |
The last observation (11), at 161 weeks, is censored. Censoring indicates that the patient did not have an event (no recurrence of aml cancer). Another subject, observation 3, was censored at 13 weeks (indicated by status=0). This subject was in the study for only 13 weeks, and the aml cancer did not recur during those 13 weeks. It is possible that this patient was enrolled near the end of the study, so that they could be observed for only 13 weeks. It is also possible that the patient was enrolled early in the study, but was lost to follow up or withdrew from the study. The table shows that other subjects were censored at 16, 28, and 45 weeks (observations 17, 6, and 9 with status=0). The remaining subjects all experienced events (recurrence of aml cancer) while in the study. The question of interest is whether recurrence occurs later in maintained patients than in non-maintained patients.
Thesurvival functionS(t), is the probability that a subject survives longer than timet.S(t) is theoretically a smooth curve, but it is usually estimated using theKaplan–Meier (KM) curve. The graph shows the KM plot for the aml data and can be interpreted as follows:
Alife table summarizes survival data in terms of the number of events and the proportion surviving at each event time point. The life table for the aml data, created using the R software, is shown.
time | n.risk | n.event | survival | std.err | lower 95% CI | upper 95% CI |
---|---|---|---|---|---|---|
5 | 23 | 2 | 0.913 | 0.0588 | 0.8049 | 1 |
8 | 21 | 2 | 0.8261 | 0.079 | 0.6848 | 0.996 |
9 | 19 | 1 | 0.7826 | 0.086 | 0.631 | 0.971 |
12 | 18 | 1 | 0.7391 | 0.0916 | 0.5798 | 0.942 |
13 | 17 | 1 | 0.6957 | 0.0959 | 0.5309 | 0.912 |
18 | 14 | 1 | 0.646 | 0.1011 | 0.4753 | 0.878 |
23 | 13 | 2 | 0.5466 | 0.1073 | 0.3721 | 0.803 |
27 | 11 | 1 | 0.4969 | 0.1084 | 0.324 | 0.762 |
30 | 9 | 1 | 0.4417 | 0.1095 | 0.2717 | 0.718 |
31 | 8 | 1 | 0.3865 | 0.1089 | 0.2225 | 0.671 |
33 | 7 | 1 | 0.3313 | 0.1064 | 0.1765 | 0.622 |
34 | 6 | 1 | 0.2761 | 0.102 | 0.1338 | 0.569 |
43 | 5 | 1 | 0.2208 | 0.0954 | 0.0947 | 0.515 |
45 | 4 | 1 | 0.1656 | 0.086 | 0.0598 | 0.458 |
48 | 2 | 1 | 0.0828 | 0.0727 | 0.0148 | 0.462 |
The life table summarizes the events and the proportion surviving at each event time point. The columns in the life table have the following interpretation:
Thelog-rank test compares the survival times of two or more groups. This example uses a log-rank test for a difference in survival in the maintained versus non-maintained treatment groups in the aml data. The graph shows KM plots for the aml data broken out by treatment group, which is indicated by the variable "x" in the data.
The null hypothesis for a log-rank test is that the groups have the same survival. The expected number of subjects surviving at each time point in each is adjusted for the number of subjects at risk in the groups at each event time. The log-rank test determines if the observed number of events in each group is significantly different from the expected number. The formal test is based on a chi-squared statistic. When the log-rank statistic is large, it is evidence for a difference in the survival times between the groups. The log-rank statistic approximately has aChi-squared distribution with one degree of freedom, and thep-value is calculated using theChi-squared test.
For the example data, the log-rank test for difference in survival gives a p-value of p=0.0653, indicating that the treatment groups do not differ significantly in survival, assuming an alpha level of 0.05. The sample size of 23 subjects is modest, so there is littlepower to detect differences between the treatment groups. The chi-squared test is based on asymptotic approximation, so the p-value should be regarded with caution for smallsample sizes.
Kaplan–Meier curves and log-rank tests are most useful when the predictor variable is categorical (e.g., drug vs. placebo), or takes a small number of values (e.g., drug doses 0, 20, 50, and 100 mg/day) that can be treated as categorical. The log-rank test and KM curves don't work easily with quantitative predictors such as gene expression, white blood count, or age. For quantitative predictor variables, an alternative method isCox proportional hazards regression analysis. Cox PH models work also with categorical predictor variables, which are encoded as {0,1} indicator or dummy variables. The log-rank test is a special case of a Cox PH analysis, and can be performed using Cox PH software.
This example uses the melanoma data set from Dalgaard Chapter 14.[3]
Data are in the R package ISwR. The Cox proportional hazards regression using R gives the results shown in the box.
The Cox regression results are interpreted as follows.
The summary output also gives upper and lower 95% confidence intervals for the hazard ratio: lower 95% bound = 1.15; upper 95% bound = 3.26.
Finally, the output gives p-values for three alternative tests for overall significance of the model:
These three tests are asymptotically equivalent. For large enough N, they will give similar results. For small N, they may differ somewhat. The last row, "Score (logrank) test" is the result for the log-rank test, with p=0.011, the same result as the log-rank test, because the log-rank test is a special case of a Cox PH regression. The Likelihood ratio test has better behavior for small sample sizes, so it is generally preferred.
The Cox model extends the log-rank test by allowing the inclusion of additional covariates.[4] This example use the melanoma data set where the predictor variables include a continuous covariate, the thickness of the tumor (variable name = "thick").
In the histograms, the thickness values arepositively skewed and do not have aGaussian-like,Symmetric probability distribution. Regression models, including the Cox model, generally give more reliable results with normally-distributed variables.[citation needed] For this example we may use alogarithmic transform. The log of the thickness of the tumor looks to be more normally distributed, so the Cox models will use log thickness. The Cox PH analysis gives the results in the box.
The p-value for all three overall tests (likelihood, Wald, and score) are significant, indicating that the model is significant. The p-value for log(thick) is 6.9e-07, with a hazard ratio HR = exp(coef) = 2.18, indicating a strong relationship between the thickness of the tumor and increased risk of death.
By contrast, the p-value for sex is now p=0.088. The hazard ratio HR = exp(coef) = 1.58, with a 95% confidence interval of 0.934 to 2.68. Because the confidence interval for HR includes 1, these results indicate that sex makes a smaller contribution to the difference in the HR after controlling for the thickness of the tumor, and only trend toward significance. Examination of graphs of log(thickness) by sex and a t-test of log(thickness) by sex both indicate that there is a significant difference between men and women in the thickness of the tumor when they first see the clinician.
The Cox model assumes that the hazards are proportional. The proportional hazard assumption may be tested using the R function cox.zph(). A p-value which is less than 0.05 indicates that the hazards are not proportional. For the melanoma data we obtain p=0.222. Hence, we cannot reject the null hypothesis of the hazards being proportional. Additional tests and graphs for examining a Cox model are described in the textbooks cited.
Cox models can be extended to deal with variations on the simple analysis.
The Cox PH regression model is a linear model. It is similar to linear regression and logistic regression. Specifically, these methods assume that a single line, curve, plane, or surface is sufficient to separate groups (alive, dead) or to estimate a quantitative response (survival time).
In some cases alternative partitions give more accurate classification or quantitative estimates. One set of alternative methods are tree-structured survival models,[5][6][7] including survival random forests.[8] Tree-structured survival models may give more accurate predictions than Cox models. Examining both types of models for a given data set is a reasonable strategy.
This example of a survival tree analysis uses the R package "rpart".[9] The example is based on 146 stage C prostate cancer patients in the data set stagec in rpart. Rpart and the stagec example are described in Atkinson and Therneau (1997),[10] which is also distributed as a vignette of the rpart package.[9]
The variables in stages are:
The survival tree produced by the analysis is shown in the figure.
Each branch in the tree indicates a split on the value of a variable. For example, the root of the tree splits subjects with grade < 2.5 versus subjects with grade 2.5 or greater. The terminal nodes indicate the number of subjects in the node, the number of subjects who have events, and the relative event rate compared to the root. In the node on the far left, the values 1/33 indicate that one of the 33 subjects in the node had an event, and that the relative event rate is 0.122. In the node on the far right bottom, the values 11/15 indicate that 11 of 15 subjects in the node had an event, and the relative event rate is 2.7.
An alternative to building a single survival tree is to build many survival trees, where each tree is constructed using a sample of the data, and average the trees to predict survival.[8] This is the method underlying the survival random forest models. Survival random forest analysis is available in the R package "randomForestSRC".[11]
The randomForestSRC package includes an example survival random forest analysis using the data set pbc. This data is from the Mayo Clinic Primary Biliary Cirrhosis (PBC) trial of the liver conducted between 1974 and 1984. In the example, the random forest survival model gives more accurate predictions of survival than the Cox PH model. The prediction errors are estimated bybootstrap re-sampling.
Recent advancements in deep representation learning have been extended to survival estimation. The DeepSurv[12] model proposes to replace the log-linear parameterization of the CoxPH model with a multi-layer perceptron. Further extensions like Deep Survival Machines[13] and Deep Cox Mixtures[14] involve the use of latent variable mixture models to model the time-to-event distribution as a mixture of parametric or semi-parametric distributions while jointly learning representations of the input covariates. Deep learning approaches have shown superior performance especially on complex input data modalities such as images and clinical time-series.
This sectiondoes notcite anysources. Please helpimprove this section byadding citations to reliable sources. Unsourced material may be challenged andremoved.(April 2021) (Learn how and when to remove this message) |
The object of primary interest is thesurvival function, conventionally denotedS, which is defined as
wheret is some time,T is arandom variable denoting the time of death, and "Pr" stands forprobability. That is, the survival function is the probability that the time of death is later than some specified timet.The survival function is also called thesurvivor function orsurvivorship function in problems of biological survival, and thereliability function in mechanical survival problems. In the latter case, the reliability function is denotedR(t).
Usually one assumesS(0) = 1, although it could be less than 1 if there is the possibility of immediate death or failure.
The survival function must be non-increasing:S(u) ≤S(t) ifu ≥t. This property follows directly becauseT>u impliesT>t. This reflects the notion that survival to a later age is possible only if all younger ages are attained. Given this property, the lifetime distribution function and event density (F andf below) are well-defined.
The survival function is usually assumed to approach zero as age increases without bound (i.e.,S(t) → 0 ast → ∞), although the limit could be greater than zero if eternal life is possible. For instance, we could apply survival analysis to a mixture of stable and unstablecarbon isotopes; unstable isotopes would decay sooner or later, but the stable isotopes would last indefinitely.
Related quantities are defined in terms of the survival function.
Thelifetime distribution function, conventionally denotedF, is defined as the complement of the survival function,
IfF isdifferentiable then the derivative, which is the density function of the lifetime distribution, is conventionally denotedf,
The functionf is sometimes called theevent density; it is the rate of death or failure events per unit time.
The survival function can be expressed in terms ofprobability distribution andprobability density functions
Similarly, a survival event density function can be defined as
In other fields, such as statistical physics, the survival event density function is known as thefirst passage time density.
Thehazard function is defined as the event rate at time conditional on survival at time
Synonyms forhazard function in different fields include hazard rate,force of mortality (demography andactuarial science, denoted by), force of failure, orfailure rate (engineering, denoted). For example, in actuarial science, denotes rate of death for people aged, whereas inreliability engineering denotes rate of failure of components after operation for time.
Suppose that an item has survived for a time and we desire the probability that it will not survive for an additional time:
Any function is a hazard function if and only if it satisfies the following properties:
In fact, the hazard rate is usually more informative about the underlying mechanism of failure than the other representations of a lifetime distribution.
The hazard function must be non-negative,, and its integral over must be infinite, but is not otherwise constrained; it may be increasing or decreasing, non-monotonic, or discontinuous. An example is thebathtub curve hazard function, which is large for small values of, decreasing to some minimum, and thereafter increasing again; this can model the property of some mechanical systems to either fail soon after operation, or much later, as the system ages.
The hazard function can alternatively be represented in terms of thecumulative hazard function, conventionally denoted or:
so transposing signs and exponentiating
or differentiating (with the chain rule)
The name "cumulative hazard function" is derived from the fact that
which is the "accumulation" of the hazard over time.
From the definition of, we see that it increases without bound ast tends to infinity (assuming that tends to zero). This implies that must not decrease too quickly, since, by definition, the cumulative hazard has to diverge. For example, is not the hazard function of any survival distribution, because its integral converges to 1.
The survival function, the cumulative hazard function, the density, the hazard function, and the lifetime distribution function are related through
Future lifetime at a given time is the time remaining until death, given survival to age. Thus, it is in the present notation. Theexpected future lifetime is theexpected value of future lifetime. The probability of death at or before age, given survival until age, is just
Therefore, the probability density of future lifetime is
and the expected future lifetime is
where the second expression is obtained usingintegration by parts.
For, that is, at birth, this reduces to the expected lifetime.
In reliability problems, the expected lifetime is called themean time to failure, and the expected future lifetime is called themean residual lifetime.
As the probability of an individual surviving until aget or later isS(t), by definition, the expected number of survivors at aget out of an initialpopulation ofn newborns isn ×S(t), assuming the same survival function for all individuals. Thus the expected proportion of survivors isS(t).If the survival of different individuals is independent, the number of survivors at aget has abinomial distribution with parametersn andS(t), and thevariance of the proportion of survivors isS(t) × (1-S(t))/n.
The age at which a specified proportion of survivors remain can be found by solving the equationS(t) =q fort, whereq is thequantile in question. Typically one is interested in themedian lifetime, for whichq = 1/2, or other quantiles such asq = 0.90 orq = 0.99.
Censoring is a form of missing data problem in which time to event is not observed for reasons such as termination of study before all recruited subjects have shown the event of interest or the subject has left the study prior to experiencing an event. Censoring is common in survival analysis.
If only the lower limitl for the true event timeT is known such thatT >l, this is calledright censoring. Right censoring will occur, for example, for those subjects whose birth date is known but who are still alive when they arelost to follow-up or when the study ends. We generally encounter right-censored data.
If the event of interest has already happened before the subject is included in the study but it is not known when it occurred, the data is said to beleft-censored.[15] When it can only be said that the event happened between two observations or examinations, this isinterval censoring.
Left censoring occurs for example when a permanent tooth has already emerged prior to the start of a dental study that aims to estimate its emergence distribution. In the same study, an emergence time is interval-censored when the permanent tooth is present in the mouth at the current examination but not yet at the previous examination. Interval censoring often occurs in HIV/AIDS studies. Indeed, time to HIV seroconversion can be determined only by a laboratory assessment which is usually initiated after a visit to the physician. Then one can only conclude that HIV seroconversion has happened between two examinations. The same is true for the diagnosis of AIDS, which is based on clinical symptoms and needs to be confirmed by a medical examination.
It may also happen that subjects with a lifetime less than some threshold may not be observed at all: this is calledtruncation. Note that truncation is different from left censoring, since for a left censored datum, we know the subject exists, but for a truncated datum, we may be completely unaware of the subject. Truncation is also common. In a so-calleddelayed entry study, subjects are not observed at all until they have reached a certain age. For example, people may not be observed until they have reached the age to enter school. Any deceased subjects in the pre-school age group would be unknown. Left-truncated data are common inactuarial work forlife insurance andpensions.[16]
Left-censored data can occur when a person's survival time becomes incomplete on the left side of the follow-up period for the person. For example, in an epidemiological example, we may monitor a patient for an infectious disorder starting from the time when he or she is tested positive for the infection. Although we may know the right-hand side of the duration of interest, we may never know the exact time of exposure to the infectious agent.[17]
Survival models can be usefully viewed as ordinary regression models in which the response variable is time. However, computing the likelihood function (needed for fitting parameters or making other kinds of inferences) is complicated by the censoring. Thelikelihood function for a survival model, in the presence of censored data, is formulated as follows. By definition the likelihood function is theconditional probability of the data given the parameters of the model.It is customary to assume that the data are independent given the parameters. Then the likelihood function is the product of the likelihood of each datum. It is convenient to partition the data into four categories: uncensored, left censored, right censored, and interval censored. These are denoted "unc.", "l.c.", "r.c.", and "i.c." in the equation below.
For uncensored data, with equal to the age at death, we have
For left-censored data, such that the age at death is known to be less than, we have
For right-censored data, such that the age at death is known to be greater than, we have
For an interval censored datum, such that the age at death is known to be less than and greater than, we have
An important application where interval-censored data arises is current status data, where an event is known not to have occurred before an observation time and to have occurred before the next observation time.
TheKaplan–Meier estimator can be used to estimate the survival function. TheNelson–Aalen estimator can be used to provide anon-parametric estimate of the cumulative hazard rate function. These estimators require lifetime data. Periodic case (cohort) and death (and recovery) counts are statistically sufficient to make nonparametric maximum likelihood and least squares estimates of survival functions, without lifetime data.
While many parametric models assume a continuous-time, discrete-time survival models can be mapped to a binary classification problem. In a discrete-time survival model the survival period is artificially resampled in intervals where for each interval a binary target indicator is recorded if the event takes place in a certain time horizon.[18] If a binary classifier (potentially enhanced with a different likelihood to take more structure of the problem into account) iscalibrated, then the classifier score is the hazard function (i.e. the conditional probability of failure).[18]
Discrete-time survival models are connected toempirical likelihood.[19][20]
The goodness of fit of survival models can be assessed usingscoring rules.[21]
The textbook by Kleinbaum has examples of survival analyses using SAS, R, and other packages.[22] The textbooks by Brostrom,[23] Dalgaard[3]and Tableman and Kim[24]give examples of survival analyses using R (or using S, and which run in R).