Review and evaluation of performance measures for survival prediction models in external validation settings

Background When developing a prediction model for survival data it is essential to validate its performance in external validation settings using appropriate performance measures. Although a number of such measures have been proposed, there is only limited guidance regarding their use in the context of model validation. This paper reviewed and evaluated a wide range of performance measures to provide some guidelines for their use in practice. Methods An extensive simulation study based on two clinical datasets was conducted to investigate the performance of the measures in external validation settings. Measures were selected from categories that assess the overall performance, discrimination and calibration of a survival prediction model. Some of these have been modified to allow their use with validation data, and a case study is provided to describe how these measures can be estimated in practice. The measures were evaluated with respect to their robustness to censoring and ease of interpretation. All measures are implemented, or are straightforward to implement, in statistical software. Results Most of the performance measures were reasonably robust to moderate levels of censoring. One exception was Harrell’s concordance measure which tended to increase as censoring increased. Conclusions We recommend that Uno’s concordance measure is used to quantify concordance when there are moderate levels of censoring. Alternatively, Gönen and Heller’s measure could be considered, especially if censoring is very high, but we suggest that the prediction model is re-calibrated first. We also recommend that Royston’s D is routinely reported to assess discrimination since it has an appealing interpretation. The calibration slope is useful for both internal and external validation settings and recommended to report routinely. Our recommendation would be to use any of the predictive accuracy measures and provide the corresponding predictive accuracy curves. In addition, we recommend to investigate the characteristics of the validation data such as the level of censoring and the distribution of the prognostic index derived in the validation setting before choosing the performance measures.


Background
Prediction models are often used in the field of healthcare to estimate the risk of developing a particular health outcome. These prediction models have an important role in guiding the clinical management of patients and monitoring the performance of health institutions [1,2]. For example, models have been developed to predict the risk of in-hospital mortality following heart valve surgery and to predict the risk of developing cardiovascular disease within the next 10 years [3,4]. Given their important role in health research, it is essential that the performance of a prediction model is evaluated in data not used for model development, using appropriate statistical methods [5,6]. This model evaluation process is generally termed 'model validation' [7,8]. The general idea of validating a prediction model is to establish that it performs well for new patients. Different types of validation process have been discussed in the literature [5][6][7][8].
The most commonly used processes include (i) splitting a single dataset (randomly or based on time) into two parts, one of which is used to develop the model and the other used for validation, (internal or temporal validation) and (ii) validating the model using a new dataset collected from a relevant patient population in different centres (external validation). Of the two approaches, external validation investigates whether a prediction model is transportable (or generalisable) to new patients.
When validating a prediction model, the predictive performance of the model is commonly addressed by quantifying: (i) the 'distance' between the observed and predicted outcomes (overall performance); (ii) the ability of the model to distinguish between low and high risk patients (discrimination); (iii) the agreement between the observed and predicted outcomes (calibration) [8]. Performance measures based on these concepts are well established for risk models for binary outcomes [1,9,10], but that is not the case for risk models for survival outcomes (survival prediction models) where censoring complicates the validation process [6].
Several performance measures have been suggested for use with survival prediction models. However, a few of these are not appropriate for use with validation data without modification. Also, some require specification of a clinically appropriate time-point or region to match the aims of the validation study. Some of these performance measures have been reviewed previously [11][12][13][14][15][16][17], although only two of the reviews were in the context of model validation. These were Hielscher et al. who reviewed 'overall performance' measures, and Schmid and Potapov who reviewed discrimination measures [12,16]. Consequently, it is still unclear which performance measures should be routinely used in practice when validating survival prediction models using external data.
A good performance measure should be unbiased in the presence of censoring in the validation data. If this were not the case, the level of censoring would affect the evaluation of model performance and a high level of censoring might lead to an over-optimistic verdict regarding the performance of the prediction model. In addition, a good measure should be straightforward to interpret and, ideally, should be easy to implement or available in widely used software.
The aim of this paper is to review all types of performance measures (overall performance, discrimination and calibration) in the context of model validation and to evaluate their performance in simulation datasets with different levels of censoring and case-mix. Where necessary, measures have been modified to allow their use with validation data and a case study is provided to describe how these measures can be estimated in practice. Recommendations are then made regarding the use of these measures in model validation.

