Boosting the concordance index for survival data - a unified framework to derive and evaluate biomarker combinations

Introduction

Recent technological developments in the fields of genomics and biomedical research have led to the discovery of large numbers of gene signatures for the prediction of clinical survival outcomes. In cancer research, for example, gene expression signatures are nowadays used to predict the time to occurrence of metastases as well as the time to progression and overall patient survival . While the importance of molecular data in clinical and epidemiological research is expected to grow considerably in the next years , the detection of clinically useful gene signatures remains a challenging problem for bioinformaticians and biostatisticians, especially when the outcome is a survival time.

After normalization and data pre-processing, the development of a new gene signature usually comprises three methodological tasks:

Task 1: Select a subset of genes that is associated with the clinical outcome.

Task 2: Derive a marker signature by finding the “optimal” combination of the selected genes.

Task 3: Evaluate the prediction accuracy of the optimal combination using future or external data.

Task 1, the selection of a clinically relevant subset of genes, is often addressed by calculating scores to rank the univariate association between the survival outcome and each of the genes . In a subsequent step, the genes with the strongest associations are selected to be included in the gene signature.

Task 2, the derivation of an optimal combination of the selected genes, is usually fulfilled by forming linear combinations of gene expression levels based on Cox regression. Due to multicollinearity problems and the high dimensionality of molecular data, a direct optimization of the Cox partial likelihood is often unfeasible . Consequently, marker combinations are often derived by combining coefficients of univariate Cox regression models , or by applying regularized Cox regression techniques (such as the Lasso or ridge-penalized regression ).

Task 3, the evaluation of prediction accuracy, is considered to be a challenging problem in survival analysis. This is because traditional performance measures for continuous outcomes (such as the mean squared error) are no longer applicable in the presence of censoring. In the literature, several approaches to address this problem exist (see, e.g., for an overview). In this article, we focus on the concordance index for time-to-event data ($C$-index ), which has become a widely used measure of the performance of biomarkers in survival studies . Briefly spoken, the $C$-index can be interpreted as the probability that a patient with a small survival time is associated with a high value of a biomarker combination (and vice versa). Consequently, it measures the concordance between the rankings of the survival times and the biomarker values and therefore the ability of a biomarker to discriminate between patients with small survival times and patients with large survival times. This strategy is especially helpful if the aim is to subdivide patients into groups with good or poor prognosis (as applied in many articles in the medical literature, e.g., ). By definition, the $C$-index has the same scale as the classical area under the curve (AUC) in binary classification: While prediction rules based on random guessing yield $C=0.5$, a perfectly discriminating biomarker combination leads to $C=1$.

Interestingly, the derivation of new gene signatures for survival outcomes via Tasks 1–3 is often addressed by combining completely different methodological approaches and estimation techniques. For example, the estimation of biomarker combinations is usually based on Cox regression and is hence carried out via the optimization of a partial likelihood criterion. On the other hand, the resulting combinations are often evaluated by using the $C$-index which has its roots in the receiver operating characteristics (ROC) methodology. This methodological inconsistency is also problematic from a practical point of view, as the marker combination that optimizes the partial log likelihood criterion is not necessarily the one that optimizes the $C$-index. In other words, if the $C$-index and therefore the discriminatory power is the evaluation criterion of interest, it may be suboptimal to use a likelihood-based criterion to optimize the marker combination. This issue is, of course, not only problematic in survival analysis but also in regression and classification. A theoretical discussion on the differences between performance measures for binary classification can, e.g., be found in .

To overcome the aforementioned inconsistencies, we propose a unified framework for survival analysis that is based on the same statistical methodology for gene selection (Task 1), derivation of an optimal biomarker combination (Task 2) and the evaluation of a new gene signature (Task 3). As will be demonstrated, all three tasks can be addressed by using the concordance index for time-to-event data as performance criterion. While the $C$-index has already been proposed for gene selection (Task 1) and the evaluation of prediction accuracy (Task 3) , the main contribution of this article is a new estimation technique that addresses the development of optimal combinations of genes (Task 2). To achieve this goal, we propose a method for finding gene combinations that is based on the gradient boosting framework . As will be shown, it is possible to use boosting to derive prediction-optimized gene combinations via direct optimization of the $C$-index. Because this new approach allows for using the $C$-index to address all three tasks, the proposed method leads to a consistent framework for the derivation of gene signatures in biomarker studies where the $C$-index is the main performance criterion.

