Recommended Words
Relatively hardcore. If you can’t get into it, you can take your time to read it later.
Main Text
This article is transcoded by 简悦 SimpRead, original address thisis.yorven.site
In biomedical research, survival analysis is a very important and common analytical method. This article provides a brief but detailed overview of the Kaplan–Meier model and Cox proportional hazards model in survival analysis, helping readers better understand survival analysis and related concepts. This article is suitable for beginners in biomedical majors as well as non-professionals interested in survival analysis.
Survival Analysis
First, let’s briefly describe the scenarios where survival analysis is used. Those already familiar can skip this part. Survival analysis is often used in the study of diseases such as cancer. For example, in a clinical trial of an anti-cancer drug, a group of cancer patients is screened and randomly divided into two groups: one group takes the trial drug, and the other group takes the control drug. After medication, the survival time of each patient from taking the drug until death is recorded. By examining whether there is a statistical difference in survival time between the two groups, the effectiveness of the trial drug can be judged. Here, death is the key observed event in the entire experiment, referred to as the event. For each patient, the exact time at which this event occurs needs to be recorded. Therefore, survival analysis can be abstractly described as studying whether there are differences in the time to occurrence of specific events under different conditions. These specific events can be death or tumor metastasis, recurrence, patient discharge, readmission, or any clearly identifiable events. Different conditions refer to different grouping bases, which can be age, gender, region, expression levels of certain genes, presence or absence of certain mutations, etc.
The following figure is from the research results mentioned by Academician Zhong Nanshan in the European Respiratory Society report on Covid-19. They analyzed the time from symptom onset to hospitalization for patients inside and outside Hubei Province. The onset of symptoms corresponds to the start event, hospitalization is the observed event just mentioned, and inside or outside Hubei Province is the grouping condition. The figure also mentions that they used the Cox model to adjust for geography, which will be covered later in this article. Readers interested in Zhong Nanshan’s report can visit this link.
Above is a brief introduction to the application scenarios of survival analysis; we believe readers now have a basic impression of survival analysis. Next, we need to discuss the issue of censored data.
Censored data literally means deleted or missing data. In survival analysis, it actually refers to data lost during the experiment or data that could no longer be tracked. For example, in a clinical trial of a tumor treatment drug, the focus is the event from taking the drug to death, but some patients may not have a death event observed during the trial—such as being unreachable, voluntarily withdrawing, requiring urgent withdrawal from the trial, or still alive when the trial ends. Such data are called censored data. Considering this type of censoring keeps the progress from the beginning until the censoring time but loses subsequent outcomes, it is referred to as right censoring. Conversely, there is also left censoring—for example, when analyzing survival time from initial disease onset to final death, some patients are known to have the disease and death time but the exact initial onset time is unknown. This type of censoring is called left censoring. There are methods to handle left censoring, which we do not discuss here (may be added later if time permits). This article mainly considers right censoring.
Before starting, let’s clarify two concepts. 1. Survival probability, denoted as Survival probability, refers to the probability that the research subject survives from the start of the study until a certain time point t. It is a function of time t, defined as S(t); 2. Hazard probability, denoted as Hazard probability, refers to the probability that the subject has survived up to but experiences the event (e.g., death) exactly at time t. It is also a function of time t, defined as H(t). The upcoming Kaplan-Meier method focuses mainly on S(t), while the Cox proportional hazards model focuses on H(t).
Now let’s get to the main topic and introduce the Kaplan–Meier method, which was jointly proposed by Kaplan and Meier in 1958. To understand the method’s details, first look at the following table (original link). In this example, death is the observed event, and the table is also called a life table.
Column A is the continuous observation time from the trial start; an asterisk (*) indicates censored data occurred at that time. Column B shows the number of subjects at risk (alive and under observation) just before the time point in column A, representing the effective population at risk after excluding deceased and censored cases. Column C shows the number of deaths exactly at the time in column A, and column D shows the number of censored cases at that time. The first row can be read as: before 0.909 years, there were 10 patients; at 0.909 years (or shortly after), one died and no censored data occurred, meaning 9 remained. Subsequently, whenever new deaths or censored data occur, a new row is added to record the time and numbers.
Let’s not yet introduce the Kaplan–Meier formula. You may first try to think about how to calculate survival probability S(t) at each time point. For example, what is the survival probability at 1.536 years, the fifth row in the table? It is easy to think that to survive by 1.536 years, the patient must survive all preceding time points 0.909, 1.112, 1.322, and 1.328. Thus, the survival probability P at 1.536 years actually equals the product of the probabilities of not dying at all these time points including 1.536:
P(survive to 1.536) = P(not dying at 0.909) * P(not dying at 1.112) * P(not dying at 1.322) * P(not dying at 1.328) * P(not dying at 1.536)
The probability of not dying at a particular time point can be estimated as 1 minus the death probability. For example:
P(not dying at 0.909) = 1 – P(dying at 0.909) = 1 – (number of deaths at 0.909) / (number at risk before 0.909) = 1 – 1/10 = 0.9
P(not dying at 1.112) = 1 – P(dying at 1.112) = 1 – (number of deaths at 1.112) / (number at risk before 1.112) = 1 – 1/9 = 0.89
Note that when censored data occur, no deaths happen at that time point, so that time point does not contribute to the cumulative survival probability. However, the total number at risk decreases, affecting the survival probabilities at following time points. For example:
P(not dying at 1.322) = 1 – P(dying at 1.322) = 1 – (number of deaths at 1.322) / (number at risk before 1.322) = 1 – 0/9 = 1
After calculating the non-death probabilities at each time point, you can multiply them sequentially to get the survival probability at each time point, i.e., the probability of surviving up to that time. The table below shows this:
Column E shows the probabilities of not dying at the respective time points, and column F shows the cumulative survival probabilities. As time increases, deaths accumulate and survival probability decreases, which is consistent with common reasoning. Also note that survival probability does not change at censored time points.
Actually, the above thinking is the main idea of the Kaplan–Meier method. Using the table, we can also express it mathematically. Suppose there are m time points, each indexed by i from 1 to m, then the survival probability S(t_i) can be represented as:
where t_i is the ith time point, n_i is the number at risk just before t_i, d_i is the number of deaths at t_i, and S(t_{i-1}) is the survival probability at the previous time point i−1.
Based on this formula, we can plot the change of survival rate as the Kaplan-Meier survival curve, shown below:
The horizontal axis is time, and the vertical axis is the cumulative survival proportion (survival probability). Plus signs indicate censored data. Generally, survival analysis compares survival between groups, so the Kaplan-Meier survival curve usually contains multiple curves, as shown below:
Different colors represent two groups of patients, showing different survival statuses over time. In the figure, the red and blue lines mostly overlap, with the red line a bit higher in the later stage, meaning the red group has a slightly higher late survival probability and slower death rate. If a group’s survival is very poor, its survival curve would descend sharply like a staircase. However, in statistics, quantitative metrics are needed to assess differences and significance rather than relying solely on visual observation.
The most direct metric is the median survival time — the survival time corresponding to a 50% survival rate, as shown here:
Different groups have different median survival times, reflecting differences in mortality risk to some extent. To compare overall survival time distributions for statistical differences, the Logrank test is generally used for survival data.
The Logrank test was first proposed by Nathan Mantel. It is a nonparametric test (Chinese translation: 对数秩检验), mainly used to compare differences in survival distributions between two samples. There are different variants of the Logrank test calculation. Here, we introduce one version, which is essentially an application of the chi-square test. The null hypothesis states that the two groups have identical survival time distributions. If the null hypothesis is rejected, we conclude there exists a statistically significant difference in survival times between groups. The expected number of deaths in a group at time t can be calculated by:
E1t = N1t * (Ot / Nt)
where E1t is the expected number of deaths at time t in group 1; N1t is the number at risk at time t in group 1 (number alive before t); Ot is the total observed deaths at time t in both groups; Nt is the total number at risk in both groups before time t. This formula reflects the logic that if survival distributions are identical, then the proportion of deaths to number at risk in the combined groups should be consistent with that in the single group. Using the total risk population and total deaths proportion from both groups, multiplying by the single group’s risk population gives the expected deaths for that group.
After obtaining expected deaths at each time point, the chi-square statistic is constructed as:
In this formula, the outer Σ sums over different groups, and the three inner Σ are sums over different time points t; j denotes group j (e.g., treatment group and control group corresponding to j=1 and j=2). ΣO_{jt} is the sum of observed deaths across all time points in group j; ΣE_{jt} is the sum of expected deaths across all time points in group j. The difference between observed and expected deaths indicates whether the actual situation matches the hypothesis: if the null hypothesis is true (survival distributions identical), differences should be small. Since differences can be positive or negative, squaring prevents cancellation. Moreover, because an absolute difference of 20 is more significant between 100 and 120 (20% difference) than between 1000 and 1020 (2%), the difference is normalized by dividing by ΣE_{jt}, converting it into a relative measure. Summing these across groups gives the X² value, which can be checked in tables to determine whether to reject the null hypothesis. A manual calculation process can be referenced here: link.
Besides the Logrank test, another common test for survival time distribution is the Breslow test (also known as the Wilcoxon test). Unlike Logrank, Breslow multiplies observed and expected death counts at each time point by a weight factor equal to the total number at risk at that time, then sums them over all time points to compute the chi-square statistic. Intuitively, as time progresses, the number at risk decreases, so the weight diminishes, meaning Breslow is more sensitive to differences in early trial stages, while Logrank, with weight 1 at all times, is more sensitive to later-stage differences. In practice, different methods can be used to explore data from multiple perspectives.
At this point, most knowledge about Kaplan–Meier methods has been covered. But another aspect may catch readers’ attention, as shown below:
You may ask: what are the light-colored areas around the survival curve? They actually represent the 95% confidence intervals corresponding to the survival curve. First, understand why there are confidence intervals for survival probabilities at each node: we assume survival probability at each time point follows a specific distribution, and each observation at a time point is akin to a random sampling. Given that S(t) as the survival probability function of time t is known, by estimating the standard error, we can estimate a Margin of Error and thus obtain a confidence interval with a certain probability.
One formula for estimating confidence intervals for survival probability data was proposed by Greenwood:
where S(t) is the survival probability function (the cumulative product of non-death probabilities explained earlier). Z_{α/2} is the α/2 quantile of the normal distribution, available from standard tables.
Besides the above confidence interval calculation, there is an exponential Greenwood formula that solves the problem where the calculated intervals may not fall within (0,1). Both formula derivations and details can be found here: link.
So far, we have covered the basics of survival analysis, including its applications, censored data explanation, survival and hazard probabilities, Kaplan–Meier curve, Logrank and Breslow tests, and confidence interval estimation. However, the explanations so far concern only univariate analysis, i.e., Kaplan-Meier and Logrank analyze a single variable—such as grouping by gender or drug type—but cannot simultaneously consider multiple factors or adjust for covariates that may affect survival. For example, to study prognosis differences of Covid-19 inside versus outside Hubei Province, imagine we suspect population differences (e.g., more elderly or children in Hubei). We cannot simply compare cases from Hubei and other provinces directly but must adjust for age and gender. This kind of adjustment cannot be done with Logrank, and here the Cox proportional hazards regression model shines.
Cox Model
Now let us explain the Cox proportional hazards regression model, formally named Proportional Hazards Regression analysis or Cox Proportional-Hazards Model, introduced by Cox in 1972. Here, we abbreviate it as the Cox model.
The Cox model is a semi-parametric model because its formula includes both parametric and nonparametric components. Simply put, parametric models have parameters of finite dimension (e.g., mean and standard deviation in normal distribution), while nonparametric models have infinitely dimensional parameters and cannot be described by a finite number of parameters. Different data sets lead to different distribution estimates in nonparametric models (e.g., decision trees, random forests). Looking at the Cox model formula clarifies why it is semi-parametric.
Here, t is survival time; x₁, x₂, …, x_p represent multiple predictor variables; b₁, b₂, …, b_p are coefficients corresponding to each variable’s effect size. h(t) is the hazard function at time t, representing risk; h₀(t) is the baseline hazard function—risk when all covariates x₁, x₂, …, x_p equal zero (have no effect). The exponential part is parametric (finite parameters b₁ to b_p), while the baseline hazard function h₀(t) is nonparametric (unspecified form, depending on data). So, Cox model is semi-parametric.
A key point is that not all survival analysis data suit Cox model analysis; certain assumptions must be met. Keep this question in mind for later discussion at the article’s end.
Now we understand the hazard function h(t) formula. Let’s interpret it. h(t) varies over time t; for a patient, risk dynamics change with time progression (e.g., tumor patients have increasing relapse and death rates over time). From the previously shown Kaplan–Meier table, survival probability decreases at later times, indicating higher risks. Moreover, different patients have different risk based on covariates such as age, gender, blood pressure, etc. For Covid-19, older age corresponds to higher fatality rates. Hence, Cox model incorporates multiple covariates (x₁, x₂, …, x_p) representing potential survival influences. The main goal is to estimate h₀(t) and coefficients b₁, b₂, …, b_p. The Cox model requires methods like maximum likelihood estimation and optimization algorithms (e.g., gradient descent) to find parameter estimates that maximize the likelihood function. Details are omitted here; interested readers can refer to: link 1, link 2, link 3.Assuming we have already calculated the appropriate ( h_0(t) ) and covariate coefficients, how do we interpret the results? We can compare the hazard ratios corresponding to different values of a covariate ( x_1 ). Here, we compare ( x_1 ) and ( x_1 + 1 ), that is, if ( x_1 ) increases by 1 unit, the hazard ratio before and after the increase is:
In the above formula, we calculate the hazard ratios for the two different values ( x_1 + 1 ) and ( x_1 ). By simplification, the hazard ratio corresponding to ( x_1 + 1 ) and ( x_1 ) is actually equal to ( \exp(b_1) ), which is ( e ) raised to the power of ( b_1 ). Simply put, if ( x_1 ) represents age, then for people aged 51 years (( x + 1 )) and 50 years (( x )), the hazard ratio of death is ( \exp(b_1) ). If ( b_1 > 0 ), then ( \exp(b_1) > 1 ), which means that an increase of 1 year in age increases the risk of death; if ( b_1 < 0 ), then ( \exp(b_1) < 1 ), indicating an increase in age reduces the risk of death; if ( b_1 = 0 ), then ( \exp(b_1) = 1 ), meaning age change has no effect on the risk of death.
Therefore, we know that for each covariate, if its coefficient is positive, it indicates that an increase in the corresponding variable increases the patient’s risk of death (or other events, such as relapse or metastasis); if its coefficient is negative, it indicates that an increase in the variable reduces the patient’s risk of death. A specific computational example can be found at link, and an example using R for calculation can be found at: link.
The following table shows a calculation result:
We can see from the table that some risk factors include age, gender, blood pressure (systolic), smoking status, serum total cholesterol, and diabetes status. After calculation with the Cox model, parameter estimates for each risk factor were obtained. For example, the parameter corresponding to age is 0.11691, which is the coefficient mentioned in the previous formula. Since it is greater than 0, it indicates that increased age increases the risk. The hazard ratio is ( \exp(0.11691) = 1.124 ), which is the number in the last column of the table. This value greater than 1 also indicates that an increase in age leads to an increased risk. For binary variables, which only take values 0 and 1 (e.g., male as 1, female as 0), the interpretation of results in the Cox model is consistent with continuous variables. If the covariate coefficient associated with gender is greater than 0, it indicates that the higher value of the variable corresponds to higher risk, meaning males have a higher risk than females. Besides focusing on coefficients, attention must also be paid to the p-value, which tests whether the parameter estimate is statistically significant. Commonly used statistical methods include the Likelihood ratio test, as well as the Wald test and score logrank statistics. A brief introduction to the Likelihood ratio test (also called likelihood ratio test) is: to determine if the introduction of a new variable is effective for the model, compare the maximum value of the likelihood function before and after adding the variable; if there is no decrease in the maximum value, it may be effective for the model, and its significance is then tested. Those interested in these three statistical methods can see link.
As mentioned before, Cox models require certain assumptions to be satisfied. At this point, you should already have an answer: the most important assumption is that the hazard ratio between any two individuals is constant over time, which is why the full name of the Cox model is the Cox proportional hazards regression model. For example, if one person’s risk of death is twice that of another person, then at any time, this person’s risk of death is always twice the other person’s. This assumption can be seen from the previous formulas, since the derived hazard ratio does not include time ( t ). This is a basic assumption that allows the Cox model to be valid. The actual situation may not always conform to this, right? Therefore, before applying the Cox model, it is best to analyze the data to see whether this assumption holds. Since the hazard ratio between individuals is fixed, the simplest example is that under any grouping, the Kaplan–Meier survival curves of the two groups should not cross. If the curves cross, it means the survival probability relationship between the two groups changes with the event, i.e., the hazard ratio changes over time, contradicting the assumption. However, in practical applications, due to smaller sample sizes, survival curves may introduce larger errors, so this judgment method may sometimes fail. A more sophisticated method is the complementary log-log plot, where the horizontal axis is ( \ln(t) ), and the vertical axis is ( \ln(-\ln(S(t))) ). After derivation:
We know that in different groups, the curves of ( \log(-\log(S(t))) ) as a function of ( t ) differ only by ( \beta X ), which does not depend on time ( t ), so the two curves should be approximately parallel or equally spaced. Reference: (Several Graphical Methods for Checking the Proportional Hazards Assumption of the Cox Model_ Yu Hongmei, link).
Besides plotting survival curves to observe data, we can also check proportional hazards via Schoenfeld residuals. For details, refer to link, as well as examining the interaction between variables and time to assess whether the hazard ratio is fixed (link), which we will not elaborate on here. For data that does not satisfy the fixed hazard ratio assumption, time-dependent covariates or stratified Cox models can be used for calculation (link).
In survival analysis, besides the Cox model, there are also some other parametric models available. Unlike the Cox model, these parametric models usually specify a possible distribution of the hazard function, such as the exponential, Weibull, and Gompertz distributions, and then further estimate the corresponding model parameters. As shown in the figure below, (a) is the exponential distribution, whose hazard ratio is constant and is rarely used in practice; (b) is the Weibull distribution, which can adjust the distribution shape using different parameters; (c) is the Gompertz distribution. Compared to the Cox model, the advantage of using these models is that the distribution curve can be inferred from parameters, providing more information, such as higher early mortality and lower late mortality; in other words, more information about the risk distribution can be obtained, while the Cox model only yields limited information such as hazard ratios and their significance. The disadvantage of these fully parametric models is obvious: the fixed distribution may not fit the actual data situations, potentially causing more errors. In practice, different models can be chosen depending on the situation. A brief introduction to the Weibull model can be found at link.
The above is a rough introduction to the main knowledge in survival analysis. Hopefully, everyone now has an intuitive understanding of most concepts in survival analysis. However, as the author’s capability is limited, there may be unclear or even incorrect statements in this article, and readers are encouraged to provide suggestions and criticisms or leave comments and discussions if there are any questions. Thank you!