Data
Two datasets, which have previously been used to develop clinical prediction models, were used as the basis of the simulation study. They differ with respect to event rates, level of censoring, types of predictors and amount of prognostic information.

Breast cancer data
This dataset contains information on 686 patients diagnosed with primary node positive breast cancer from the German Breast Cancer Study [18]. The outcome of interest is recurrence-free survival time and the dataset contains 299 (44%) events. The median follow-up time is 3 years. The predictors are age, number of positive lymph nodes, progesterone receptor concentration, tumour grade (1-3), and hormone therapy (yes/no). These data have been analysed previously by Sauerbrei and Royston and their Model III was used as the basis for simulation [19]. That is, the continuous predictors were all transformed using fractional polynomials (FPs) and tumour grade was dichotomised (1/2-3). Number of positive lymph nodes and progesterone receptor concentration were each modelled using one FP term whereas age was modelled using two FP terms.

Hypertrophic cardiomyopathy data
This dataset contains information on a retrospective cohort of 1504 patients with hypertrophic cardiomyopathy (HCM) from a single UK cardiac centre [20]. The outcome of interest is sudden cardiac death or an equivalent event, i.e., a composite outcome and the dataset contains just 84 (6%) events. The median follow-up time is over 6 years. The predictors of interest are age, maximal wall thickness, left atrial diameter, left ventricular outflow gradient, family history of sudden cardiac death (yes/no), non-sustained ventricular tachycardia and unexplained syncope (yes/no). The prediction model produced by O'Mahony et al. was used as the basis for simulation [20]. In particular, maximal wall thickness was modelled using linear and quadratic terms.

Prediction models for survival outcomes
Prediction models for survival outcomes are commonly developed using the Cox proportional hazards model and hence this model was used in the simulations [5,21]. The Cox model models the hazard h(t|x) at time t as a product of a nonparametric baseline hazard h 0 (t) and the exponential of the prognostic index η = β 1 x 1 + … + β p x p = β T x. The latter is a linear combination of p predictor values with regression coefficients β 1 , …, β p providing weights. The predictive form of this model can be written in terms of the survival function as where S(t|x) is the probability of surviving beyond time t given predictors x, and S 0 (t) is the baseline survivall function at time t, where S 0 (t) = exp[−∫ 0 t h 0 (u)du]. To make predictions at a specific time-point τ, one requires estimatesβ andŜ 0 τ ð Þ.

Performance measures for survival prediction models
Measures were selected for investigation on the basis of their performance in previous reviews [11][12][13][14][15][16], their ease of interpretation, and their availability, or ease of implementation, in statistical software. The selected performance measures are now described in the context of model validation. All measures were implemented in Stata using either built-in or user-written routines [22].

Measures of overall performance
Six measures of 'overall performance' were selected, of which four are based on predictive accuracy and two on explained variation [14]. These 'R 2 -type' measures typically take values between 0 and 1, though negative values may be possible in validation data if the prediction model is out-performed by the null model that has no predictors. This issue is discussed later. The measures based on predictive accuracy were Graf et al's R 2 measure and its integrated counterpart [23], Schemper and Henderson's R 2 measure [24], and a modified version of the latter based on Schmid et al. [25]. The measures based on explained variation were Kent and O'Quigley's R 2 PM [26], and Royston and Sauerbrei's R 2 version of their separation statistic D [27]. Nagelkerke's R 2 measure was not considered due to its known poor performance in the presence of censoring [26,28].

Graf et al's R 2
BS and R 2

