Concordance for prognostic models with competing risks

1Clinical decision-making and cost-effectiveness analyses often rely on prognostic models that quantify a subject's absolute risk of a disease event of interest over time. However, study populations increasingly consist of elderly individuals with varying degrees of co-morbidity who are likely to experience one of several disease endpoints other than the endpoint of main interest. As an example, prediction of coronary heart disease (CHD) events in elderly subjects is complicated by the fact that subjects may die from other causes prior to the observation of the disease event of interest.

2It is well known that the naive application of standard survival analysis leads to bias and risk over-estimation if competing risks are present and that specialized methods are needed. A key quantity for medical decision-making in the presence of competing risks is the absolute risk of the event of interest over time as quantified by its (covariate-dependent) cumulative incidence function. Thus, regression models are particularly attractive when they provide subject-specific estimates of the absolute risks based on a set of covariates.

3Several measures for quantifying the accuracy of prognostic models have been adapted from the standard survival setting with only one failure cause to competing risks. Measures include prediction error curves, time-dependent sensitivity, specificity, and area under the receiver operating characteristic (ROC) curve (AUC), and reclassification methods. For survival data, the concordance index is a frequently reported measure of discrimination and we have previously presented a simple adaptation of Harrell's concordance estimator to the competing risks setting.

4In the present paper, we motivate and formally define a cause-specific concordance index in the presence of competing risks. Notably, the proposed concordance index depends only on the cumulative incidence function of the event of interest. We clarify the relation of the concordance to time-dependent AUC measures and discuss a possible alternative definition. We then study estimation of a truncated concordance index in the presence of right-censoring. We introduce an inverse probability of censoring weighted (IPCW) estimator and demonstrate its consistency and asymptotic normality if the working model for the censoring distribution is correctly specified. The empirical bias and mean-square error as well as coverage of asymptotic and bootstrap confidence intervals are examined in a simulation study. Finally, we illustrate the methods for an example of coronary risk prediction in older woman using data from the Rotterdam Study.

5Competing risks data without censoring are given by pairs of data where is the time to the event and is the event type. For the purpose of discussing the definition and estimation of the cause-specific concordance index it is sufficient to assume that there are only two competing events. Thus, for simplicity of presentation we let denote the event of interest and the occurrence of any competing event. In applications, it may be important to model all competing events separately.

6The concordance index is defined for any prognostic score depending on baseline variables which can be used to order subjects with respect to the risk of an event of type 1. For example, could be a single baseline marker or the linear predictor of a regression model for the event of interest derived on a training data set. We assume that higher values of are associated with higher risks of the event of interest.

7To motivate our definition with an example, assume that is the time to death and that a specific treatment were available which prevented death due to the event of interest but would not affect death from other causes. The immediate benefit from such a treatment would be greatest for subject with a high risk of dying from the event of interest early, less for individuals dying from the event of interest late, and negligible for subjects with a low risk of experiencing the event of interest at all (i.e. those likely to die from competing causes). Consequently, for a random pair of subjects and, the first subject would be in greater need of treatment than the second subject if they experienced the event of interest and the second subject experienced the event of interest later ( and ) or not at all. In these cases, the ranking of the risk marker for the pair of subjects is concordant if. Pairs of individuals where both experience the competing event are not comparable as neither of them would be in need of treatment.

8To formally define the concordance probability for the event of interest, we assume an independent test set of i.i.d. realizations of from the joint distribution of the marker and the competing risks outcome and define for any randomly chosen pair of subjects from this distribution. The concordance probability for the competing event,, is defined analogously.

9Define the cumulative incidence function for the event of interest as and the improper random variable as. has a distribution function equal to for and a point mass of at. As an associate editor pointed out, can be written in terms of leading essentially to the standard definition of concordance for survival data:.

10We also note that. Thus, and we can rewrite as According to (2.2), the cause-specific concordance for event 1 depends on and the marginal distribution of the marker but not on the cumulative incidence function of the competing event. This feature is not obvious in formula (2.1) but desirable when the aim is to assess the discriminative ability of a marker for.

11In Appendix B of supplementary material available at Biostatistics online, we illustrate the properties of the concordance probability for a single marker and competing risks outcomes simulated according to cause-specific proportional hazards models with constant baseline hazards. The illustrations suggest that to achieve a high concordance, the marker needs to be strongly associated with an increased cause-specific hazard of the event of interest but only weakly or, even better, reversely associated with the cause-specific hazard of the competing event. This can be explained by the fact that the overall effect of a covariate on the cumulative incidence function of the event of interest depends on both cause-specific baseline hazards and both cause-specific hazard ratios.

