STATISTICS IN MEDICINE. VOL. 10, 733-738 (1991)
A NOTE ON THE CALCULATION OF EXPECTED SURVIVAL, ILLUSTRATED BY THE SURVIVAL OF LIVER TRANSPLANT PATIENTS BIRTHE LYKKE THOMSEN AND NIELS KEIDING Statistical Research Unit, Faculty of Medicine, University of Copenhagen, Blegdamsvej 3, DK-2200 Copenhagen N . , Denmark
AND DOUGLAS G. ALTMAN Medical Statistics Laboratory. Imperial Cancer Research Fund, P.O. Box 123, Lincoln's Inn Fieldr, London WC2A 3PX. U.K.
SUMMARY Based on historical regression analyses of survival one may calculate expected survival curves for patients in a study group, based on their observed covariates. This note discusses this calculation and illustrates it on a comparison of the survival of a certain group of liver transplant patients to their prognosis under conservative treatment.
INTRODUCTION With the increased availability of published regression analyses of survival with specific diseases it has become possible to make historical comparisons between the observed survival in a (possibly fairly limited) study group and the survival that was to be expected from the published regression analyses. One area where this procedure has been used is transplantation, where randomized studies are ethically impossible, so that the best approximation to an answer to the question 'what would have happened to these patients if conservative treatment had been preferred? will be via comparison with historical controls. We discuss below some issues of the exact choice of adjustment procedure, in the context of a study on liver transplantation.'
THE DIRECT ADJUSTED SURVIVAL CURVE Neuberger et al.' reported on the survival of the first 29 patients with primary biliary cirrhosis (PBC) who had a liver transplant in the Cambridge/King's College Hospital programme. Their Kaplan-Meier survival curve is shown in Figure 1. The Kaplan-Meier curve may be regarded as an estimate of the mean survival function S,(t) = ZSi(t)/29,where Si(t)is the survival function of patient i. In the present situation an estimate was available of Si(t) under the hypothetical assumption that patient i had continued on conservative treatment rather than having been transplanted. 0277-671 5/9 1/05073346$05.00 0 1991 by John Wiley & Sons, Ltd.
Received March 1989 Revised August 1990
734
B. L. THOMSEN, N. KEIDING AND D. G. ALTMAN
i
12
21
60
months Figure 1. Observed survival’ (- . -)of 29 PBC patients who had a liver transplantation,compared with the distribution expected from an earlier clinical trial of PBC patients,’ calculatedas direct adjusted survival curve (+-) and accordingto Neuberger et al.’ (- - -)
A recent Cox regression analysis of a randomized clinical trial2 of PBC patients allowed the calculation of the estimated survival probability for time t for patient i with covariates zi as fi(t)= f(t,zi). Accordingly, the obvious estimate of the mean survival function S,(t) of the 29 transplant patients with covariate vectors zl,. . . , zZ9 is given as the mean of the individual estimated survival curves
c
q t ) = 291 291 f(t,q), -
i=
and this was termed the ‘direct adjusted survival curve’ by Makuch3 and Gail and B ~ a rMarkus ;~ et al.’ used it in the liver transplantation context. Neuberger et al. used a somewhat different adjustment procedure, averaging log(-log f(t, zf)) rather than g(t, zi). The direct adjusted survival curve and the approach by Neuberger et al. were both very recently discussed by Bonsel et ~ l . , ~ who also suggested a third method, albeit not explained in their paper, which ostensibly takes the staggered entry directly into account. Figure 1 shows the adjusted survival curves; Neuberger et al.’s suggestion represents a survival curve with less spread than the direct adjusted survival curve. In the Cox model, with standard notation, the hazard function of patient i is 4(t) = &(t)exP(Pzi)-
Thus log( - logf(t, Zi)) =
pz, + logA,(t);
where A,(t) is Breslow’s’ estimator of the cumulative hazard A,@) = j:=, I,(u)du, and bzi can be considered to be a prognostic index for a patient with covariate vector zi. The averaging suggested by Neuberger et al. produces an adjusted survival curve of S,(t) = exp(- A,(t)e@f)
which is also the predicted survival curve of a patient with prognostic index Pz equal to the average prognostic index of the patients, in particular of a patient with all covariates equal to the
A NOTE ON THE CALCULATION OF EXPECTED SURVIVAL
735
a
0.8-
r ._ D
-
,
0.6-
0 L
.
a
2"
.-
O.'-
2 , ' 7 rA
0.2-
months
b
0.0 -
4
12
21
60
months
Figure 2. Survival curves for two 'samples' consisting of the two most extreme patients (- - -) and of two 'middle' patients (-), respectively, (a) direct adjusted, and (b) adjusted according to Neuberger et al.'
average Z of the covariate values. Though they did not derive their adjustment this way, they thus arrived at the same result as if they had disregarded the sample variation. As elaborated below, this explains the above noticed difference in spread of the two adjustments. At first sight there may seem to be little difference in principle between the two averaging operations suggested above and it may require some adaptation to appreciate that s(t)accounts for the sample variability whereas SN(t)does not. To make the point clearer, we consider two 'samples', each consisting of two patients. Sample A contains the two patients with the lowest (5.12) and highest (8.36)prognostic indices, having an average of 6-74. Sample B contains the two patients with prognostic indices closest to this average, 6.66 and 6.67. Figure 2(a) shows that the direct adjusted survival curves qt)of samples A and B are very different, reflecting the large difference in sample variation between samples A and B. In contrast, Figure 2(b) shows that averaging the prognostic indices in the Neuberger et al. fashion, the adjusted survival curves &(t) are (almost) identical, depending only on the (almost) identical average prognostic indices. The difference between the survival curves estimated by the two methods depends on the variability in
736
B. L. THOMSEN, N. KEIDING AND D. G. ALTMAN
the sample of the prognostic indices. In the present example, only PBC patients with a bad prognosis were offered a liver transplant, implying a smaller variability of the prognostic indices than often encountered in Cox analyses. Recent extensive investigations in population heterogeneitysv9show that the survival distribution of the population of patients cannot have the same time-dependence of the hazard as that of the individual patients, because, briefly put, the frail die first, leading to relative reduction of population mortality over time as compared to the pattern of any individual patient. The intensity L,(t) epi of S,(t) is proportional to the underlying intensity L,(t) and thus to any individual death intensity; in contrast $(t) has intensity
where the second factor is decreasing in t, thus capturing the just mentioned aspect of the heterogeneous sample. Besides the prediction variation of an individual’s survival time around its expectation and the sample variation among the new individuals as expressed by the variation among their covariate values, a third source of variation is important: the estimation variation of the parameters of the model. It is, however, beyond the scope of this note to discuss derivation of a numerical test - or, equivalently, confidence intervals - of the difference between the observed and the direct adjusted survival curve. THE ‘EXPECTED SURVIVAL CURVE A different concept of the survival curve to be expected for a study population if individuals are assumed to be exposed to some known ‘standard’ death rates has been extensively used in connection with the grouped-time ‘actuarial’ approach to survival and was recently put into the context of continuous-time survival ana1y~is.l~ We briefly indicate the latter version in the adaptation relevant here. If the mortality of the patients followed that of the PBC patients in the clinical trial,’ the average hazard of the patients still under observation at time t would be L*(t) =
Yi(t)Lo(t)exP(B’zi)/Y(t)
where Yi(t)indicates that patient i is still at risk at time t and Y(t)= XYi(t). Defining the cumulative hazard as A*(t) = f t = o L*(u)du, the survival function S*(t) = exp{ - A*(t)]
is a continuous-time version of the ‘expected survival rate’” or ‘expected survival curve’.’’ It has a desirable property:I3 let the standard Nelson-Aalen estimator of A(t) be defined by
where T , < T2 < , . . are the times of (observed) deaths. Then A*(t) - A(t) has expectation zero, under the null hypothesis that patient i has hazard L,(t). Therefore S*(t) represents, under the null hypothesis, an expected survival curve. Outside of the null hypothesis it is not so easy to interpret S*(t), containing as it does information on study group mortality, through Yi(t),as well as
737
A NOTE ON THE CALCULATION OF EXPECTED SURVIVAL
0
4
12
28
60
months Figure 3. Expected survival according to Neuberger et al.' (- - -), as direct adjusted survival curve (-) survival curve'" (- . . -)
and as 'expected
standard mortality, through l , ( t ) .S*(t) therefore cannot be recommended as an expected survival function, despite its wide use for this purpose. A similar difficulty outside of the null hypothesis was recently discussed for the well-known person-years method. l4 S*(t) is, however, useful for inference on excess mortality, as follows. In the particular case where the study group has an additive excess mortality a(t) over the standard, that is, patient i has hazard a(t) + ,Ii([), one may obtain an estimate of the corresponding 'survival' function
as the so-called relative survival function s(t)/S*(t).Andersen and Vaeth' showed that results on unbiasedness, consistency and asymptotic normality are available. In the present example S*(t) is even steeper than S,(t) (see Figure 3). Our explantion is that mortality after transplantation depends on different risk factors from mortality under medical treatment, so that the order in which patients actually die is different from what it 'should be' according to l i ( t ) .This also implies that the simple additive model for the hazard is inadequate in the present example. Whatever the explanation, the case documents the drawback of using the 'expected survival rate' S * ( t ) rather than the immediately interpretable 'direct adjusted survival curve' Q t ) . Our conclusion is that for the specific purpose of calculating 'what the survival would have been if some given individual-specific mortality rates were operating', the direct adjusted survival curve should be preferred. ACKNOWLEDGEMENTS
This study was supported by the Danish Natural Science Council grant no. 11-7575. REFERENCES
1. Neuberger, J., Altman, D. G., Christensen, E., Tygstrup, N. and Williams, R. 'Use of a prognostic index in evaluation of liver transplantation for primary biliary cirrhosis', Transplantation, 4, 71 3-716 (1986).
738
B. L. THOMSEN, N. KEIDING AND D. G . ALTMAN
2. Christensen, E., Neuberger, J., Crowe, J., Altman, D. G., Popper, H., Portmann, B., Doniach, D..
Ranek, L., Tygstrup, N. and Williams, R. ‘Beneficial effect of azathioprine and prediction of prognosis in primary biliary cirrhosis. Final results of an international trial’, Gastroenterology, 89, 1084-1091 (1985). 3. Makuch, R. W. ‘Adjusted survival curve estimation using covariates’, Journal of Chronic Diseases, 3, 437-443 (1982). 4. Gail, M. H. and Byar, D. B. ‘Variance calculations for direct adjusted survival curves with applications to testing for no treatment effect’, Biometrical Journal, 28, 587-599 (1986). 5. Markus, B. H., Dickson, E. R., Grambsch, P. M., Fleming, T. R., Mazzaferro, V., Klintmalm, G. B. G.,
Wiesner, R. H., van Thiel, D. H. and Starzl, T. E. ‘Efficacy of liver transplantation in patients with primary biliary cirrhosis’, The New England Journal of Medicine, 320, 1709-1713 (1989). 6. Bonsel, G. J., Klompmaker, I. J., van’t Veer, F., Habbema, J. D. F. and Slooff, M. J. H. ‘Use of prognostic models for assessment of value of liver transplantation in primary biliary cirrhosis’, “he Lancet, 335, 493-497 (1990). 7. Breslow, N. Contribution to the discussion of D. R. Cox ‘Regression models and life tables’, Journal of the Royal Statistical Society, Series B, 34, 216217 (1972). 8. Hougaard, P. ‘Life table methods for heterogeneous populations: Distributions describing the heterogeneity’, Biometrika, 71, 75-83 (1984). 9. Aalen, 0.0. ‘Heterogeneity in survival analysis’, Statistics in Medicine, 7 , 1121-1137 (1988). 10. Ederer, F., Axtell, L. M. and Cutler, S. J. ‘The relative survival rate: A statistical methodology’, National Cancer Institute Monographs, 6, 101-121 (1961). 11. Lee, E. T. Statistical Models for Survival Data Analysis, Lifetime Learning Publications, Belmont, California, 1980, Chapter 4. 12. Hill, C., Laplanche, A. and Rezvani, A. ‘Comparison of the mortality of a cohort with the mortality of a reference population in a prognostic study’, Statistics in Medicine, 4, 295-302 (1985). 13. Andersen, P. K. and Vaeth, M. ‘Simple parametric and nonparametric models for excess and relative mortality’, Biometrics, 45, 523-535 (1989). 14. Keiding, N. and Vaeth, M. ‘Calculating expected mortality’, Statistics in Medicine, 5, 327-334 (1986).
A NOTE ON THE CALCULATION OF EXPECTED SURVIVAL, ILLUSTRATED BY THE SURVIVAL OF LIVER TRANSPLANT PATIENTS BIRTHE LYKKE THOMSEN AND NIELS KEIDING Statistical Research Unit, Faculty of Medicine, University of Copenhagen, Blegdamsvej 3, DK-2200 Copenhagen N . , Denmark
AND DOUGLAS G. ALTMAN Medical Statistics Laboratory. Imperial Cancer Research Fund, P.O. Box 123, Lincoln's Inn Fieldr, London WC2A 3PX. U.K.
SUMMARY Based on historical regression analyses of survival one may calculate expected survival curves for patients in a study group, based on their observed covariates. This note discusses this calculation and illustrates it on a comparison of the survival of a certain group of liver transplant patients to their prognosis under conservative treatment.
INTRODUCTION With the increased availability of published regression analyses of survival with specific diseases it has become possible to make historical comparisons between the observed survival in a (possibly fairly limited) study group and the survival that was to be expected from the published regression analyses. One area where this procedure has been used is transplantation, where randomized studies are ethically impossible, so that the best approximation to an answer to the question 'what would have happened to these patients if conservative treatment had been preferred? will be via comparison with historical controls. We discuss below some issues of the exact choice of adjustment procedure, in the context of a study on liver transplantation.'
THE DIRECT ADJUSTED SURVIVAL CURVE Neuberger et al.' reported on the survival of the first 29 patients with primary biliary cirrhosis (PBC) who had a liver transplant in the Cambridge/King's College Hospital programme. Their Kaplan-Meier survival curve is shown in Figure 1. The Kaplan-Meier curve may be regarded as an estimate of the mean survival function S,(t) = ZSi(t)/29,where Si(t)is the survival function of patient i. In the present situation an estimate was available of Si(t) under the hypothetical assumption that patient i had continued on conservative treatment rather than having been transplanted. 0277-671 5/9 1/05073346$05.00 0 1991 by John Wiley & Sons, Ltd.
Received March 1989 Revised August 1990
734
B. L. THOMSEN, N. KEIDING AND D. G. ALTMAN
i
12
21
60
months Figure 1. Observed survival’ (- . -)of 29 PBC patients who had a liver transplantation,compared with the distribution expected from an earlier clinical trial of PBC patients,’ calculatedas direct adjusted survival curve (+-) and accordingto Neuberger et al.’ (- - -)
A recent Cox regression analysis of a randomized clinical trial2 of PBC patients allowed the calculation of the estimated survival probability for time t for patient i with covariates zi as fi(t)= f(t,zi). Accordingly, the obvious estimate of the mean survival function S,(t) of the 29 transplant patients with covariate vectors zl,. . . , zZ9 is given as the mean of the individual estimated survival curves
c
q t ) = 291 291 f(t,q), -
i=
and this was termed the ‘direct adjusted survival curve’ by Makuch3 and Gail and B ~ a rMarkus ;~ et al.’ used it in the liver transplantation context. Neuberger et al. used a somewhat different adjustment procedure, averaging log(-log f(t, zf)) rather than g(t, zi). The direct adjusted survival curve and the approach by Neuberger et al. were both very recently discussed by Bonsel et ~ l . , ~ who also suggested a third method, albeit not explained in their paper, which ostensibly takes the staggered entry directly into account. Figure 1 shows the adjusted survival curves; Neuberger et al.’s suggestion represents a survival curve with less spread than the direct adjusted survival curve. In the Cox model, with standard notation, the hazard function of patient i is 4(t) = &(t)exP(Pzi)-
Thus log( - logf(t, Zi)) =
pz, + logA,(t);
where A,(t) is Breslow’s’ estimator of the cumulative hazard A,@) = j:=, I,(u)du, and bzi can be considered to be a prognostic index for a patient with covariate vector zi. The averaging suggested by Neuberger et al. produces an adjusted survival curve of S,(t) = exp(- A,(t)e@f)
which is also the predicted survival curve of a patient with prognostic index Pz equal to the average prognostic index of the patients, in particular of a patient with all covariates equal to the
A NOTE ON THE CALCULATION OF EXPECTED SURVIVAL
735
a
0.8-
r ._ D
-
,
0.6-
0 L
.
a
2"
.-
O.'-
2 , ' 7 rA
0.2-
months
b
0.0 -
4
12
21
60
months
Figure 2. Survival curves for two 'samples' consisting of the two most extreme patients (- - -) and of two 'middle' patients (-), respectively, (a) direct adjusted, and (b) adjusted according to Neuberger et al.'
average Z of the covariate values. Though they did not derive their adjustment this way, they thus arrived at the same result as if they had disregarded the sample variation. As elaborated below, this explains the above noticed difference in spread of the two adjustments. At first sight there may seem to be little difference in principle between the two averaging operations suggested above and it may require some adaptation to appreciate that s(t)accounts for the sample variability whereas SN(t)does not. To make the point clearer, we consider two 'samples', each consisting of two patients. Sample A contains the two patients with the lowest (5.12) and highest (8.36)prognostic indices, having an average of 6-74. Sample B contains the two patients with prognostic indices closest to this average, 6.66 and 6.67. Figure 2(a) shows that the direct adjusted survival curves qt)of samples A and B are very different, reflecting the large difference in sample variation between samples A and B. In contrast, Figure 2(b) shows that averaging the prognostic indices in the Neuberger et al. fashion, the adjusted survival curves &(t) are (almost) identical, depending only on the (almost) identical average prognostic indices. The difference between the survival curves estimated by the two methods depends on the variability in
736
B. L. THOMSEN, N. KEIDING AND D. G. ALTMAN
the sample of the prognostic indices. In the present example, only PBC patients with a bad prognosis were offered a liver transplant, implying a smaller variability of the prognostic indices than often encountered in Cox analyses. Recent extensive investigations in population heterogeneitysv9show that the survival distribution of the population of patients cannot have the same time-dependence of the hazard as that of the individual patients, because, briefly put, the frail die first, leading to relative reduction of population mortality over time as compared to the pattern of any individual patient. The intensity L,(t) epi of S,(t) is proportional to the underlying intensity L,(t) and thus to any individual death intensity; in contrast $(t) has intensity
where the second factor is decreasing in t, thus capturing the just mentioned aspect of the heterogeneous sample. Besides the prediction variation of an individual’s survival time around its expectation and the sample variation among the new individuals as expressed by the variation among their covariate values, a third source of variation is important: the estimation variation of the parameters of the model. It is, however, beyond the scope of this note to discuss derivation of a numerical test - or, equivalently, confidence intervals - of the difference between the observed and the direct adjusted survival curve. THE ‘EXPECTED SURVIVAL CURVE A different concept of the survival curve to be expected for a study population if individuals are assumed to be exposed to some known ‘standard’ death rates has been extensively used in connection with the grouped-time ‘actuarial’ approach to survival and was recently put into the context of continuous-time survival ana1y~is.l~ We briefly indicate the latter version in the adaptation relevant here. If the mortality of the patients followed that of the PBC patients in the clinical trial,’ the average hazard of the patients still under observation at time t would be L*(t) =
Yi(t)Lo(t)exP(B’zi)/Y(t)
where Yi(t)indicates that patient i is still at risk at time t and Y(t)= XYi(t). Defining the cumulative hazard as A*(t) = f t = o L*(u)du, the survival function S*(t) = exp{ - A*(t)]
is a continuous-time version of the ‘expected survival rate’” or ‘expected survival curve’.’’ It has a desirable property:I3 let the standard Nelson-Aalen estimator of A(t) be defined by
where T , < T2 < , . . are the times of (observed) deaths. Then A*(t) - A(t) has expectation zero, under the null hypothesis that patient i has hazard L,(t). Therefore S*(t) represents, under the null hypothesis, an expected survival curve. Outside of the null hypothesis it is not so easy to interpret S*(t), containing as it does information on study group mortality, through Yi(t),as well as
737
A NOTE ON THE CALCULATION OF EXPECTED SURVIVAL
0
4
12
28
60
months Figure 3. Expected survival according to Neuberger et al.' (- - -), as direct adjusted survival curve (-) survival curve'" (- . . -)
and as 'expected
standard mortality, through l , ( t ) .S*(t) therefore cannot be recommended as an expected survival function, despite its wide use for this purpose. A similar difficulty outside of the null hypothesis was recently discussed for the well-known person-years method. l4 S*(t) is, however, useful for inference on excess mortality, as follows. In the particular case where the study group has an additive excess mortality a(t) over the standard, that is, patient i has hazard a(t) + ,Ii([), one may obtain an estimate of the corresponding 'survival' function
as the so-called relative survival function s(t)/S*(t).Andersen and Vaeth' showed that results on unbiasedness, consistency and asymptotic normality are available. In the present example S*(t) is even steeper than S,(t) (see Figure 3). Our explantion is that mortality after transplantation depends on different risk factors from mortality under medical treatment, so that the order in which patients actually die is different from what it 'should be' according to l i ( t ) .This also implies that the simple additive model for the hazard is inadequate in the present example. Whatever the explanation, the case documents the drawback of using the 'expected survival rate' S * ( t ) rather than the immediately interpretable 'direct adjusted survival curve' Q t ) . Our conclusion is that for the specific purpose of calculating 'what the survival would have been if some given individual-specific mortality rates were operating', the direct adjusted survival curve should be preferred. ACKNOWLEDGEMENTS
This study was supported by the Danish Natural Science Council grant no. 11-7575. REFERENCES
1. Neuberger, J., Altman, D. G., Christensen, E., Tygstrup, N. and Williams, R. 'Use of a prognostic index in evaluation of liver transplantation for primary biliary cirrhosis', Transplantation, 4, 71 3-716 (1986).
738
B. L. THOMSEN, N. KEIDING AND D. G . ALTMAN
2. Christensen, E., Neuberger, J., Crowe, J., Altman, D. G., Popper, H., Portmann, B., Doniach, D..
Ranek, L., Tygstrup, N. and Williams, R. ‘Beneficial effect of azathioprine and prediction of prognosis in primary biliary cirrhosis. Final results of an international trial’, Gastroenterology, 89, 1084-1091 (1985). 3. Makuch, R. W. ‘Adjusted survival curve estimation using covariates’, Journal of Chronic Diseases, 3, 437-443 (1982). 4. Gail, M. H. and Byar, D. B. ‘Variance calculations for direct adjusted survival curves with applications to testing for no treatment effect’, Biometrical Journal, 28, 587-599 (1986). 5. Markus, B. H., Dickson, E. R., Grambsch, P. M., Fleming, T. R., Mazzaferro, V., Klintmalm, G. B. G.,
Wiesner, R. H., van Thiel, D. H. and Starzl, T. E. ‘Efficacy of liver transplantation in patients with primary biliary cirrhosis’, The New England Journal of Medicine, 320, 1709-1713 (1989). 6. Bonsel, G. J., Klompmaker, I. J., van’t Veer, F., Habbema, J. D. F. and Slooff, M. J. H. ‘Use of prognostic models for assessment of value of liver transplantation in primary biliary cirrhosis’, “he Lancet, 335, 493-497 (1990). 7. Breslow, N. Contribution to the discussion of D. R. Cox ‘Regression models and life tables’, Journal of the Royal Statistical Society, Series B, 34, 216217 (1972). 8. Hougaard, P. ‘Life table methods for heterogeneous populations: Distributions describing the heterogeneity’, Biometrika, 71, 75-83 (1984). 9. Aalen, 0.0. ‘Heterogeneity in survival analysis’, Statistics in Medicine, 7 , 1121-1137 (1988). 10. Ederer, F., Axtell, L. M. and Cutler, S. J. ‘The relative survival rate: A statistical methodology’, National Cancer Institute Monographs, 6, 101-121 (1961). 11. Lee, E. T. Statistical Models for Survival Data Analysis, Lifetime Learning Publications, Belmont, California, 1980, Chapter 4. 12. Hill, C., Laplanche, A. and Rezvani, A. ‘Comparison of the mortality of a cohort with the mortality of a reference population in a prognostic study’, Statistics in Medicine, 4, 295-302 (1985). 13. Andersen, P. K. and Vaeth, M. ‘Simple parametric and nonparametric models for excess and relative mortality’, Biometrics, 45, 523-535 (1989). 14. Keiding, N. and Vaeth, M. ‘Calculating expected mortality’, Statistics in Medicine, 5, 327-334 (1986).
