I am trying to build a two-stage model where I classify my data and split it into two datasets. I then fit a Cox Proportional Hazard model to each dataset and make predictions and can calculate the C-index.

I then want to combine my results in "one-bag" and get the combined results. My problem is however that I cannot figure out a way to "manually" calculate the C-index without using the regression, as I now have two regressions I cannot use the build-in functions in the packages.

Do you guys have any suggestions on how to do this?

I am using the Pysurvival package as I am using the DeepSurv method, using build-in functions in the lifelines library is therefore not an option.


1 Answer 1


The C-index is simply the fraction of pairs of comparable cases in which the predicted and observed order of events agrees. With a single Cox model and at most one event per individual, the pairs of comparable cases are all pairs of cases with events and all event/censored pairs in which the event time is earlier than the censoring time.

With two separate Cox models* you are in a situation similar to a stratified Cox model, with two separate baseline hazards. With separate baseline hazards, the relative hazards for individuals in different strata differ over time. The R survival package thus treats cases in separate strata as non-comparable for concordance estimates.

The procedure most compatible with standard practice would be to combine the two separate concordances as their average, weighted by the numbers of comparable pairs in the two models. That doesn't, however, evaluate anything across your two models. To combine the two separate models, you have to fix either the time of comparison or the survival probability. The following assumes that you have survival probability estimates (or at least survival-ranking estimates) from your two-subset modeling that provides predictions for all cases.

Fixing the survival probability would be conceptually close to concordance; for example, you could compare the orders of median predicted survival times pairwise against the observed survival times. (Mean survival times aren't useful in Cox models in the typical case where survival probabilities don't drop all the way to 0.) With a Cox model you might end up having non-comparable pairs of median predicted times, as some combinations of baseline hazards and covariate values might not get all the way down to a median survival time. It would make sense to do this over a range of survival-probability quantiles, not just the median, to see whether you're getting anything interpretable from combining the models.

Fixing a particular survival time is standard for evaluating model calibration (agreement between predicted and observed survival probabilities), rather than the discrimination between cases evaluated by concordance. Consider whether you might actually be more interested in calibration rather than discrimination.

For an approach to concordance/discrimination related to fixing survival times, you could consider time-dependent receiver operating characteristic (ROC) curves as proposed by Heagerty and Zheng. For their "incident/dynamic" time-dependent ROCs, you construct an ROC curve for your model predictions at each event time, using the case with the event to evaluate sensitivity and the other cases still at risk to evaluate specificity. Integrating appropriately weighted areas under those ROC curves through to a final time of interest gives a measure related to concordance. This method should be possible with your modeling, provided that your model produces some continuous measure of predicted event-time ranking that applies across the 2 models. The code for the risksetROC package in R could provide hints for implementation, as I'm not aware of an implementation in Python.

Having a single unstratified model would overcome the above difficulties, as then the rank order of linear-predictor values provides the rank order of predicted survival times (without having to get predicted survival times).

As you are working in Python, consider the lifelines package for survival work. The package author is making a lot of progress toward providing Python survival-analysis functionality that has long been available in R and its predecessors S/S-Plus. The documents include some succinct but very clear explanations of survival analysis. The package has functions like predict_median and predict_percentile to get predictions from your models, and a survival-probability calibration function. The concordance_index utility function calculates concordance directly from lists of predicted and observed event orders, taking censoring of observed times into account. (You'll have to handle the non-comparability of predicted times yourself in setting up the data.)

*I'm putting aside my concerns about the premature dichotomization that can come from your using a fixed classification to split into two separate datasets at an early stage. Please make sure that you have a good justification for that, rather than building a more flexible model that evaluates all the cases together.

  • $\begingroup$ How do I however get the predicted survival times @EdM? I am using pysurvival, and I can get the predicted survival curves at different times but at which time should I take the predicted survival time? Will that be the area under the survival curve? (I am doing both a dichotomized model and a model which evaluates all cases together.) $\endgroup$
    – andKaae
    Commented May 22, 2021 at 9:05
  • 1
    $\begingroup$ @andKaae I've expanded the answer substantially. The standard practice is to treat strata with separate baseline hazards separately and then average them. What I suggest is the closest I can think of for getting something close to a concordance when trying to combine across two models. Whether this is useful or interpretable will depend a lot on the data and your models. $\endgroup$
    – EdM
    Commented May 22, 2021 at 15:10
  • $\begingroup$ My problem however is that I cannot use lifelines as I am working with non linear CPH models, such as DeepSurv. So I need to find a way to manually calculate the Concordance index by merging the predictions of the two groups. If that cant be done maybe simply just averaging the two c-indexes. $\endgroup$
    – andKaae
    Commented May 30, 2021 at 13:24
  • $\begingroup$ @andKaae I added a paragraph about another approach that might be helpful, based on a weighted average of AUCs from ROCs estimated at each event time. That's a concordance approach related to fixing (a set of) survival times. I'm not aware of a Python implementation. Even if you can't use an R package or the Python lifelines directly, you could use the code to help design your own. $\endgroup$
    – EdM
    Commented May 30, 2021 at 17:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.