15. Survival Analysis

Though the Cutler-Ederer estimate is still in use [14–18], it has been largely replaced as a method for estimation of survival curves by the Kaplan-Meier estimate. Nonetheless it is useful as an introduction to survival curve estimation.

In the Cutler-Ederer or actuarial estimate [12], the assumption is that the deaths and losses are uniformly distributed over a set of fixed-length intervals. On the average, this means that one half the losses will occur during the first half of each interval. The estimate for the probability of surviving the j ^th interval, given that the previous intervals were survived, is , where The λ _j losses are assumed to be at risk, on the average, one half the time and thus should be counted as such. These conditional probabilities are then multiplied together to obtain an estimate, Ŝ(t), of the survival function at time t.

The Kaplan-Meier Estimate relaxes the assumption of events distributed uniformly across fixed length intervals. Using the time of death, observations can be ranked. This is a useful improvement, since in a clinical trial with staggered entry of participants and censored observations, survival data will be of varying degrees of completeness.

As a very simple example, suppose that 100 participants were entered into a study and followed for 2 years. One year after the first group was started, a second group of 100 participants was entered and followed for the remaining year of the trial. Assuming no losses to follow-up, the results might be as shown in Table 15.1. For Group I, 20 participants died during the first year and of the 80 survivors, 20 more died during the second year. For Group II, which was followed for only 1 year, 25 participants died. Now suppose the investigator wants to estimate the 2-year survival rate. The only group of participants followed for 2 years was Group I. One estimate of 2-year survival, S(2), would be or 0.60. Note that the first-year survival experience of Group II is ignored in this estimate. If the investigator wants to estimate 1 year survival rate, S(1), she would observe that a total of 200 participants were followed for at least 1 year. Of those, 155 (80 + 75) survived the first year. Thus, or 0.775. If each group were evaluated separately, the survival rates would be 0.80 and 0.75. In estimating the 1-year survival rate, all the available information was used, but for the 2-year survival rate the 1-year survival experience of Group II was ignored.

Another procedure for estimating survival rates is to use a conditional probability. For this example, the probability of 2-year survival, S(2), is equal to the probability of 1-year survival, S(1), times the probability of surviving the second year, given that the participant survived the first year, pr(2|1). That is, . In this example, . The estimate for pr(2|1) is 60/80 = 0.75 since 60 of the 80 participants who survived the first year also survived the second year. Thus, the estimate for or 0.58, which is slightly different from the previously calculated estimate of 0.60.

Kaplan and Meier [13] described how this conditional probability strategy could be used to estimate survival curves in clinical trials with censored observations. Their procedure is usually referred to as the Kaplan-Meier estimate, or sometimes the product-limit estimate, since the product of conditional probabilities leads to the survival estimate. This procedure assumes that the exact time of entry into the trial is known and that the exact time of the event or loss of follow-up is also known. For some applications, time to the nearest month may be sufficient, while for other applications the nearest day or hour may be necessary. Kaplan and Meier assumed that a death and loss of follow-up would not occur at the same time. If a death and a loss to follow-up are recorded as having occurred at the same time, this tie is broken on the assumption that the death occurred slightly before the loss to follow-up.

In this method, the follow-up period is divided into intervals of time so that no interval contains both deaths and losses. Let p _j be equal to the probability of surviving the j ^th interval, given that the participant has survived the previous interval. For intervals labeled j with deaths only, the estimate for p _j , which is , is equal to the number of participants alive at the beginning of the j ^th interval, n _j , minus those who died during the interval, δ _j , with this difference being divided by the number alive at the beginning of the interval, i.e. . For an interval j with only l _j losses, the estimate is one. Such conditional probabilities for an interval with only losses would not alter the product. This means that an interval with only losses and no deaths may be combined with the previous interval.

