How to compare (repeated) cross-validated survival curves

Question

I want to cross-validate survival curves, and i am not sure I am doing it correctly. The problem is as follows:

I have around 50 patients and their survival times and their (say clinical) measurements. The goal is to use the measurements to find a model that can predict whether a new patient will have low or high risk "to survive". Because of the low number of patients, the goal is just to make sure that the (cross-validated, not single) model can "significantly" distinguish between the low and the high risk group.

For now i did the following: I applied a repeated 10-fold CV, i.e. i repeat 50 times the following procedure:

• Remove the n-th fold from the data.
• Apply the feature selection to the 9 remaining folds [*].
• Apply the training method to the 9 remaining folds to obtain a model.
• Use the model to predict on the n-th fold, i.e. I get a risk for each patient in the n-th fold.
• Gather all the predictions over all folds to have an "independent" prediction on all the training set.
• As the model yields a risk for each patient, I can now take the median of all risks to subdivide the all patients into low and high risk groups. These are now (I hope) cross-validated survival curves. Usually I would apply the log-rank test to get a significance, whether the two curves differ.

But as a single cross-validation is not stable, I did repeat the above procedure 50 times. Then I obtain 50 cross-validated survival curves. My main question is:

• How to compare these 50 curves correctly, using some kind of statistics to yield a "significance"?

One way would be not to generate the cross-validated survival after each repeat, but instead just gather all the 50 (cross-validated) risks for each patient and take the mean of it. With these means i could create two new survival curves and take their log-rank statistics. Which sounds somehow reasonable.

Unluckily i did not find many references for this problem ,the closest is the paper from Simon, Subramanian et al, "Using cross-validation to evaluate predictive accuracy of survival risk classifiers based on high-dimensional data", see https://www.ncbi.nlm.nih.gov/pubmed/21324971.

But it seems that this paper tells me that my approach is not correct, and would be just a "measure" of spread (The paper seems not deal with repeated cross-validation, unluckily).

To obtain a significance (H_0=the groups do not differ at all times) for the two survival groups, the paper tells me I should do a permutation test by repeating these two steps say 1000 times:

• Permute the survival data (time+event) and the variables (this amounts to simply permute the survival data, as my table only consists of (time, event, clinical measurements))
• Apply the full procedure above to get a log-rank statistics.

All the log-rank statistics will now yield a distribution, and i can now determine in which percentile my original, not permuted log-rank statistics will be. Using an alpha level of 5%, this would mean 97.5% percent of all permuted log-rank values should be on the left (or right) of my unpermuted value.

My questions:

a) Do I have to use the second method (with does not involve repeating cross-validation but permutation) to get a correct result? Or can I work with the repeated cross-validation procedure I alreay have (And, maybe also quite important: Did I correctly understand the procedure in the paper?)

b) Do i have to use the log-rank statistics in the permutation test? Or can I use Harrell's c-index instead?

c) As the 10-fold CV takes several hours (using random forests), is there any kind of shortcut?

d) If my original approach is not correct, can I re-use the 50 repeats i already computed to get a (good) approximation of the "true" p-value?

e) Is there an R package that can help me in any way with the above permutation test?

Thanks a lot.

[*] It would be better to again apply cross-validation in this step, but I ignore this, hoping it will not make a large difference.

Answer 1

If you want a regression/continuous-risk model, then you might be better off going back to the original R CoxBoost package rather than trying to interface with the mlr package for cross-validation and fight with cross-validated survival curves. The CoxBoost package provides the cv.CoxBoost() function to choose the optimal number of boosting steps and the CoxBoost() function to produce a Cox survival model from all the data in a way that includes feature selection along with their (effectively penalized) regression coefficients, providing a linear predictor of log-hazard. Many contributors to this site would prefer to see such continuous-risk models, as classification schemes typically lose information and hide crucial (often unconscious) assumptions.

Of note, the CoxBoost package also provides an estimPVal() function to provide a permutation-based test of significance of the coefficients in the model you build. This should be adequate for your purpose, if you run with enough permutations. Note that a permutation-based approach was also proposed in the paper by Simon et al that you cite for testing significance of differences between cross-validated survival curves without requiring repeated CV. If you are still worried about the stability of a single CV, note that it's only used in the above approach to select the optimal number of boosting steps, and you could do cv.CoxBoot() 50 times and combine the 50 estimates of mean partial log likelihood against number of steps (mean.logplik) to choose the best number of steps.

To speed up calculations, take advantage of the parallel-processing capabilities provided in cv.CoxBoost() and estimPVal(). I know of no other shortcut, other than to use a cloud computing service with parallel processing if it's too slow on your desktop machine.

You should be including accepted standard clinical features for the disease of interest in addition to your 1700 additional predictors. An advantage of CoxBoost() is that it allows for such "mandatory covariates," and with estimPVal() it provides estimates of whether the additional predictors improve upon what the standard clinical features provide. That's critical in this type of work, as there's a risk that new predictors might seem exciting but then end up only to be proxies for standard clinical features. You want to catch that possibility early on before you waste too much of your career on a dead end.

In your application, however, you may be in some trouble as you only have about 40 events (50 cases, 20% censored). That would limit traditional Cox modeling without overfitting to about 4 predictors, and a best-fitting penalized model might end up with only about 4 effective degrees of freedom when penalization is considered along with the number of included predictors. So your ability to find a new set of predictors that adds to standard clinical features may be limited.

If for some reason you need to work with classification and Kaplan-Meier curves, the paper by Simon et al that you cite is probably state of the art for practical cross-validation of survival data in that context. You could do a lot worse than to follow these recommendations of the Chief of the Biometric Research Branch at the US National Cancer Institute and his colleagues. Your approach for a single run of K-fold CV is close to the one they propose, except that the choice of the number of boosting steps should be repeated within each fold as they recommend. (The general rule in model validation is that all aspects of the model building need to be repeated in each fold of CV or, in other validation approaches, for each bootstrap sample.)

You do seem to understand their approach. Permutation-based tests like they propose aren't that hard to code; in R, the boot() function along with boot.ci() allows for permutation testing as well as bootstrap resampling. The survival data would be the "data" to permute, and the matrix of predictors would be incorporated (without permutation) into the "statistic" function that you write to produce the log-rank value or other statistic you wish to examine. The Simon et al paper does consider time-dependendent ROC curves in cross-validation; as the C-index is equivalent to the area under an ROC curve I suppose that you could use that instead. Harrell, however, recommends elsewhere on this site against using the C-index to compare among survival models, as he finds it not sensitive enough to discriminate. I haven't thought through whether that objection might be overcome in a cross-validation context.