12Finally, it is important to discuss modifications of definition (2.1) for tied data. For example, it may happen that. Depending on the application, it may then be sensible to count such pairs with a weight of: To simplify notation, we use definition (2.1) as the basis for our further developments.

13We motivated our definition of concordance with a specific treatment for the event of interest which does not affect the competing event. In this situation, a case subject with has a larger immediate benefit from treatment than a control subject with or as subjects experiencing a competing event have no benefit from treatment at all. However, in other situations, the treatment may affect both event types and then subjects with and with and would not be comparable. Here, it would be more relevant to distinguish cases with from those who haven not had any event up to that time point, i.e. those with. This leads to an alternative definition of concordance: Of note, also depends on the cumulative incidence function of the competing event. Thus, it might be less suitable if the main goal is to assess the relevance of a marker or a prognostic model for predicting the absolute risk of the event of interest alone, and we will not pursue it further. However, could be valuable for assessing joint models for the cumulative incidence of both competing events.

14The proposed concordance measures are closely related to measures of the time-dependent AUC which have been proposed to assess discrimination for competing risks data at a fixed time point. We review these measures in Appendix C of supplementary material available at Biostatistics online, and show that can be written as a weighted average of one proposed time-dependent AUC measure over time. This is in analogy with a similar result for survival analysis without competing risks and supports the use of as a global summary measure of performance.

15We now generalize the concordance index defined in Section 2.1 in two ways. First, we replace the simple prognostic score by a more general prediction model for the risk of event until time, i.e. estimates of, which can be obtained by combining cause-specific hazards models, by fitting a Fine and Gray regression model or by direct binomial regression. Secondly, to include the typical application where individuals have a limited duration of follow-up we define a truncated version of the concordance index. Following and, we define The parameter quantifies the ability of the model to correctly rank events of interest up to time and to discriminate them from competing events. The truncation is necessary to enable estimation of from right-censored data with a limited follow-up duration.

16As before, we can write the truncated concordance (2.2) as a functional of and of the marginal distribution of the predictor values of a pair of individuals: If we introduce notation for the order of the predicted risks at time for a pair of individuals, then

17To indicate the end of follow-up for subject we introduce a subject-specific censoring time. Thus, we observe only where, and. We also use the following notation: The event-free survival probability conditional on the covariate is then given by We allow the censoring distribution to depend on the covariates but assume throughout that is conditionally independent of given. This implies where is the conditional probability of being uncensored at time. Noting we also have for:

18An asymptotically biased estimate of is given by where is an indicator for the order of predicted risks at time. This can be interpreted as the proportion of definitely concordant pairs amongst evaluable pairs, i.e. pairs for which one individual experiences the event of interest and concordance can be decided based on the observed (potentially censored) data.

19This estimate evaluated at the time corresponding to the maximum follow-up duration is a direct adaptation of Harrel's for survival data to the competing risks context and has been previously defined in. While simple, a major problem of this estimator is that by ignoring non-evaluable pairs without any correction, bias is introduced. It is well known that Harrel's depends on the censoring distribution and shares this limitation.

20We derive an IPCW estimate for based on a working model for. Let be a time point where We assume that the model is correctly specified and that for all there exists a uniformly consistent, weakly asymptotically linear estimator sequence with influence function. This implies: For example, we can specify a Cox regression model and use the estimate where is the Breslow estimator of the baseline hazard function and the maximum partial likelihood estimator of the regression coefficients. If the censoring is conditionally independent of the competing risks outcome given the predictors and the Cox model correctly specified, then condition (3.3) is satisfied (see, e.g. ). As an alternative, we could assume that the censoring is independent of the competing risks outcome and the predictors. If this assumption is correct, then the Kaplan–Meier estimate for the censoring distribution satisfies (3.3).

21Based on we construct the weights and define an IPCW estimate of:

22If the working model is correctly specified and is a consistent estimator of then is a consistent estimator of for all. Furthermore, if (3.3) is satisfied, then is asymptotically linear and converges in distribution to a normal random variable with mean 0.

23If the working model is correctly specified and is a consistent estimator of then is a consistent estimator of for all. Furthermore, if (3.3) is satisfied, then is asymptotically linear and converges in distribution to a normal random variable with mean 0.

24A proof of the lemma is given in Appendix D of supplementary material available at Biostatistics online. For the case of independent censoring, supplementary material available at Biostatistics online also presents an explicit formula for the influence function and a consistent estimator of the asymptotic variance. The proposed concordance estimator has been implemented with the function cindex of the R package pec. Example code is provided in Appendix A of supplementary material available at Biostatistics online.

25A simulation study was performed to assess bias and root mean square error (RMSE) of the proposed IPCW estimator and coverage of asymptotic and bootstrap confidence intervals. Simulations were for a single prognostic marker and a parameter-free time-independent model. The covariate was simulated to follow a standard normal distribution.