Table 15.2 presents the survival experience for this example as a lifetable. Each row in the lifetable indicates the time at which a death or an event occurred. One or more deaths may have occurred at the same time and they are included in the same row in the lifetable. In the interval between two consecutive times of death, losses to follow-up may have occurred. Hence, a row in the table actually represents an interval of time, beginning with the time of a death, up to but not including the time of the next death. In this case, the first interval is defined by the death at 0.5 months up to the time of the next death at 1.5 months. The columns labeled n _j , δ _j , and l _j correspond to the definitions given above and contain the information from the example. In the first interval, all 20 participants were initially at risk, one died at 0.5 months, and later in the interval (at 0.6 months) one participant was lost to follow-up. In the second interval, from 1.5 months up to 3.0 months, 18 participants were still at risk initially, two deaths were recorded at 1.5 months and one participant was lost at 2.0 months. The remaining intervals are defined similarly. The column labeled is the conditional probability of surviving the interval j and is computed as (n _j  − δ _j )/n _j or (20 − 1)/20 = 0.95, (18 − 2)/18 = 0.89, etc. The column labeled Ŝ(t) is the estimated survival curve and is computed as the accumulated product of the (0.85 = 0.95 × 0.89, 0.79 = 0.95 × 0.89 × 0.93, etc).

The graphical display of the next to last column of Table 15.2, Ŝ(t), is given in Fig. 15.6. The step function appearance of the graph is because the estimate of S(t), Ŝ(t) is constant during an interval and changes only at the time of a death. With very large sample sizes and more observed deaths, the step function has smaller steps and looks more like the usually visualized smooth survival curve. If no censoring occurs, this method simplifies to the number of survivors divided by the total number of participants who entered the trial.

Because Ŝ(t) is an estimate of S(t), the true survival curve, the estimate will have some variation due to the sample selected. Greenwood [19] derived a formula for estimating the variance of an estimated survival function which is applicable to the Kaplan-Meier method. The formula for the variance of Ŝ(t), denoted V[Ŝ(t)] is given by where n _j and δ _j are defined as before, and K is the number of intervals. In Table 15.2, the last column labeled V[Ŝ(t)] represents the estimated variances for the estimates of S(t) during the six intervals. Note that the variance increases as one moves down the column. When fewer participants are at risk, the ability to estimate the survival experience is diminished.

Other examples of this procedure, as well as a more detailed discussion of some of the statistical properties of this estimate, are provided by Kaplan and Meier [13]. Computer programs are available [20–23] so that survival curves can be obtained quickly, even for very large sets of data.

The Kaplan-Meier curve can also be used to estimate the hazard rate, λ, if the survival curve is exponential. For example, if the median survival time is estimated as T_M, then 0.5 = S(T_M) = exp(−λT_M) and thus as an estimate of λ. Then the estimate for S(t) would be exp (−t) In comparison to the Kaplan-Meier, another parametric estimate for S(t) at time t _j , described by Nelson [24], is where δ _i is the number of events in the i ^th interval and n _i is the number at risk for the event. While this is a straightforward estimate, the Kaplan-Meier does not assume an underlying exponential distribution and thus is used more than this type of estimator.

One possible comparison between groups is to specify a time t* for which survival estimates have been computed using the Kaplan-Meier [13] method. At time t*, one can compare the survival estimates Ŝ _C (t *) and Ŝ _I (t *) using the statistic where V[Ŝ _C (t *)] and V[Ŝ _I (t *)] are the Greenwood estimates of variance [19]. The statistic Z(t*) has approximately a normal distribution with mean zero and variance one under the null hypothesis that . The problem with this approach is the multiple looks issue described in Chap. 16. Another problem exists in interpretation. For example, what conclusions should be drawn if two survival curves are judged significantly different at time t* but not at any other points? The issue then becomes, what point in the survival curve is most important.

For some studies with a T year follow-up, the T year mortality rates are considered important and should be tested in the manner just suggested. Annual rates might also be considered important and, therefore, compared. One criticism of this suggestion is that the specific points may have been selected post hoc to yield the largest difference based on the observed data. One can easily visualize two survival curves for which significant differences are found at a few points. However, when survival curves are compared, the large differences indicated by these few points are not supported by the overall survival experience. Therefore, point-by-point comparisons are not recommended unless a few points can be justified and specified in the protocol prior to data analysis.