IBS
The R 2 BS measure proposed by Graf et al. is based on quantifying prediction error at a time-point τ using a quadratic loss function [21]. Specifically, R BS is the prediction error at time τ for the prediction model and I(T > τ) is the individual survival status at this timepoint. Similarly, BS(τ) is the prediction error for the null model at the same time-point, and is based on the survival function Ŝτ ð Þ from the null model. The integrated version, R IBS 2 (τ), is defined in a similar way to R BS 2 (τ) but involves integrating both BS(t|x) and BS(t) over the range [0, τ].
The calculation of prediction errors for both of these models in validation data requires estimates of the corresponding baseline survival functions. This, however, is rarely provided by model developers [6]. One solution might be to estimate these survival functions by refitting the Cox model with the PI as the sole predictor in the validation data. This is the approach that we took when calculating R 2 BS and R 2 IBS , and the R 2 SH and R 2 S measures described below.
Schemper and Henderson's R 2 SH and Schmid et al's R 2

S
The R 2 measure proposed by Schemper and Henderson (denoted here by R 2 SH ) is similar to Graf et al's R 2 IBS but is based on an absolute loss function [24]. This loss function was chosen to reduce the impact of poorly predicted survival probabilities, which are likely to occur in the right tail of the survival distribution. Specifically, is the prediction error at time τ for the prediction model and W(τ) = 1/∫ 0 τ f(t)dt is a weight function to compensate for the measure being defined only on (0, τ). Similarly, D(τ) is the prediction error for the null model. Schmid et al. prove that Schemper and Henderson's estimator of D(τ|x) and D(τ) is not robust to model misspecification and suggest an improved estimator [25]. We estimated a summary measure, denoted by R 2 S , based on this estimator. Kent [26]. That is, seeks to quantify the proportion of variation in the outcome explained by the predictors in the prediction model, where σ ϵ 2 ≅ π 2 /6 is the variance of the error term in an equivalent Weibuill model [13].
This measure does not use the observed survival times directly in its calculation and instead relies on the prediction model being correctly specified. As a result, R 2 PM could be misleading if an apparent strong relationship between the outcomes and predictors in development data is not reproduced in validation data. To overcome this deficiency, we suggest re-calibrating the prediction model to the validation dataset before calculation of R 2 PM . This procedure will tend to reduce the value of R 2 PM and is described later.

Royston and Sauerbrei's R 2 D
Royston and Sauerbrei's R 2 D is similar to R 2 PM but is based on the authors' own D statistic, a measure of prognostic separation described later. That is, [27]. The ratio D 2 /κ 2 is an estimator of Var(η), provided that η is Normally distributed.

Measures of discrimination
Four measures of discrimination were selected, of which three are based on concordance and one on prognostic separation. Discrimination measures assess how well a model can distinguish between low-and high-risk patients, and concordance measures in particular quantify the rank correlation between the predicted risk and the observed survival times. Concordance measures usually take values between 0.5 and 1, where a value of 0.5 indicates no discrimination and a value of 1 indicates perfect discrimination. The selected concordance measures were those of Harrell [29], Uno et al. [30], and Gönen and Heller [31], and the selected prognostic separation measure was Royston and Sauerbrei's D statistic [27].

Harrell's C H
The concordance probability is the probability that of a randomly selected pair of patients (i,j), the patient with the shorter survival time has the higher predicted risk. Formally, where η i and η j are the prognostic indices for patients i and j, and T i and T j are the corresponding survival times. Harrell's estimator C H considers all usable pairs of patients for which shorter time corresponds to an event and estimates C H as the proportion of these pairs for which the patient with the shorter survival time has the higher predicted risk [31]. A modified version of this estimator, C H (τ), restricts the calculation to include just those patient pairs where T i < τ and may provide more stable estimates [29,30]. This truncated version may also be preferred if one were primarily interested in the discrimination of a prediction model over a specified period, for example within 5 years [20].

