# How to compare (repeated) cross-validated survival curves

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.

• How many predictors are you trying to select among, and how many of the patients had events versus having been censored? Do you have time-to-event information, or just whether an event occurred within a particular time frame (i.e., is this survival modeling or yes/no modeling as with logistic regression)? – EdM Mar 7 '18 at 1:46
• Around 1700 predictors, so its high dimensional and rughly 20% are censored. I have the time-to-event information in days, not intervals, so its survival, not binary. I am using Cox-boosting (surv.cv.CoxBoost in mlr) for this. – complexM Mar 7 '18 at 9:40
• Is there some reason why you have to use boosted random forests for this? It seems like the type of problem that could be handled by an elastic net approach much more expeditiously. If you must use boosting, what is the nature of your base learners? Also, the way you word the question makes it sound like you might be treating this more like a classification (high/low risk) rather than a continuous-risk model. If so, that might not be the most useful approach. – EdM Mar 7 '18 at 15:57
• Yes, there is a simple reason: cv.CoxBoost showed the best performance, and the base learners are Cox models. I tried glmnet (should be elastic nets), they did perform much worse for my dataset. It maybe sounds like a classification problem beacuse I am used to think in these terms, but I explicitly want this to be a regression/risk model, even if the subsequent step is forming two groups. Do you have some suggestions what to change? Because I am quite ok with the results, I just dont know how to "prove" the null-hypothesis that the two resulting groups differ. – complexM Mar 7 '18 at 16:31

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.

• Thank you a lot for your very detailed answer. I really appreciate the time you took to reply to each of my questions. I will need some time to digest the answers (and test the approach you propose) and will write back soon. Thanks again! – complexM Mar 8 '18 at 13:19
• First of all, thanks again for answering and sorry for the long delay on my side. I had to test and think it through. For my purposes I decided to stay with the cross-validation approach, mainly because I am not so much interested in a model and its coefficients, but in its predictions and its generalization ability. For this end I need some internal validation, which is why I kept the cross-validation approach (I could have used bootstrapping as well, but personally I like CV more). – complexM Mar 31 '18 at 10:15
• So I removed the mlr wrapper and the feature selection step to see what happens, and used either CoxBoost or cv.CoxBoost directly, but the results were not as good as with the feature selection step (and the mlr wrapper, which should actually not really alter the results, except that the comparison would not be perfectly fine as the folds of the CV are not the very same). Maybe, as you said, this has directly to do with the dof of the 'data'. – complexM Mar 31 '18 at 10:15
• Thus for now I will follow the approach from Simon et al., and will also use their time-dependend ROC curves instead of the C-index. The latter indeed was not sensitive enough, i.e. from the cross-validated Kaplan-Meier curves (and the subsequent permutation test) a very clear difference was visible, but the C-index did not capture this. Again thanks for all your help! – complexM Mar 31 '18 at 10:15