One summary measure of survival experience is the time at which 50% of the cohort has had the event. One common and easy way to estimate the median survival time is from the Kaplan-Meier curve. (See for example, Altman [1].) This assumes that the cohort has been followed long enough so that over one-half of the individuals have had the event. Confidence intervals may be computed for the median survival times [25]. If this is the case, we can compare the median survival times for intervention and control M _I and M _C , respectively. This is most easily done by estimating the ratio of the estimates M _I /M _C . A ratio larger than unity implies that the intervention group has a longer median survival and thus a better survival experience. A ratio less than unity would indicate the opposite.

We can estimate 95% confidence intervals for M _I /M _C by where the standard deviation, SD, of M _I /M _C is computed as for cases where the survival curves are approximately exponential, and O _I  = the total number of events in the intervention group (i.e., ∑δ _i ) and O _C  = the total number of events in the control group.

Because of the limitations of comparison of point-by-point estimates, Gehan [26] and Mantel [27] originally proposed statistical methods to assess the overall survival experience. These two methods were important steps in the development of analytical methods for survival data. They both assume that the hypothesis being tested is whether two survival curves are equal, or whether one is consistently different from the other. If the two survival curves cross, these methods should be interpreted cautiously. Since these two original methods, an enormous literature has developed on comparison of survival curves and is summarized in several texts [8–11]. The basic methods described here provide the fundamental concepts used in survival analysis.

Mantel [27] proposed the use of the procedure described by Cochran [28] and Mantel and Haenszel [29] for combining a series of 2 × 2 tables. In this procedure, each time, t _j , a death occurs in either group, a 2 × 2 table is formed as follows:

The entry a _j represents the observed number of deaths at time t _j in the intervention group and c _j represents the observed number of deaths at time t, in the control group. At least a _j or c _j must be non-zero. One could create a table at other time periods (that is, when a _j and c _j are zero), but this table would not make any contribution to the statistic. Of the n _j participants at risk just prior to time t _j , a _j  + b _j were in the intervention group and c _j  + d _j were in the control group. The expected number of deaths in the intervention group, denoted E(a _j ), can be shown to be and the variance of the observed number of deaths in the intervention group, denoted as V(a _j ) is given by These expressions are the same as those given for combining 2 × 2 tables in the Appendix of Chap. 17. The Mantel-Haenszel (MH) statistic is given by and has approximately a chi-square distribution with one degree of freedom, where K is the number of distinct event times in the combined intervention and control groups. As an asymptotic approximation, the (signed) square root of MH, can be compared to a standard normal distribution [30, 31].

Application of this procedure is straightforward. First, the times of events and losses in both groups are ranked in ascending order. Second, the time of each event, and the total number of participants in each group who were at risk just before the death (a _j  + b _j , c _j  + d _j ) as well as the number of events in each group (a _j , c _j ) are determined. With this information, the appropriate 2 × 2 tables can be formed.

The Mantel-Haenszel statistic can be computed from these eight 2 × 2 tables (Table 15.4) or directly from Table 15.3. The term since there are only two deaths in the intervention group. Evaluation of the term or . The value for is computed as

This term is equal to 2.21. The computed statistic is MH = (2 − 4.89)²/2.21 = 3.78. This is not significant at the 0.05 significance level for a chi-square statistic with one degree of freedom. The MH statistic can also be used when the precise time of death is unknown. If death is known to have occurred within an interval, 2 × 2 tables can be created for each interval and the method applied. For small samples, a continuity correction is sometimes used. The modified numerator is where the vertical bars denote the absolute value. For the example, applying the continuity correction reduces the MH statistic from 3.76 to 2.59.

Gehan [26] developed another procedure for comparing the survival experience of two groups of participants by generalizing the Wilcoxon rank statistic. The Gehan statistic is based on the ranks of the observed survival times. The null hypothesis, S _I (t) = S _C (t), is tested. The procedure, as originally developed, involved a complicated calculation to obtain the variance of the test statistic. Mantel [32] proposed a simpler version of the variance calculation, which is most often used.