Methods

A prediction rule for $T$ will be formed as a linear combination

Concordance index

Our proposed framework to derive and evaluate biomarker combinations is based on the concordance index (“$C$-index”) which is a general discrimination measure for the evaluation of prediction models . It can be applied to continuous, ordinal and dichotomous outcomes . For time-to-event outcomes, the $C$-index is defined as

During the last decades, the $C$-index has gained enormous popularity in biomedical research; for example, searching for the terms “concordance index” and “c-index” in PubMed resulted in 1156 articles by the time of writing this article. Generally, a value of $C$ close to 1 indicates that the marker $\eta$ is close to a perfect discriminatory power, while a marker that does not perform better than chance results in a value of 0.5. For example, the famous Gail model for the prediction of breast cancer is estimated to yield a value of $C=0.67$ .

Being a flexible discrimination measure, the $C$-index is especially useful for selecting and ranking genes from a pre-processed set of high-dimensional gene expression data (Task 1 described in the Introduction). In other words, Task 1 can be addressed by computing the $C$-index (and hence the marginal discriminatory power) for each individual gene or biomarker, where only those genes with the highest $C$-index are incorporated into the yet-to-derive optimal combination (Task 2). Although there exist various other ways to rank genes and select the most influential ones, the $C$-index has been demonstrated to be especially advantageous for this task .

An estimator of the $C$-index for survival data is given by

Boosting the concordance index

The core of our proposed framework to address Tasks 1 – 3 is the derivation of a prediction-optimized linear combination of genes that is optimal w.r.t. to the $C$-index for time-to-event data. Our approach will be based on a component-wise gradient boosting algorithm that uses the $C$-index as optimization criterion.

Gradient boosting algorithms are generally based on a loss function $\rho(T,\eta)$ that is assumed to be differentiable with respect to the predictor $\eta\equiv\eta(X)$. The aim is then to estimate the “optimal” prediction function

by using gradient descent techniques. As the theoretical mean in (6) is usually unknown in practice, gradient boosting algorithms minimize the empirical risk $\mathcal{R}:=\sum_{i=1}^{n}\rho(t_{i},\eta(x_{i}))$ over $\eta$ instead.

Setting $\mathcal{R}=-\widehat{C}_{\text{Uno}}(T,\eta)$, however, is unfeasible because $\widehat{C}_{\text{Uno}}(T,\eta)$ is not differentiable with respect to $\eta_{i}$ and can therefore not be used in a gradient boosting algorithm. To solve this problem, we follow the approach of Ma and Huang and approximate the indicator function in (7) by the sigmoid function $K(u)=1/(1+\exp(-u/\sigma))$. Here, $\sigma$ is a tuning parameter that controls the smoothness of the approximation (details on the choice of $\sigma$ will be given in the Numerical Results section). Replacing the indicator function in (7) by its smoothed version results in the smoothed empirical risk function

By definition, the smoothed empirical risk $-\widehat{C}_{\text{smooth}}(T,\eta)$ is differentiable with respect to the predictor $\eta_{i}$. Its derivative is given by

In the next step of the gradient boosting algorithm, the derivative in (10) is iteratively fitted to a set of base-learners. Typically, an individual base-learner (simple regression tool, e.g., a tree or a simple linear model) is specified for each marker. To ensure that the estimate of the optimal predictor $\eta^{*}$ is a linear combination of the components of $X$, we will apply simple linear models as base-learners (cf. ). In other words, each base-learner is a simple linear model with one component of $X$ as input variable. Consequently, there are $p$ base-learners, which will be denoted by $b_{l}$, $l=1,\dots,p$. Each base-learner refers to one component of $X$ and therefore to one marker (or gene).

The component-wise gradient boosting algorithm for the optimization of the smoothed $C$-index is then given as follows:

Initialize the estimate of the marker combination $\hat{\eta}^{}$ with offset values. For example, set $\hat{\eta}^{}=\mathbf{0}$, leading to $\hat{\beta}_{l}^{}=0$ for all components $l=1,\dots,p$. Choose a sufficiently large maximum number of iterations $m_{\text{stop}}$ and set the iteration counter $m$ to 1.

Compute the negative gradient vector by using formula (10) and evaluate it at the marker combination $\hat{\eta}^{[m-1]}$ of the previous iteration:

Fit the negative gradient vector $U^{[m]}$ separately to each of the components of $X$ via the base-learners $b_{l}(\cdot)$:

