Comparison of 2 Methods for Calculating Adjusted Survival Curves From Proportional Hazards Models

Adjusted survival curves are often included in published research articles to present the anticipated survival of 2 or more comparison groups in the circumstance of balanced covariate risk factors between groups. The most widely used method for generating adjusted survival curves from Cox proportional hazards models is referred to as the mean of covariates method.¹ This method is simple to implement but has been considered by some^{1 - 2} to be both mathematically and conceptually problematic. Critics point to the availability of a better method for calculating such curves called the corrected group prognosis method.^{1 - 4}

A search of the literature reveals several examples of studies that present adjusted survival curves^{5 - 9} or estimates of risk-adjusted survival.^{10 - 11} Some of these studies^{5 - 7} explicitly mention that the suboptimal mean of covariates method was used, while others^{8 - 9} fail to describe the method used to derive adjusted survival curves.

In this statistical report, we provide brief descriptions of the mean of covariates and corrected group prognosis methods and then apply the 2 methods to calculate risk-adjusted survival curves for patients with and without diabetes undergoing cardiac catheterization.¹² We discuss the discrepant results from these analyses and provide access to programs for the implementation of the corrected group prognosis method in 3 statistical packages.

In the widely used mean of covariates method, mean values for covariates are inserted into the survival function of the proportional hazards model. For continuous variables like age, the mean value among study patients is used. For dichotomous covariates, a value between 0 and 1 is used to reflect the proportion of patients in the database with the condition (eg, 0.17, if 17% have the condition).

In the corrected group prognosis method, survival curves are first calculated for each of the unique combinations of covariates in a database, based on the coefficients from a single proportional hazards model developed on the entire database. A weighted average of these individual curves is then calculated, with weights proportional to the number of individuals at each level of covariates. For models that contain a large number of covariates, there may be many covariate combinations for which individual curves must be calculated. For the statistical case report that we present 2419 survival curves had to be calculated and averaged to yield the final adjusted curves that are presented in the "Results" section. We have developed programs that readers can download from the Internet¹³ for implementing the corrected group prognosis method in SAS, STATA, and S-Plus.

We recently published an article that presents adjusted survival curves after cardiac catheterization for individuals with and without diabetes undergoing cardiac catheterization between January 1, 1995, and December 31, 1996.¹² The purpose of this analysis was to provide information on the adjusted survival experience of individuals with diabetes, while controlling for age, sex, and comorbidities. Details of that study's methods are provided elsewhere.^{12 ,14}

In our analysis, we used a proportional hazards analysis to compare survival up to 3 years after cardiac catheterization by diabetes status while controlling for the other clinical variables that were significant predictors of survival (P≤.05). The proportional hazards assumption was evaluated graphically and found to be appropriate. Risk-adjusted survival curves were then plotted from the proportional hazards model using first the mean of covariates method then the corrected group prognosis method.

Table 1 presents the clinical characteristics of the 11 468 patients studied along with hazard ratios for all variables that were statistically significant in the proportional hazards model. Patients with diabetes generally had a higher prevalence of various clinical risk variables.

The unadjusted survival estimates to 1044 days of follow-up (the last day on which a death occurred) were 86.8% and 93.2% for patients with and without diabetes, respectively (P<.001). The unadjusted hazard ratio for diabetes was 2.0 (95% confidence interval [CI],1.7-2.3), and with adjustment, the hazard ratio decreased to 1.2 (95% CI, 1.0-1.4).

To graphically represent this decrease, we first applied the mean of covariates method and found that the adjusted curves fell above the unadjusted curves for patients with and without diabetes (Figure 1). Adjusted survival to 1044 days was 94.1% and 94.9%, respectively, according to the mean of covariates method.

Figure 2 presents the results of the corrected group prognosis analysis. In this instance, the adjusted survival curves are more appropriately positioned between the unadjusted curves, with survival at 1044 days of 91.3% and 92.4%, respectively.

The marked discrepancy between the results of these 2 statistical analyses undoubtedly relates to the previously described limitations of the mean of covariates method.^{1 - 2} These include the assignment of mean covariate values between 0 and 1 for dichotomous variables (eg, 0.07) that are meaningless at the individual level and the recognition that the method calculates the hazard for a hypothetical average individual rather than a population-averaged value.

A simple explanation for the discrepant results is that the averaging of covariate values occurs at different parts of the proportional hazards survival function. In the mean of covariates method, the averaging occurs within the function's exponent, whereas it is actual survival curves that are averaged in the corrected group prognosis method. We present an analogous example for clarification: the mean of the numbers 1, 1, and 4 is equal to 2. If we then put this mean value in the exponent of the base number 10 (ie, 10²), the result is 100. In contrast, if we calculate the mean of 10¹, 10¹, and 10⁴, the result is 3340. The first calculation represents averaging within an exponent—analogous to the mean of covariates method. The second example corresponds to the type of averaging that occurs in the corrected group prognosis method (ie, averaging of expressions that include exponents).

To further explore the circumstances under which distortion of curves is likely to be greatest, we conducted sensitivity analyses for which we controlled for only 1 covariate at a time. This allowed us to determine that the distortion of the mean of covariates method is greatest when the covariate(s) being controlled for is prevalent in the database and when the hazard ratio(s) associated with the covariates is large.

This statistical report provides an important message to researchers, readers, and journal editors. Researchers should consider abandoning the mean of covariates method despite its relative simplicity and the availability of macros for its implementation in some statistical packages. The limitations of the method are described,^{1 - 2} and this case report demonstrates how misleading curves can arise.

Readers of the medical literature need to recognize that caution is necessary when interpreting adjusted survival curves that are presented without corresponding unadjusted curves. Caution is also required when articles fail to mention the method used to calculate adjusted curves or when authors explicitly state that the mean of covariates method was used. Journal editors and reviewers should also take note of this issue when reviewing papers that present adjusted survival curves.

Although the corrected group prognosis method is more complex, we provide an Internet site¹³ from which programs can be downloaded for use in SAS, STATA, or S-Plus. The programs are easy to use and can be applied either to a sample database provided on the Web site or to other databases. We hope that these programs and our report will heighten awareness regarding this important statistical issue.