The N _I observations from the intervention group and the N _C observations from the control group must be combined into a sequence of N _C  + N _I observations and ranked in ascending order. Each observation is compared to the remaining N _C  + N _I  − 1 observation and given a score U _i which is defined as follows:

The survival outcome for the i ^th participant will certainly be larger than that for participants who died earlier. For censored participants, it cannot be determined whether survival time would have been less or greater than the i ^th observation. This is true whether the i ^th observation is a death or a loss. Thus, the first part of the score U _i assesses how many deaths definitely preceded the i ^th observation. The second part of the U _i score considers whether the current, i ^th , observation is a death or a loss. If it is a death, it definitely precedes all later ranked observations regardless of whether the observations correspond to a death or a loss. If the i ^th observation is a loss, it cannot be determined whether the actual survival time will be less than or greater than any succeeding ranked observation, since there was no opportunity to observe the i ^th participant completely.

Table 15.5 ranks the 40 combined observations (N _C  = 20, N _I  = 20) from the example used in the discussion of the Mantel-Haenszel statistic. The last 19 observations were all censored at 12 months of follow-up, 7 in the control group and 12 in the intervention group. The score U ₁ is equal to the zero observations that were definitely less than 0.5 months, minus the 39 observations that were definitely greater than 0.5 months, or U ₁ = −39. The score U ₂ is equal to the one observation definitely less than the loss at 0.6 months, minus none of the observations that will be definitely greater, since at 0.6 months the observation was a loss, or U ₂ = 1. U ₃ is equal to the one observation (0.5 months) definitely less than 1.0 month minus the 37 observations definitely greater than 1.0 month giving U ₃ = 36. The last 19 observations will have scores of 9 reflecting the nine deaths which definitely precede censored observations at 12.0 months.

The Gehan statistic, G, involves the scores U _i and is defined as where W = Σ U _i , (for U _i ’s in control group only) and The G statistic has approximately a chi-square distribution with one degree of freedom [26, 32]. Therefore, the critical value is 3.84 at the 5% significance level and 6.63 at the 1% level. In the example, W = −87 and the variance V(W) = 2,314.35. Thus, G = (−87)²/2,314.35 = 3.27 for which the p-value is equal to 0.071. This is compared with the p-value of 0.052 obtained using the Mantel-Haenszel statistic.

The Gehan statistic assumes the censoring pattern to be equal in the two groups. Breslow [33] considered the case in which censoring patterns are not equal and used the same statistic G with a modified variance. This modified version should be used if the censoring patterns are radically different in the two groups. Peto and Peto [34] also proposed a version of a censored Wilcoxon test. The concepts are similar to what has been described for Gehan’s approach. However, most software packages now use the Breslow or Peto and Peto versions.

The general methodology of comparing two survival curves using this methodology has been further evaluated [35–40]. These two tests by Mantel-Haenzel and Gehan, can be viewed as a weighted sum of the difference between observed number of events and the expected number at each unique event time [7, 40]. Consider the previous equation for the logrank test and rewrite the numerator as where

and w _j is a weighting factor. The test statistic W ² /V(W) has approximately a chi-square distribution with one degree of freedom or equivalently has approximately a standard normal distribution. If w _i  = 1, we obtain the Mantel-Haenszel or logrank test. If w _i  = n _j /(N + 1), where N = N _C  + N _I or the combined sample size, we obtain the Gehan version of the Wilcoxon test. Tarone and Ware [40] pointed out that the Mantel-Haenszel and Gehan are only two possible statistical tests. They suggested a general weight function w _i  = [n _j /(N + 1)] ^θ where 0 < θ < 1. In particular, they suggested that θ = 0.5. Prentice [38] suggested a weight where d _i  = (a _i  + c _i ) which is related to the product limit estimator at t _j as suggested by Peto and Peto [34]. Harrington and Fleming [35] generalize this further by suggesting weights for ρ > 0.