Uno et al's C U
In the presence of censoring C H and C H (τ) are biased, even under independent censoring, as they ignore patient pairs where the shorter observed time is censored [15,32]. Due to this deficiency, Uno et al. [30] proposed a modified estimator C U (τ) that uses weightings based on the probability of being censored. Furthermore, like C H (τ), the calculation may also be restricted to include just those patient pairs where T i < τ. Uno et al. found that their estimator was reasonably robust to the choice of τ, but noted that the standard error of the estimate could be quite large if τ were chosen such that there was little follow-up or few events beyond this time point [29].

Gonen and Heller's C GH
Gönen and Heller proposed an alternative estimator C GH based on a reversed definition of concordance [30], which is the probability that of a randomly selected pair of patients (i, j), the patient with the higher predicted risk has the shorter survival time. To avoid bias caused by censoring, their estimator is a function of the model parameters and the predictor distribution and assumes that the proportional hazards assumption holds. As with R 2 PM , C GH does not use the observed event and censoring times in its calculation and relies on the prediction model being correctly specified [15]. Therefore, we suggest re-calibrating the prediction model to the validation dataset before calculating C GH .

Royston and Sauerbrei's D
The D statistic is a discrimination measure that quantifies the observed separation between subjects with low and high predicted risk [27]. Specifically, D estimates κσ, where σ is the standard deviation of the prognostic index and κ ¼ ffiffiffiffiffiffiffiffi 8=π p . The scale factor κ enables D to be interpreted as the log hazard ratio that compares two equal-sized risk groups defined by dichotomising the distribution of the patient prognostic indices at the median value.

Measures of calibration
One calibration measure was selected, the commonly used calibration slope proposed by van Houwelingen [33], which is based on the analogous measure for binary outcomes [34,35].

Calibration slope
The calibration slope is simply the slope of the regression of the observed survival outcomes on the predicted prognostic index [33]. It is estimated by fitting a Cox model to new survival outcomes with the predicted prognostic index, η, as the sole predictor in the model Values of α̂1 close to 1 suggest that the prediction model is well calibrated. Moderate departures from 1 indicate that some form of model re-calibration may be necessary. In particular, α̂1≪1 suggests over-fitting in the original data with predictions that may be too low for low risk patients or too high for high risk patients.
A brief summary of the performance measures is given in Table 1.