Select the component $l^{*}$ that best fits the negative gradient vector according to the least squares criterion, i.e., select the base-learner $b_{l^{*}}$ defined by

Update the marker combination $\hat{\eta}$ for this component:

where sl is a small step length $(0<\text{sl}\ll 1)$. For example, if $\text{sl}=0.1$, only 10% of the fit of the base-learner is added to the current marker. This procedure shrinks the effect estimates towards zero, effectively increasing the numerical stability of the update step .

As only the base learner $\hat{b}_{l^{*}}$ was selected, only the effect of component $l^{*}$ is updated ($\hat{\beta}_{l^{*}}^{[m]}=\hat{\beta}_{l^{*}}^{[m-1]}+\text{sl}\cdot\hat{b}_{l^{*}}^{[m]}(x_{l^{*}})$) while all other effects stay constant ($\hat{\beta}_{l}^{[m]}=\hat{\beta}_{l}^{[m-1]}$ for $l\neq l^{*}$).

Stop if $m=m_{\text{stop}}$. Else increase $m$ by one and go back to step (2).

By construction, the proposed boosting algorithm automatically estimates the optimal linear biomarker combination that maximizes the smoothed $C$-index. The principle behind the proposed algorithm is to minimize the empirical risk $\mathcal{R}=-\widehat{C}_{\text{smooth}}(T,\eta)$ by using gradient descent in function space, where the function space is spanned by the base-learners $b_{l}$, $l=1,...,p$. In other words, the algorithm iteratively descents the empirical risk by updating $\hat{\eta}^{[m]}$ via the best fitting base-learner. Because the base-learners are simple linear models (each containing only one possible biomarker as predictor variable) and because the update in step (5) of the algorithm is additive, the final solution $\hat{\eta}^{[m_{\text{stop}}]}$ effectively becomes a linear combination of these markers.

The two main tuning parameters of gradient boosting algorithms are the stopping iteration $m_{\text{stop}}$ and the step length sl. In the literature it has been argued that the choice of the step length is of minor importance for the performance of boosting algorithms . Generally, a larger step length leads to faster convergence of the algorithm. However, it also increases the risk of overshooting near the minimum of $\mathcal{R}$. In the following sections we will use a fixed step-length of $\text{sl}=0.1$, which is a common recommendation in the literature on gradient boosting (and which is also the default value in the R package mboost ). The stopping iteration $m_{\text{stop}}$ is considered to be the most important tuning parameter of boosting algorithms . The optimal value of $m_{\text{stop}}$ is usually determined by using cross-validation techniques . Small values of $m_{\text{stop}}$ reduce the complexity of the resulting linear combination $\hat{\eta}^{[m_{\text{stop}}]}$ and avoid overfitting via shrinking the effect estimates. In case of boosting the $C$-index, however, overfitting is less problematic as the predictive performance of $\eta$ is not related to the actual size of the coefficients but to the concordance of the rankings between marker values and the observed survival times. As a result, the stopping iteration $m_{\text{stop}}$ in this specific case is less relevant and can be also specified by a fixed large value (e.g., $m_{\text{stop}}=50000$).

Regarding the boosting algorithm for the smoothed $C$-index, an additional tuning parameter is given by the smoothing parameter $\sigma$. While too large values of $\sigma$ will lead to a poor approximation of the indicator functions in (7), too small values of $\sigma$ might overfit the data (and might therefore result in a decreased prediction accuracy). Details on how to best choose the value of $\sigma$ will be given in the Numerical Results section.

The boosting algorithm presented above is implemented in the add-on software package mboost of the open source statistical programming environment R . The specification of the new Cindex() family and a short description of how to apply the algorithm in practice are given in the Appendix.

Evaluation

After having applied the $C$-index to select the most influential genes (Task 1), and after having used the proposed boosting algorithm to combine the selected genes (Task 2), a final challenge is to evaluate the prediction accuracy of the resulting gene combination (Task 3). Since the $C$-index was used for Tasks 1 and 2, it is also a natural criterion to evaluate the derived marker combination in Task 3. As argued before, it is advantageous from both a methodological perspective as well as from a practical one to use the same criterion for estimation and evaluation of a biomarker combination.