All of these methods give different weights to the various parts of the survival curve. The Mantel-Haenszel or logrank statistic is more powerful for survival distributions of the exponential form where λ _I (t) = θ λ _C (t) or S _I (t) = {S _C (t)} ^θ where θ ≠ 1 [32]. The Gehan type statistic [26], on the other hand, is more powerful for survival distributions of the logistic form S(t,θ) = e ^t + θ /(1 + e ^t + θ ). In actual practice, however, the distribution of the survival curve of the study population is not known. When the null hypothesis is not true, the Gehan type statistic gives more weight to the early survival experience, whereas the Mantel-Haenszel weights the later experience more. Tarone and Ware [40] indicate other possible weighting schemes could be proposed which are intermediate to these two statistics [35, 40]. Thus, when survival analysis is done, it is certainly possible to obtain different results using different weighting schemes depending on where the survival curves separate, if they indeed do so. The logrank test is the standard in many fields such as cancer and heart disease. The condition λ _I (t) = θ λ _C (t) says that risk of the event being studied in the intervention is a constant multiple of the hazard λ _C (t). That is, the hazard rate in one arm is proportional to the other and so the logrank test is best for testing proportional hazards. This idea is appealing and is approximately true for many studies.

There has been considerable interest in asymptotic (large sample) properties of rank tests [37, 39] as well as comparisons of the various analytic methods [36]. While there exists an enormous literature on survival analysis, the basic concepts of rank tests can still be appreciated by the methods described above.

Earlier, we discussed using an exponential model to summarize a survival curve where the hazard rate λ determines the survival curve. If we can assume that the hazard rate is reasonably constant during the period of follow-up for the intervention and the control group, then comparison of hazard rates is a comparison of survival curves [1]. The most commonly used comparison is the ratio of the hazards, R = λ _I /λ _C . If the ratio is unity, the survival curves are identical. If R is greater than one, the intervention hazard is greater than control so the intervention survival curve falls below the standard curve. That is, the intervention is worse. On the other hand, if R is less than one, the control group hazard is larger, the control group survival curve falls below the intervention curve, and intervention is better.

We can estimate the hazard ratio by comparing the ratio of total observed events (O) divided by expected number of events (E) in each group; that is, the estimate of R can be expressed as That is, O _I = Σa _i , O _C = Σ b _i , E _I = ΣE(a _i ), and E _C =ΣE(b _i ). Confidence intervals for the odds ratio R are most easily determined by constructing confidence intervals for the log of the odds ratio ln R [41]. The 95% confidence interval for ln R is to where K = (O _I  − E _I )/V and V is the variance as defined in the logrank or Mantel-Haenszel statistics. (That is, V equals V(a _i ).) We then connect confidence intervals for ln R to confidence intervals for R by taking antilogs of the upper and lower limit. If the confidence interval excludes unity, we could claim superiority of either intervention or control depending on the direction. Hazard ratios not included in the interval can be excluded as likely outcome summaries of the intervention. If the survival curves have relatively constant hazard rates, this method provides a nice summary and complements the Kaplan-Meier estimates of the survival curves.

Previous chapters have discussed the rationale for taking stratification into account. If differences in important covariates or prognostic variables exist at entry between the intervention and control groups, an investigator might be concerned that the analysis of the survival experience is influenced by that difference. In order to adjust for these differences in prognostic variables, she could conduct a stratified analysis or a covariance type of survival analysis. If these differences are not important in the analysis, the adjusted analysis will give approximately the same results as the unadjusted.

Three basic techniques for stratified survival analysis are of interest. The first compares the survival experience between the study groups within each stratum, using the methods described in the previous section. By comparing the results from each stratum, the investigator can get some indication of the consistency of results across strata and the possible interaction between strata and intervention.

The second and third methods are basically adaptations of the Mantel-Haenszel and Gehan statistics, respectively, and allow the results to be accumulated over the strata. The Mantel-Haenszel stratified analysis involves dividing the population into S strata and within each stratum j, forming a series of 2 × 2 tables for each K _j event, where K _j is the number of events in stratum j. The table for the i ^th event in the j ^th stratum would be as follows:

The entries a _ij , b _ij , c _ij , and d _ij are defined as before and Similar to the non-stratified case, the Mantel-Haenszel statistic is which has a chi-square distribution with one degree of freedom. Analogous to the Mantel-Haenszel statistic for stratified analysis, one could compute a Gehan statistic W _j and V(W _j ) within each stratum. Then an overall stratified Gehan statistic is computed as which also has chi-square statistic with one degree of freedom.

If there are many covariates, each with several levels, the number of strata can quickly become large, with few participants in each. Moreover, if a covariate is continuous, it must be divided into intervals and each interval assigned a score or rank before it can be used in a stratified analysis. Cox [42] proposed a regression model which allows for analysis of censored survival data adjusting for continuous as well as discrete covariates, thus avoiding these two problems.

One way to understand the Cox regression model is to again consider a simpler parametric model. If one expresses the probability of survival to time t, denoted S(t), as an exponential model, then S(t) = e ^−λt where the parameter, λ, is called the force of mortality or the hazard rate as described earlier. The larger the value of λ, the faster the survival curve decreases. Some models allow the hazard rate to change with time, that is λ = λ(t). Models have been proposed [43–45] which attempt to incorporate the hazard rate as a linear function of several baseline covariates, x ₁, x ₂, …, x _p that is, λ(x ₁, x ₂, …, x _p ) = b ₁ x ₁ + b ₂ x ₂ + … + b _p x _p . One of the covariates, say x ₁, might represent the intervention and the others, for example, might represent age, sex, performance status, or prior medical history. The coefficient, b ₁, then would indicate whether intervention is a significant prognostic factor, i.e., remains effective after adjustment for the other factors. Cox [42] suggested that the hazard rate could be modeled as a function of both time and covariates, denoted λ(t, x ₁, x ₂, …, x _p ). Moreover, this hazard rate could be represented as the product of two terms, the first representing an unadjusted force of mortality λ ₀(t) and the second the adjustment for the linear combination of a particular covariate profile. More specifically, the Cox proportional hazard model assumes that That is, the hazard λ(t, x ₁, x ₂, …, x _n ) is proportional to an underlying hazard function λ ₀(t) by the specific factor exp (b ₁ x ₁ + b₂ x ₂ …). From this model, we can estimate an underlying survival curve S ₀(t) as a function of λ ₀(t). The survival curve for participants with a particular set of covariates X, S(t,x) can be obtained as . Other summary test statistics from this model are also used. The estimation of the regression coefficients b ₁,b ₂, …, b _p is complex, requiring non-linear numerical methods, and goes beyond the scope of this text. Many elementary texts on biostatistics [1, 3, 5, 46] or review articles [7] present further details. A more advanced discussion may be found in Kalbfleish and Prentice [10] or Fleming and Harrington [9]. Programs exist in many statistical computing packages which provide these estimates and summary statistics to evaluate survival curve comparisons [20–23]. Despite the complexity of the parameter estimation, this method is widely applied and has been studied extensively [47–55]. Pocock, Gore, and Kerr [52] demonstrate the value of some of these methods with cancer data. For the special case where group assignment is the only covariate, the Cox model is essentially equivalent to the Mantel-Haenszel statistic.

One issue that is sometimes raised is whether the hazard rates are proportional over time. Methods such as the Mantel-Haenszel logrank test or the Cox Proportional Hazards model are optimal when the hazards are proportional [9]. However, though there is some loss of power, these methods perform well as long as the hazard curves do not cross, even if proportionality does not hold [56]. When the hazards are not proportional, which intervention is better depends on what time point is being referenced. If a significant difference is found between two survival curves using the Mantel-Haenszel logrank test or the Cox Proportional Hazards model when the hazards are not proportional, the two curves are still significantly different. For example, time to event curves are shown in Chap. 18. Figure 18.​2a shows three curves for comparison of two medical devices with best medical or pharmacologic care. These three curves do not have proportional hazards but the comparisons are still valid and in fact the two devices demonstrate statistically significant superiority over the best medical care arm. The survival curves do not cross although are close together in the early months of follow-up.

The techniques described in this chapter as well as the extensions or generalizations referenced are powerful tools in the analysis of survival data. Perhaps none is exactly correct for any given set of data but experience indicates they are fairly robust and quite useful.