Case study to illustrate the performance measures
A case study is now presented using the breast cancer data in order to describe how the performance measures may be evaluated in a validation setting. The dataset was split randomly into two parts with two thirds of the data used for model development and one third used for model validation. A Cox model was fitted to the development data using the same predictors as in Sauerbrei and Royston's Model III [19] and the predicted prognostic index was calculated for all patients in the validation data using the estimated regression coefficients β̂. The values of all performance measures are shown in Table 2 with 95% confidence intervals estimated using the bootstrap techniques based on 200 bootstrap samples. For those measures that require specification of a time-point τ, 3 years was deemed to be clinically appropriate. This was also the median follow-up time.
The estimated prediction errors used to estimate R 2 BS and R 2 IBS are shown in Fig. 1a. The errors for the prediction and null models are similar for the first 12 months after which the superiority of the prediction model is evident. The corresponding prediction errors used to estimate R 2 SH appear similar in shape although the magnitude of the errors are larger due to the use of an absolute loss function (Fig. 1b). The prediction errors used to estimate R 2 S are almost indistinguishable from those used to estimate R 2 SH (results not shown). R 2 IBS , R 2 SH and R 2 S were estimated after averaging the prediction errors over the first 3 years. As expected R 2 SH and R 2 S are very similar, and both are slightly larger than R 2 IBS . R 2 BS was estimated using just the prediction errors at 3 years in Fig. 1a. Its value is larger than that of R 2 IBS as the separation between the prediction errors is close to its maximum at this time-point.
To estimate R 2 PM in the valids first re-calibrated for reasons explained earlier. That is, the Cox model h tj ηð Þ ¼ h 0 t ð Þ exp α ηð Þ was fitted to the validation data, where ηî s the predicted prognostic index calculated using the regression coefficients estimated in the development data. R 2 PM was then estimated using α̂2 Var ηð Þ rather than Var ηð Þ . No such re-calibration is required to estimate R 2 D since, unlike R 2 PM , it uses the observed survival outcomes in the validation data. The values of R 2 PM and R 2 D are very similar and noticeably larger than the measures based on predictive accuracy ( Table 2). We note that a naïve calculation (without re-calibration) of R 2 PM would have produced a much larger value of 0.292, which would have provided an over-optimistic quantification of the model's predictive performance.
The concordance measures C H and C U were estimated using all usable pairs in which shorter time corresponds to an event, and C GH was estimated after re-calibrating the prediction model to the validation data as described above. The values of these 3 measures are all reasonably similar and would lead to similar conclusions in practice ( Table 2). A naïve estimation of C GH (without recalibration) would have produced a much larger value of 0.696. Restricting the estimation of C H and C U by censoring survival times in the validation data at 3 years produces slightly higher values for both measures, suggesting that the risk model has slightly better discrimination when considering survival over just the first 3 years. The D statistic suggests that the prediction model provides a reasonably high amount of prognostic separation (Table 2). Specifically, if one were to form two risk groups of equal size in the validation data, then the corresponding hazard ratio would be exp(0.998) = 2.71. The calibration slope estimate of 0.76 (equal to the α̂estimated during the re-calibration process above) suggests that the prediction model has been slightly over-fitted. We note that, in practice, one can detect and adjust for model over-fitting during model development.
The selected measures all provide useful information. R 2 IBS, R 2 SH, and R 2 S, provide a summary measure quantifying the improvement in predictive accuracy offered by the prediction model over the null model. The R 2 BS measure is more appropriate if one is interested in predictive accuracy at a specific time-point, which is sometimes the case in practice. The prediction error curves provide additional insight into the performance of the prediction model at different time-points. R 2 PM and R 2 D , which both quantify explained variation, produced very similar values though calculation of R 2 PM required a re-calibration of the prediction model. The concordance measures C H , C U and C GH produced similar estimates, though calculation of C GH required the prediction model to be re-calibrated. Additionally, if required, the calculation of C H and C U can be restricted which may be appropriate if one wishes to quantify the discrimination of a prediction model before a specified time-point. Finally, the D statistic produces an intuitive quantification of prognostic separation and the calibration slope provides a succinct indication of whether the prediction model is over-fitted or not.

Evaluation of the performance measures using simulation
Following the case study, a simulation study was performed to investigate how robust the measures are with respect to both censoring and the characteristics of the validation data. The simulation design is now described.

Simulation scenarios
The simulation study is based on the breast cancer and HCM datasets described earlier. For both datasets, development and validation datasets were generated by simulating new outcomes based on a true model and combining these with the original predictor values. Models were fitted in the development data and the performance measures estimated in the validation data. Measures which require a choice of time-point τ (including C H and C U ) used 3 years for the breast cancer data and 5 years for the HCM data. These values were chosen as they are close to the respective median follow-up times in the original datasets and are conventional choices for survival data. In practice, the choice of time-point would be clinically motivated and based on the underlying research question.
The performance measures were investigated over a range of scenarios to mimic real situations. For all simulations, validation data were constructed to have one of three different risk profiles (denoted low, medium, and Range: 0 to ∞ Interpretation: log hazard ratio between two equal sized prognostic groups fromed by dichotomising the PI at its median..

Available in Stata and easy to implement in other software
Calibration Cal Slope Regression slope of the PI and assesses the agreement between the observed and predicted survival..
Range: −∞ to ∞ Interpretation: a value of 1 suggests perfect calibration and a value much lower than 1 suggest overfitting.
Easy to implement in any software high). The use of different risk profiles reflect the fact that, in practice, the characteristics of the patients in the development and validation data may differ [36]. In particular, the event rate for patients in the validation data may be higher or lower than that for patients in the development data due to differences in case-mix. Four levels of random censoring were considered for the validation datasets (0, 20, 50, and 80%) which combined with the risk profiles, results in a total of 12 validation scenarios for each clinical dataset. The development datasets were generated with no censoring. 5,000 pairs of development and validation datasets were generated for each scenario.