When it comes to the $C$-index, two additional points have to be taken into consideration: First, as the task is to obtain the most precise estimation for the discriminatory power, it is no longer necessary to use the smoothed version $\widehat{C}_{\text{smooth}}$ (which was included for numerical reasons in the boosting algorithm) for evaluation. Consequently, we propose to apply the original estimator $\widehat{C}_{\text{Uno}}$ for evaluating biomarker combinations in Task 3. Second, when applying the estimator $\widehat{C}_{\text{Uno}}$ to the observations in a test sample, a natural question is how to calculate the Kaplan-Meier estimator $\hat{G}_{n}^{L}(t)$ of the unconditional survival function of $T_{\text{cens}}$. In principle, there are three possibilities for the calculation of $\hat{G}_{n}^{L}(t)$: The Kaplan-Meier estimator can be computed from either the test or from the training data, or, alternatively, from the combined data set containing all observations in the learning and test samples. Following the principle that all estimation steps should be carried out prior to Task 3, we will base computation of the Kaplan-Meier estimator on the learning data.

Numerical results

We first investigated the performance of our approach based on simulated data. The aim of our simulation study was:

To analyze if the proposed framework is able to select a small amount of informative markers from a much larger set of candidate variables.

To check if gradient boosting is able to derive the optimal combination $\eta$ of the selected markers, and to compare its performance to competing Cox-based estimation schemes.

To investigate the effect of the smoothing parameter $\sigma$ that controls the smoothness inside the sigmoid function, as well as potential effects of the sample size and the censoring rate on the performance of our approach.

The simulated survival times are generated via a log-logistic distribution for accelerated failure time (AFT) models . They are based on the model equation $\log(T)=\mu+\phi W$, where $T$ is the survival time, $\mu$ and $\phi$ are location and scale parameters, and $W$ is a noise variable, following a standard logistic distribution. As a result, the density function for realizations $t_{i}$ can be written as

For Task 1, we first pre-selected a subset of $p^{*}$ predictors from the $p=1000$ available markers. We ranked the predictors based on their individual values of $\hat{C}_{\text{Uno}}$ and included only the $p^{*}=\{5,10,30\}$ best-performing markers in the boosting algorithm. The results suggest that the $C$-index is clearly able to identify markers that are truly related to the outcome: Although all predictors had a relatively high pairwise correlation ($\rho=0.5$), the four informative markers had a selection probability of 98.5% for $p^{*}=5$ (99% for $p^{*}=10$ and 99.5% for $p^{*}=30$).

To find the optimal combination $\eta$ of the pre-selected markers (Task 2), we applied the proposed boosting approach on training samples with size $n=100$. The resulting coefficients for $p^{*}=5$ and smoothing parameter $\sigma=0.1$ are presented in Figure 1. The boosting algorithm seems to be able to derive the optimal combination of the pre-selected markers, as the structure displayed by the coefficients is essentially the same as the one of the underlying true combination $\eta_{\mu}$. The discriminatory power of the resulting biomarker does not depend on the absolute size of the coefficients: As the $C$-index is based solely on the concordance between biomarker and survival time, what matters in practice is the relative size of the coefficients. As can be seen from Figure 1, the estimated positive effect for $x_{1}$ is larger than the one for $x_{2}$. On the other hand, the negative effect of $x_{4}$ is correctly estimated to be more pronounced than the the one of $x_{3}$. The coefficient of the falsely selected marker is on average close to zero.

In a third step, we evaluated the performance of the resulting optimized marker combinations (Task 3) on separate test samples. The resulting estimates $\hat{C}_{\text{Uno}}$ for different simulation settings are presented in Table 1. The highest discriminatory power (median $\hat{C}_{\text{Uno}}=0.763$, range = 0.559–0.819) can be observed for $p^{*}=5$, which is closest to the true number of informative markers. We compared the performance of our proposed algorithm to penalized Cox regression approaches such as Cox-Lasso and Cox regression with ridge-penalization – see Figure 2. The proposed boosting approach clearly outperforms the competing estimation schemes, supporting our view that applying traditional Cox regression might be suboptimal if the discriminatory power is the performance criterion of interest. We additionally computed the optimal $C$-index resulting from the true marker combination $\eta_{\mu}$ with known coefficients. The values of the true $C$-index on the test samples are on average only slightly better than the ones of boosting the concordance index (median $\hat{C}_{\text{Uno}}=0.778$ – see Table 1).