26Conditional on, uncensored competing risks data was assumed to follow cause-specific Cox-exponential models: This was implemented by simulating latent exponentially distributed event times and and then setting and for and for. We consider two competing risks scenarios. In scenario CR1, we set and. In scenario CR2, we set and. As truncation time points, we used the median and the 75% quantile of the marginal distribution of. The truncation time points and corresponding true values of were determined by simulation based on a large uncensored data set of size.

27Censoring times were drawn under a third Cox-exponential model: The observed time was obtained as the minimum of and and considered as right-censored if. We repeated the simulations for independent censoring and covariate-dependent censoring. For each truncation time point the values of were found by simulation such that the expected proportion of right-censored event times amongst observations with was 25%, 50%, or 75%, respectively. For each of the scenarios, we report results for sample sizes 250 and 1000 averaged across 1000 simulated data sets.

28In each simulated data set, we computed three different estimators of: the naive estimator, the IPCW estimator based on the marginal Kaplan–Meier estimate of the censoring distribution and the IPCW estimator based on a Cox regression model for the censoring distribution.

29Bias and root mean squared errors for chosen as the -quantile of are shown in Table 1. Table 2 shows the associated coverage of percentile bootstrap (all estimators) and asymptotic Wald-type confidence intervals ( only) and contrasts empirical standard errors with average bootstrap and asymptotic standard errors for the estimate. Corresponding results for chosen as the median are shown in Appendix E of supplementary material available at Biostatistics online.

30Based on Tables 1 and 2, we draw the following conclusions. The naive estimator can be biased and this can lead to insufficient coverage.The ICPW estimator can also be biased if the censoring depends on the covariate. In some cases, has a smaller bias but for high rates of censoring it can do worse than even though the censoring depends on the covariate.Coverage of confidence intervals for and was generally close to the nominal 95% except for some scenarios with high censoring rates of 75%. Both average asymptotic and bootstrap standard errors closely resembled the empirical standard errors of.

31The naive estimator can be biased and this can lead to insufficient coverage.

32The ICPW estimator can also be biased if the censoring depends on the covariate. In some cases, has a smaller bias but for high rates of censoring it can do worse than even though the censoring depends on the covariate.

33Coverage of confidence intervals for and was generally close to the nominal 95% except for some scenarios with high censoring rates of 75%. Both average asymptotic and bootstrap standard errors closely resembled the empirical standard errors of.

34Specialist medical societies recommend initiation of preventive treatment for CHD based on a subjects’ predicted 10-year risk for CHD. To accurately predict the absolute risk of CHD in older people, prognostic models for CHD need to account for the competing risk of non-CHD death. In this section, we revisit the example of on coronary risk prediction based on data of elderly women from the Rotterdam Study, a prospective, population-based cohort of elderly subjects living in a suburb area of Rotterdam, the Netherlands.

35We analyzed data from 10 years of follow-up of 4144 women aged between 55 and 90 years who were free of CHD at baseline. During that follow-up period, 389 women experienced a CHD event and 921 women died without prior CHD event. Only 41 women of those event-free had less than 10 years of follow-up. We randomly split the data set into a training data set (2763 women with 249 CHD events) and a validation data set (1381 women with 140 CHD events).

36Using the training set, we estimated the parameters of a Fine–Gray regression model which included the “traditional” baseline risk factors for CHD: age, treatment for high blood pressure (yes versus no), systolic blood pressure (separate slopes depending on whether the subject was on blood pressure treatment or not), diabetes mellitus, log-transformed total cholesterol to HDL cholesterol ratio, and smoking status (current versus never or former smoker). All these risk factors were associated with an increased CHD risk and, except for diabetes, all reached conventional significance. We also investigated the role of age, the strongest predictor variable, as a simple marker for CHD.

37Concordance estimates were obtained for these models in the validation set. The dependence of the censoring distribution on the covariates was investigated with a Cox regression model which yielded no trends and non-significant Wald tests for all variables. Thus, all IPCW estimates of concordance were based on the marginal Kaplan–Meier estimator for the censoring distribution from the validation set.

38The left panel of Figure 1 shows the discrimination ability of the two Fine–Gray models for varying time horizons between 1 and 10 years. The model including all risk factors shows higher discriminative ability compared with the model based on age alone. For both models, concordance estimates stabilize after about 2.5 years of follow-up and remain fairly stable though slightly decreasing. The decrease may occur because earlier events are easier to predict than later events.

39The right panel of Figure 1 shows time-dependent ROC curves at time 10 years for the two models. For this graph, cases were defined as subjects with and, controls as subjects with or. Estimation was also based on IPCW-weighting as implemented in the R package.