Generating new survival and censoring times
Survival times were generated using the Weibull distribution as below where η and γ are the observed regression prognostic indices and shape parameter respectively (both used here as the proxy of the true values) and u is a uniformly distributed random variable on (0, 1). For the breast cancer data, the prognostic indices and shape parameter were obtained by fitting a Weibull proportional hazards model using the same predictors as in Sauerbrei and Royston's Model III [19]. For the HCM data, the prognostic indices were based on the regression coefficients estimated by O'Mahony et al. [20] and just the shape parameter was estimated using a Weibull model with the prognostic index specified as an offset.
To introduce random censoring, additional Weibull distributed censoring times were simulated using t c = (−log(u)/ λ) 1/γ where different choices of the scalar λ were used to give different proportions of censoring. A subject was considered to be censored if their censoring time was shorter than their survival time.

Generating validation data with different risk profiles
The three different risk profiles were created in the validation data, by first splitting the patients into tertile risk groups based on their true prognostic index η. It is assumed that the first tertile group consists of low-risk patients, the second medium risk, and the third high- a) low risk profile: 80% of the patients were sampled (without replacement) from the low-risk group, 50% from the medium-risk group, and 20% from the high-risk group; b) medium risk profile: 50% of the patients were sampled from the low-risk group, 50% from the medium-risk group, and 50% from the high-risk group; c) high risk profile: 20% of the patients were sampled from the low-risk group, 50% from the medium-risk group, and 80% from the high-risk group.
This sampling procedure was performed before generating each validation dataset and resulted in validation datasets that were half the size of the original datasets. In contrast, no sampling of patients was performed when generating the development datasets; all patients were used.
The prognostic indices were approximately normally distributed for all risk profiles and for both datasets. There was slight skewness, particularly in the low and medium risk profile datasets. For example, the (average) skewness in the HCM datasets was 0.8 (low risk profile), 0.4 (medium) and 0.1 (high). There was a similar trend in the breast cancer datasets, although the values were lower. The variance was largest for the medium profile datasets which was to be expected considering the sampling scheme. This suggests that the medium profile datasets contained more prognostic information than the low and high profile datasets. Finally, there was more prognostic information in the breast cancer data, as evidenced by the wider range of the corresponding prognostic indices. Table 3 shows the mean values of the overall performance measures over 5000 simulations for the breast cancer data, for the four levels of censoring and three risk profiles. The three summary measures based on predictive accuracy (R 2 IBS , R 2 SH and R 2 S ) produced very similar values and were all unaffected by censoring. The values of these measures were highest for the medium risk profile simulations, where the patient characteristics were essentially the same in the development and validation samples, and lowest for the low risk profile simulations. Variability increased with increasing censoring, as expected, and was highest for R 2 IBS . This can clearly be seen in Fig. 2 which shows the distribution of the values of the performance measures over the 5000 simulations (a few negative values were deleted to aid clarity). R 2 BS , evaluated at 3 years, was also unaffected by censoring and achieved higher values in the medium risk profile simulations. R 2 BS also produced some negative values (4%) when censoring was 80%. The two measures based on explained variation (R 2 PM and R 2 D ) produced similar values that were twice as large as the values obtained for R 2 IBS , R 2 SH and R 2 S . R 2 PM was unaffected by censoring but R 2 D increased slightly as censoring increased. The relationships between the various overall performance measures are shown in Fig. 3 for the medium risk profile scenario. In particular, there was excellent agreement between the R 2 SH and R 2 S measures which weakened as censoring increased (ρ = 0.54 for 80% censoring). Also, there was good agreement between R 2 PM and R 2 D which seemed little affected by censoring (ρ = 0.95). Very similar relationships were seen for the low and high risk scenarios (results not shown). Table 4 shows the mean values of the discrimination and calibration measures for the breast cancer data. The Harrell and Uno c-indices were estimated twice, first using all usable patient pairs (C H and C U ) and second by restricting the calculations by censoring times greater than 3 years (C H (3) and C U (3)). For 0% censoring, the C H and C GH mean values were very similar, and the C H and C U estimates were identical by definition. C H tended to increase as censoring increased, whereas C U and particularly C GH were little affected. C GH was the least variable of these three estimates. The variability in C H and C U was similar except for 80% censoring where the variability in C U was far larger (Fig. 4). The increased variability was probably due to large values of the weights caused by the high degree of censoring. The mean value of C H (3) (and C U (3)) was slightly larger than that for C H which suggests that the models were better able to discriminate within the first 3 years compared to across the whole follow-up period. Both C H (3) and C U (3) were relatively stable with respect to censoring, and the variability of both measures was similar. The calibration slope and particular the D statistic showed a slight tendency to increase with censoring. The relationships between the discrimination measures are shown in Fig. 5 for the medium risk profile scenario. In particularly, there was reasonable agreement between the concordance measures and D. The strong relationship between C H and C H (3) for 80% censoring is explained by the fact that there were few observed failure times above 3 years with this level of censoring. Very similar relationships were seen for the low and high risk breast cancer scenarios (results not shown).