To evaluate the possible effects of different sample sizes and censoring rates we modified the mean censoring time leading to approximate censoring rates of 30% and 70% and generated training samples of size $n=\{50,200,500\}$. Results are included in Table 1. As expected, the $C$-index resulting from our framework increases as censoring rates become small (median $\hat{C}_{\text{Uno}}=0.820$, range = 0.736–0.858) and decreases in settings with a large proportion of censored observations (median $\hat{C}_{\text{Uno}}=0.668$, range = 0.421–0.776). However, the same effect can be observed for the true $C$-index resulting from the true marker combination $\eta_{\mu}$ (30% censoring $\hat{C}_{\text{Uno}}=0.830$, 70% censoring $\hat{C}_{\text{Uno}}=0.690$). For larger training samples, the variance of the coefficient estimates decreases (see Figure 1). As a result, the discriminatory power resulting from our boosting algorithm improves (for $n=500$, median $\hat{C}_{\text{Uno}}=0.778$, range = 0.614–0.818) and gets nearly as good as the true $C$-index ($\hat{C}_{\text{Uno}}=0.781$). This finding further confirms the ability of our approach to find the most optimal marker combination possible – see Figure 2. Note that also the true $C$-index differs slightly between the different sample sizes, as the training sample enters in $\hat{C}_{\text{Uno}}$ via the Kaplan-Meier estimator $\hat{G}_{n}^{L}(t)$.

To investigate the effect of the smoothing parameter inside the sigmoid function, we additionally applied our boosting procedure for every simulation setting with different values of $\sigma$. The resulting estimates $\hat{C}_{\text{Uno}}$ are presented in Table 2. Compared to the effects of the sample size or the number of pre-selected markers $p^{*}$, the smoothing parameter $\sigma$ only seems to have a minor effect on the performance of our algorithm. In light of these empirical results, we recommend to apply a fixed small value (e.g., $\sigma=0.1$, which is also the default value in the Cindex() family for the mboost package – see the Appendix).

For both approaches to fit penalized Cox regression (Cox lasso, Cox ridge), we applied the R add-on package penalized . In order to evaluate $\hat{C}_{\text{Uno}}$, we used the UnoC() function implemented in the survAUC package .

Applications to predict the time to distant metastases

In the next step, we further analyzed the performance of our gradient boosting algorithm in two applications to estimate and evaluate the optimal combination of pre-selected biomarkers. All markers are used to predict the time to distant metastases in breast cancer patients. As in the simulation study, we compared the results of our proposed algorithm to Cox regression with lasso and ridge penalization. Additionally, we considered four competing boosting approaches for survival analysis, which do not directly optimize the $C$-index. The first is classical Cox regression, estimated via gradient boosting, while the other three are parametric accelerated failure-time (AFT) models assuming a Weibull, log-normal or log-logistic distribution . For all boosting approaches (Weibull AFT boosting, loglog-AFT boosting and Cox boosting) we used the corresponding pre-implemented functions of the mboost package. To ensure comparability, we used the same linear base-learners as described above for all boosting approaches.