40Table 3 shows the estimated concordance for predicting CHD and non-CHD death, respectively, during the 10 years follow-up in the validation data. The Fine and Gray model for non-CHD death used the same covariates as the model for CHD. Age alone is a strong predictor for non-CHD death in this elderly population but the multiple Fine and Gray model did not substantially improve concordance. This is not surprising as most additional covariates are established CHD-specific risk factors which would not be expected to strongly affect non-CHD death (except for other deaths related to the cardiovascular system).

41Table 3 also displays time-dependent AUC measures at years in the validation data using two different definitions of controls. AUCs with controls defined as subjects with or (consistently with our concordance definition) were slightly higher but comparable with the concordance whereas AUCs using only subjects with as controls were substantially higher. This could be explained by the fact that age is a strong predictor of both CHD and non-CHD death which hampers discrimination of CHD events from non-CHD deaths.

42We have presented a formal definition of the concordance probability for prognostic models in the presence of competing risks. Like the concordance probability for survival or binary data, it provides a simple overall numeric measure of discrimination. To deal with right-censored data, we derived an IPCW estimator of the truncated concordance probability and established consistency and asymptotic normality under mild assumptions. Asymptotic properties of the proposed estimator rely on the assumption that the censoring distribution is correctly specified and conditionally independent of the competing risks process given covariates. In many applications, it will be reasonable to assume that the censoring mechanism does not depend on covariates and then the marginal Kaplan–Meier estimate of the censoring distribution can be used. However, if for some reason the design or conduct of a clinical study introduced a dependence between the follow-up duration and covariates which also affect the competing risks process (any component, including the cause-specific hazards of competing events), then it is recommended to use a working regression model for the censoring distribution in order to avoid biased conclusions.

43The fact that we estimate a truncated version of concordance rather than the unconstrained concordance probability should not be seen as a limitation of our approach. Indeed it is impossible to assess the performance of prognostic models beyond the maximum follow-up duration without strong and untestable assumptions. If we assume independent censoring, our estimator is defined if we truncate at any time before or at the largest observed censoring or event time. As for the Kaplan–Meier estimator, the effect of censoring on the variability of the IPCW estimator is increasing with increasing truncation time. However, it is difficult to recommend a general purpose truncation time, in particular because the truncation time point influences the interpretation of the concordance probability. To avoid unstable results in practical applications, we recommend that the analyst develops an appropriate model for the censoring distribution, e.g. the Kaplan–Meier estimator or a Cox model, and then investigates the predicted probabilities of being uncensored at the candidate truncation times. Multiple truncation time points can be evaluated and it can be useful to compare discrimination ability at different truncation time points. A model which is good at discriminating patients with an early failure time (e.g. after surgery) from others may not be good at discriminating subsequent failure times amongst patients who survive a first high risk period.

44As discussed in Appendix C of supplementary material available at Biostatistics online, our approach is related to time-dependent sensitivity, specificity, and ROC curves for competing risks and the concordance can be written as a weighted average of the time-dependent for incident cases and controls defined as observations with or. Thus, the concordance serves as an overall summary of discrimination whereas the time-dependent AUC measures discrimination of the event status at one specific time point.

45We assumed that the prognostic model was derived on an independent training data set and only in this setting are asymptotic or bootstrap confidence intervals for the truncated concordance readily available. Clearly, independent training data are not always available, and even if they are, a joint analysis of all data will be more efficient. However, some form of internal cross-validation is needed to develop and assess a prognostic model with a single data set.

46It is important to emphasize that our definition of concordance assesses a prognostic model for the absolute risk of the event of interest in the presence of competing risks. In line with earlier work, we regard this risk as crucial for medical decision-making in the competing risks setting. However, in many instances explicit consideration of competing events will also be important and modeling the entire competing risks multi-state process will provide further insights. As an example, our illustration of concordance for a single marker (presented in supplementary material available at Biostatistics online) shows that discrimination of prognostic models for the event of interest is hampered if covariates affect both cause-specific hazards with regression coefficients of the same sign, especially if there is a strong association with the competing risk or if the baseline competing hazard is high. This indicates that to achieve high discrimination ability one needs predictors which are only weakly or, even better, reversely associated with the cause-specific hazard of the competing event. Moreover, in settings where all competing events are of similar importance, joint accuracy criteria for the entire competing risks multi-state process are needed and their development is an important area for future research.

47Supplementary material is available at http://biostatistics.oxfordjournals.org.

48M.W. was supported by the Wellcome Trust and the Li Ka Shing Foundation—University of Oxford Global Health Programme. Funding to pay the Open Access publication charges for this article was provided by the Wellcome Trust.