Simulation results
The results for the overall performance measures for the HCM data can be seen in Table 5 and Fig. 6. The mean values were all lower than the corresponding values for the breast cancer data. In particular, the predictive accuracy values were considerably lower due to the relatively low number of events (5%) that occurred before 5 years. In addition there were many negative R 2 IBS and R 2 BS values (11 and 9% respectively). R 2 D was affected by censoring in the low and medium risk profile simulations which may be explained by skewness in the prognostic index [6]. For example, the prognostic index was most skewed in the low risk profile HCM simulations, which is where greatest effect of censoring was observed. The relationships between the overall performance measures were similar, and often slightly stronger, than those seen in the breast cancer simulations (results not shown).
The results for the discrimination and calibration measures for the HCM data can be seen in Table 6 and Fig. 7. As with the overall performance measures, the discrimination values were all lower than the corresponding values for the breast cancer data. C H was again badly affected by censoring. In addition, D, like R 2 D , was also affected by censoring in the low and medium risk profile simulations. Also notable is the increased variability of the C H (5) and C U (5) measures compared to their unrestricted counterparts. This again is due to the relatively  Fig. 2 Box plots showing the distribution of the overall performance measures for 3 risk profiles (low, medium and high) and 4 levels of censoring (0, 20, 50 and 80%) for the breast cancer data over 5000 simulations low number of events that occurred before 5 years and the consequent low number of patient pairs used to estimate both measures. Again, the relationships between the discrimination measures were similar to those seen in the breast cancer simulations (results not shown). In particular, there was excellent agreement between C GH and D (ρ = 0.99).