Desmedt et al. collected a data set of 196 node-negative breast cancer patients to validate a 76-gene expression signature developed by Wang et al. . The signature, which is based on Affymetrix microarrays, was developed separately for estrogen-receptor (ER) positive patients (60 genes) and ER-negative patients (16 genes). In addition to the expression levels of the 76 genes, four clinical predictor variables were considered (tumor size, estrogen receptor (ER) status, grade of the tumor and patient age). The data are publicly available on GEO (http://www.ncbi.nlm.nih.gov/geo, accession number GSE 7390).

Similar to Wang et al. , we used the time from diagnosis to distant metastases as primary outcome and considered the 76 genes together with the four clinical predictors. Observed metastasis-free survival ranged from 125 days to 3652 days, with 79.08% of the survival times being censored.

The main results of our analysis are presented in Figure 3. As expected, the unified framework to estimate and evaluate the optimal marker signature based on the $C$-index is not only methodologically more consistent than the Cox and AFT approaches, but also leads to to marker signatures that show a higher discriminatory power on external or future data (median $\hat{C}_{\text{Uno}}=0.736$, range = 0.467–0.854). As discussed in the methodological section, it is crucial to evaluate the discriminatory power on external data: the estimated $C$-index on the training sample was more than 35% higher (median $\hat{C}_{\text{Uno}}=0.986$) and hence extremely over-optimistic .

Considering the interpretation of the resulting coefficient estimates for the clinical predictors, it is crucial to note that boosting methods for the $C$-index and the AFT models yield biomarker combinations $\eta^{*}$ where larger values indicate longer predicted survival times. On the other hand, classical Cox regression models rely on the hazard; higher values are hence associated with smaller survival times. If this is taken into account, results from the different approaches were in fact very similar. Both age of the patients and size of the tumor had a negative effect on the time to recurrence for all seven approaches. The same holds true for the tumor grade poor differentiation which resulted in a negative effect compared to good differentiation and intermediate differentiation. A positive ER status, on the other hand, was associated with a larger metastasis-free survival in all approaches. Regarding the coefficients of the 76 genes, results from our approach to boost the $C$-index were highly correlated to the ones of the other four boosting approaches (which rely on the same base-learners) – absolute correlation coefficients computed from the 100 subsamples ranged from 0.77 to 0.99. Also coefficients resulting from the popular ridge-penalized Cox regression were highly correlated with the ones from our approach – absolute correlation coefficients ranged from 0.47 to 0.84.

Breast cancer data by van de Vijver et al.

The second data set consists of $144$ lymph node positive breast cancer patients that was collected by the Netherlands Cancer Institute . The data set, which is publicly available as part of the R add-on package penalized , was used by van de Vijver et al. to validate a 70-gene signature for metastasis-free survival after surgery developed by van’t Veer et al. . In addition to the expression levels of the 70 genes, the data set contains five clinical predictor variables (tumor diameter, number of affected lymph nodes, ER status, grade of the tumor and patient age). Observed metastasis-free survival times ranged from $0.055$ months to $17.660$ months, with $67\%$ of the survival times being censored.

Resulting values of the $C$-index of the new approach and the six considered competitors are presented in Figure 3. The improvement from applying the proposed unified framework compared to boosting the Cox proportional hazard model or applying ridge-penalized Cox regression was much less pronounced than in the previous data set. However, on average, boosting the $C$-index still led to the best combination of markers regarding the discriminatory power (median $\hat{C}_{\text{Uno}}=0.662$, range = 0.257–0.836). Interestingly, as in the previous data set, the lasso penalized Cox regression was clearly outperformed by the ridge-penalized competitor (which has been suggested for this specific data set by van Houwelingen et al. ). Furthermore, the ridge-penalized approach performed at least as good as the considered boosting approaches (except the new approach to boost the $C$-index). As in the previous data set, we again additionally evaluated the $C$-index on the training sample in order to assess the resulting over-optimism. As expected, the estimated $C$-index on the training sample was extremely biased (median $\hat{C}_{\text{Uno}}=0.973$).

The resulting coefficients for the clinical predictors were again comparable for the seven different approaches. A positive ER status was associated with a larger metastasis-free survival for all seven approaches, the same also holds true for the age of the patient. On the other hand, the size of the tumor, the number of affected lymph nodes and a poor tumor grade resulted for all approaches in a negative effect on the survival time. Regarding the coefficients of the 79 genes, the highest correlation could again be observed for the boosting algorithms: Absolute correlation coefficients obtained from the 100 subsamples ranged from 0.64 to 0.95. Correlation between coefficients resulting from our approach to boost the $C$-index and the ones from ridge-penalized Cox regression was slightly lower, it ranging from 0.30 to 0.82.

Discussion

In this article we have proposed a framework for the development of survival prediction rules that is based on the concordance index for time-to-event data ($C$-index). As the $C$-index is an easy-to-interpret measure of the accuracy of survival predictions (based on clinical or molecular data), it has become an important tool in medical decision making. Generally, the focus of the $C$-index is on measuring the “discriminatory power” of a prediction rule: It quantifies how well the rankings of the survival times and the values of a biomarker (or marker combinations) in a sample agree. In particular, the $C$-index is methodologically different from measures that evaluate how well a prediction rule is “calibrated” (i.e., from measures that quantify “how closely the predicted probabilities agree numerically with the actual outcomes” ). Specifically, prediction rules that are well calibrated do not necessarily have a high discriminatory power (and vice versa).

While several authors have proposed the use of the $C$-index for feature selection and the evaluation of molecular signatures , the main contribution of this paper is a new approach for the derivation of marker combinations that is based directly on the $C$-index. Consequently, when using the proposed method, it is no longer necessary to rely on traditional methods such as Cox regression – which focus on the derivation of well-calibrated prediction rules instead of well-disciminating prediction rules and may therefore be suboptimal when the optimization of the discriminatory power is of main interest.

Conceptually, our approach is in direct line with recent articles by Ma and Huang , Wang and Chang and Schmid et al. who developed a set of algorithms for the optimization of discrimination measures for binary outcomes (such as the area under the curve (AUC) and the partial area under the curve and (PAUC)). Because the $C$-index is in fact a summary measure of a correspondingly defined AUC measure for time-to-event data , our optimization technique relies on similar methodological concepts, such as the application of boosting methods and the use of smoothed indicator functions.

A possible future extension of our approach might be to include the task of selecting the most influential genes in the proposed boosting algorithm. While our simulation study and the breast-cancer examples were based on the pre-selection of genes, the proposed boosting method could also be applied directly to high-dimensional molecular data, so that Tasks 1 and 2 are effectively combined. This can be accomplished by optimizing the stopping iteration so that only a (low-dimensional) subset of the candidate genes is incorporated in the resulting marker combination (“early stopping”, cf. ). Further research is warranted on the issues of early stopping and automated feature selection in the case of boosting the concordance index for survival data.

The results of our empirical analysis suggest that the new approach is competitive with state-of-the-art methods for the derivation of marker combinations. As demonstrated in the Numerical Results section, the resulting marker combinations are not only easy to compute and have a meaningful interpretation but can also lead to a higher discriminatory power of the resulting gene signatures.

Acknowledgments

The authors thank Sergej Potapov for his help with the analysis of the breast cancer data. The work of Matthias Schmid and Andreas Mayr was supported by Deutsche Forschungsgemeinschaft (DFG) (www.dfg.de), grant SCHM 2966/1-1. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

References

Appendix

The concept of boosting was first introduced in the field of machine learning . The basic idea is to boost the accuracy of a relatively weak performing classificator (termed “base-learner”) to a more accurate prediction via iteratively applying the base-learner and aggregating its solutions. Generally, the concept of boosting leads to a drastically improved prediction accuracy compared to a single solution of the base-learner . This basic concept was later adapted to fit statistical regression models in a forward stagewise fashion . One of the main advantages of this approach is the interpretability of the final solution, which is basically the same as in any other statistical model . This can not be achieved with competing machine learning algorithms as Support Vector Machines or Random Forests . Specifically, the boosting approach can be used to develop prediction rules for survival outcomes . Although there exist also likelihood-based approaches for boosting , we will focus here on gradient-based boosting as it is the better fitting approach for boosting the distribution-free $C$-index.

The most flexible implementation of gradient boosting is the mboost add-on package for the Open Source programming environment R . The mboost package contains a large variety of different pre-implemented base-learners and loss functions, that can be combined by the user via different fitting functions. For a tutorial on the how to apply the package for practical data analysis, see .

To apply gradient boosting to optimize linear biomarker combinations w.r.t. the $C$-index in the version of Uno et al. , it is necessary to specify the newly developed Cindex() family inside the glmboost() function.

The Cindex family object includes the sigmoid function $K(u)=1/(1+\exp(-u/\sigma))$ as approximation of the indicator functions in the estimated $C$-index. The sigmoid function is evaluated inside the R functions approxGrad() and approxLoss(), which are part of the Cindex object. The weights

are computed via the internal function compute_weights() for both the empirical risk

(implemented in the risk() function) as well as for the negative gradient

(implemented in the ngradient() function).

Those different functions that define the optimization problem are finally plugged into the mboost specific Family() function to build a new boost$\_$family. Details on how to implement user-specific families in mboost are presented in the Appendix of . The complete Cindex object is then given as follows:

Application

We will briefly demonstrate how to apply the Cindex family in practice to derive the optimal combination of pre-selected biomarkers. We will use the van de Vijver et al. data set of $144$ lymph node positive breast cancer patients that was also considered in the main article. The data set is publicly available as part of the R add-on package penalized . The 70-gene signature for metastasis-free survival after surgery was originally developed by van’t Veer et al. .

We first split the data set in 100 training observations and 44 test observations. To ensure better readability of the code we do not carry out stratified subsampling but just use the first 100 patients as training sample. Model fitting is carried out by the glmboost() function of the mboost package. As linear models are the default base-learners for glmboost(), no additional base-learner has to be specified. As appropriate family object we specify the Cindex family described above.

For evaluating the discriminatory power of the resulting prediction on test data, we use the UnoC() function of the survAUC package . It implements the unbiased estimator $\widehat{C}_{\text{Uno}}$, as proposed by Uno et al. .