Discussion
The aim of this research was to review some of the promising performance measures for evaluating prediction models for survival outcomes, modify them if necessary for use with external validation data, and perform a simulation study based on two clinical datasets in order to make practical recommendations.
Measures based on predictive accuracy quantify the predictive ability of the prediction model, relative to a null model with no predictors, on a percentage scale and can be readily communicated to health researchers. The measures investigated in this study (R 2 IBS , R 2 BS , R 2 SH and R 2 S ) may be estimated for any survival prediction model provided that both the prognostic index and baseline survival function are available, although R 2 SH also implicitly assumes that the model is correctly specified [12]. If the baseline survival function is not available, which is usually the case in practice [6], then one approach might be to   estimate it using the validation data. This is a pragmatic choice as the baseline survival function is rarely presented in practice by model developers. An alternative, arguably better, approach would be to estimate the baseline survival function for the prediction model (with covariates), but not the null model, using the development data. This alternative approach was investigated in the case study and produced very similar results (not shown). A negligible difference in the baseline survival function is also reported in [6] when it was estimated using development and validation data separately. Similarly, for predictions from a null model, which are required for these measures, we suggest using the Kaplan-Meier estimate from the validation data. Again for these measures, a choice of time-point is also required since the summary measures (R 2 IBS , R 2 SH and R 2 S ) are estimated over a specified range and R 2 BS is estimated at a specified time-point. In practice, the choice of time-point will be guided by the clinical research question and the length of follow-up. For  [12].
The measures based on explained variation (R 2 PM and R 2 D ) may be estimated for any proportional hazards model provided that the prognostic index is available, although we suggest that the prediction model is re-calibrated to the validation data before calculation of R 2 PM to ensure that the survival times in the validation data are used in its calculation. Both measures provided very similar values in our simulations. R 2 PM was robust to censoring, but R 2 D tended to increase with censoring if the prognostic index was skewed.  Concordance measures are routinely used in practice since the concept of correctly ranking patient pairs can be readily communicated to health researchers [5]. C H and C U can be estimated for any survival prediction model that is able to rank patients. In addition, the calculation of C U may also be restricted to a specified range of time, which may be useful to match with clinical aims or to compare concordance across different datasets. The calculation of C H may also be restricted though it is not clear how often this is done in practice [37]. C GH has a similar interpretation to the other concordance measures but requires that the model is correctly specified. As with R 2 PM , we suggest that the prediction model is re-calibrated to the validation data before calculation of C GH to ensure that the survival times in the validation data are used. Harrell's C H , in its unrestricted form, is probably the most used concordance measure in practice [5]. However, it was affected by censoring, which is a finding noted by others [16]. Specifically, C H tended to increase for moderate to high levels of censoring, which is not an uncommon scenario with medical data, and is therefore likely to give an overoptimistic view of a prediction model's discriminatory ability. Therefore, it cannot be recommended in such scenarios. In contrast, both C U and C GH were reasonably stable in the presence of censoring. C GH was the less variable of the two measures as a consequence of it being model-based [16]. The restricted versions of C H and C U were little affected by censoring but care needs to be taken when selecting the time-point to ensure that the time period contains a reasonable number of events.
The remaining discrimination measure D has an appealing interpretation as it can be communicated as a (log) relative risk between low and high risk groups of patients. It requires that the proportional hazards assumption holds and that the prognostic index is normally distributed. As with R 2 D it may be affected by censoring if the prognostic index is skewed [27]. The sole calibration measure under investigation, the calibration slope, was robust to censoring. It assumes that the proportional hazards assumption holds although more general approaches are described by van Houwelingen [33,38].

Conclusions
Harrell's C H is routinely reported to assess discrimination when survival prediction models are validated [5]. However, based on our simulation results, we recommend that C U is used instead to quantify concordance when there are moderate levels of censoring. Alternatively, C GH could be considered, especially if censoring is very high, but we suggest that the prediction model is re-calibrated first. The restricted version of C H may also be used provided that the time-point is chosen carefully. We also recommend that D is routinely reported to assess discrimination since it has an appealing interpretation, although the distribution of the prognostic index would need to be checked for normality. 'Overall performance' measures are perhaps under used in practice. Our recommendation would be to use any of the predictive accuracy measures and provide the corresponding predictive accuracy curves. In particular, R 2 SH and R 2 S have relatively low variability. The calibration slope is a useful measure of calibration and recommended to report routinely. In addition, one could investigate calibration graphically, for example by comparing observed and predicted survival curves for groups of patients [6]. Finally, we also recommend to investigate the characteristics of the validation data such as the level of censoring and the distribution of the prognostic index derived in the validation setting before choosing the performance measures.