Random forests for survival data: which methods work best and under what conditions?

1 Introduction

Survival analysis – the statistical analysis of event times like deaths or disease recurrence – is important for many studies in medicine, as well as in manufacturing and engineering [1]. Parametric analysis of survival data is well established [1] and can easily achieve the dual goals of estimating a probability of survival beyond some fixed time t and predicting the survival time for a subject with a particular set of covariates. However, as with any other parametric statistical analysis, the validity of the results depends, sometimes critically, on how well the proposed model approximates the true distribution of responses. Random forests for regression and classification [2] are able to detect both direct and interacting signals from covariates, and they can estimate complex regression surfaces from arbitrary response distributions without requiring that a model be specified. These remarkable features have attracted researchers to adapt random forests to handle the complexities imposed by survival data, which are usually censored and highly skewed. Random forests for survival data and related techniques have been developed that produce estimated survival distributions and/or predicted survival times using ensembles of trees that are built and processed in varying ways [3], [4], [5], [6], [7], [8], [9], [10], [11], [12], [13], [14], [15].

Given the growing literature on survival trees and forests, knowing which techniques perform well in terms of point prediction and/or survival function estimation is important. However, few systematic comparisons of the performance of survival trees and forests have been conducted.

An earlier study by Radespiel-Tröger et al. [16] compared seven splitting criteria for building individual survival trees (not forests). They concluded, based on an application to a single real dataset, that splits using an adjusted log-rank statistic produce the best results using the integrated Brier score (IBS) and a version of R ². Later, Shimokawa et al. [17] compared nine older (proposed before 2000) splitting criteria for individual survival trees through a simulation study. They compared four data-generating mechanisms, each at four censoring rates, and evaluated the performance of each splitting criterion using IBS and a form of R ². They found that all but two methods performed comparably well in general at all levels of censoring, and recommended certain splitting criteria over others depending on the data setting (i.e., censoring rate and shape of hazard rate).

Moradian et al. [10] compared seven methods for building survival forests against two standard, non-forest methods using simulation studies and real datasets. The authors tested these against two benchmarks: a Cox model and the Kaplan–Meier (K-M) estimator. They found that their proposed splitting criteria tended to rank among the best performers according to the integrated absolute error and IBS metrics, with much more variable rankings according to the C-index metric. Our simulation study, in contrast, is much more extensive in that it tests various factors that affect survival forest performance.

To our knowledge, to date, no one has undertaken an in-depth comparison of the many alternative forest-based methods for survival analysis. In light of the fact that they are based on a variety of algorithms that focus on producing estimates of different population quantities, these methods may perform quite differently for different tasks. The lack of systematic comparison leaves a knowledge void for both practitioners and researchers. With this paper, we aim to fill part of this gap and to supply both practitioners and statistics and machine learning researchers with advice on which techniques should be considered for the analysis of survival data with particular properties. We systematically investigated various factors that influence survival forest performance (forest construction method, censoring, sample size, distribution of the response, structure of the linear predictor, and presence of correlated or noisy covariates via an extensive simulation study with well-defined performance metrics. In particular, we studied 11 methods that have recently been proposed in the literature and identified 6 top performers. We found that all the factors that we investigated have significant impact on the methods’ relative accuracy of point predictions of survival times and survival function estimates.

We begin in Section 2 with a review of more recently proposed methods for building survival forests. In Section 3, we describe in detail the design of our simulation study as well as the metrics we used to evaluate performance. In Section 4, we discuss our results and recommendations for practitioners. In Section 5, we offer explanations for the observed differences in some methods’ relative performance in certain settings. We conclude with our main contributions, our study’s limitations, and future research opportunities in Section 6.

2 Survival forest methods review

Our review of the literature revealed many methods for building forest-based ensembles for survival data (hereafter “survival forests”). For our comparative study, we focused on more recent methods because they are generally favoured in other papers’ comparisons. We also restricted our focus to methods designed for estimating survival functions and point prediction. Therefore, we excluded methods focused on estimating treatment effects (e.g., [18, 19]). Table 1 lists all such methods proposed since 2000 that we were able to find. The first five methods are based on the forest construction algorithm proposed by Ishwaran et al. [6] and differ only in the splitting criterion used (thus we name the methods based on the splitting criterion). The other methods, proposed in mostly subsequent papers, are characterized by alternative forest construction algorithms that differ in ways beyond the splitting criterion. We describe all these methods below.

The most popular method, random survival forests (RSFs), proposed by Ishwaran et al. [6], adapts the random forest algorithm [2] to survival data. The main algorithmic adaptations in random survival forests are the splitting criterion used to recursively partition the bootstrap samples and the statistic that summarizes the observations in each terminal node of each tree in the ensemble.

The RSF algorithm begins by drawing B bootstrap samples from the training data. The data excluded from each sample represent a test set, the “out-of-bag” (OOB) data, for each sample. A survival tree is then constructed using each bootstrap sample. At each tree node, M covariates are randomly selected as candidates upon which to base a split. The split is made on the candidate variable that results in the greatest survival differences, according to some splitting criterion, between the observations in each child node. Successive splits are then made on each of these child nodes until the tree is grown to full size, ensuring that each terminal node contains no fewer than one unique failure. The observations in each terminal node are then used to compute the Nelson-Aalen estimate of the cumulative hazard function (CHF) for that node. Therefore, for a new individual, each tree produces an estimated CHF. These estimates can then be averaged to obtain an ensemble CHF, which serves as the final RSF prediction for this individual [6].

Ishwaran et al. [6] proposed several splitting criteria for constructing trees within RSFs. The log-rank splitting criterion, which finds the split that maximizes the log-rank statistic, has since become perhaps the standard splitting criterion. Ishwaran et al. [6] proposed a standardized variant of the log-rank splitting criterion, the log-rank score statistic.

Other statistics have been proposed by other authors as splitting criteria for RSFs. Schmid et al. [9] proposed the C-index, hypothesizing that better performance might be achieved when the same method is used for both splitting and for measuring error (C-index is the standard method of quantifying prediction error of RSFs, which we describe in more detail in Section 3.3.1). Due to the suboptimal performance of the log-rank splitting criterion when survival or hazard functions cross, Moradian et al. [10] proposed an L ₁ splitting criterion that is based on the integrated absolute difference between each child node’s estimated survival distributions. Boström et al. [11, 12] used the Brier score – another commonly used error metric (see Section 3.3.2) – as a splitting criterion. We refer to the forests based on the three splitting criteria discussed in this paragraph as RSF-C-index, L ₁, and RSF-Brier.

Random forests have a bias towards variables with many possible splits, as identified by Breiman et al. [20]. To address this issue, Hothorn et al. [4] proposed a new type of decision tree, termed a “conditional inference” tree, which can be amalgamated into conditional inference forests (CIFs). CIF-based trees differ from RSF-based trees in several key ways. First, after randomly selecting M covariates (similarly to the RSF algorithm), a null hypothesis of independence between each candidate variable and the survival outcome is tested using permutation tests [21]. The candidate variable with the strongest association to the survival outcome is chosen for splitting. If there is no statistical evidence of an association between the survival outcome and any of the candidate variables, the node becomes a terminal node. Second, the CIF differs from the RSF not just in the structure of individual trees but also in the aggregation mechanism. Specifically, while RSFs use an equal weighting scheme to compute the ensemble CHF, CIFs adopt a nearest neighbour weighting scheme to aggregate observation times from each tree, yielding an ensemble survival function.

The k-fold RIST algorithm [7] is based on the extremely randomized trees (ERTs) algorithm [22]. To start, the algorithm generates B ERTs. For each node split, a random subset of covariates is selected (like with RSFs); for each of these covariates, a random split point is proposed, and the one that best optimizes a specific criterion is chosen for the actual split. The ERTs are then aggregated to estimate the survival distribution for each censored observation (given the covariates). Using these distributions to impute B sets of random survival times with which to replace the censored survival times, B new datasets are created and B new ERTs are fit. This process of model fitting and imputation is recursively repeated k ≥ 1 times, where k is selected by the user. The survival distribution for a new individual is estimated via the ensemble survival function based on the ERTs in the kth iteration.

The RotSF algorithm produces an ensemble of trees, where each tree is based on a bootstrapped sample and a perturbed version of the matrix of explanatory variables [8]. The perturbations are achieved by applying principal components rotations to randomly selected subsets of explanatory variables separately for each tree. The standard RSF algorithm is then applied to the resulting dataset.

ORSF also transforms the matrix of explanatory variables used for building trees [13]. At each split, Cox proportional hazards models are fit using three different amounts of regularization. These fitted models produce three estimates of the hazard function for each individual. These estimated hazard functions – which are proportional to exponentiated linear combinations of the variables selected for that split – act as three candidate splitting variables. The log-rank splitting criterion is then used in the usual manner to determine whether the node should be split on any of these variables and, if so, where. ORSF requires selecting the tuning parameter, α, which governs the elastic net penalty [13]; Zou and Hastie [23] found that the choice of α = 0.05, indicating an elastic ridge-type penalty, led to the “best” performance according to their error metrics. Like CIF, ORSF produces a weighted K-M estimate of the conditional survival function instead of a CHF, and it uses sampling without replacement instead of bootstrap sampling.

Li and Bradic [15] proposed censored regression forests based on a simple relationship among the conditional distributions of the observed times (the minimum of the survival and censoring times), the survival times, and the censoring times. Their algorithm involves first estimating the conditional distribution of the observed times, given the covariates, using either generalized random forests (CRF-GRF; [14]) or quantile regression forests (CRF-QRF; [5]). They then estimate the conditional distribution of the censoring times given the covariates (via, for example, the K-M estimate). Conditional quantiles of the distribution of survival times are then estimated from the quotient of these two estimated conditional distributions.

3 Methods

In this section, we describe the design of our simulation study, which we used to compare survival forest methods, and the error metrics we used to evaluate each method’s performance.

We use the following notation throughout. Let T and C be the survival time and censoring time, respectively, of an individual. We assume that T and C are conditionally independent given the covariates. Let Y = min(T, C) be the observed time, and let δ = 1 ( T ≤ C ) denote the censoring indicator. Let X = (X ₁, …, X _p ) denote a p-dimensional vector of covariates. Let S(t) and h(t) denote the survival and hazard functions, respectively, of T. We assume in all cases that the covariates affect the conditional distribution of survival times only through a linear predictor, Xβ , where β is a vector of parameters that modifies the effect of X . We describe the linear predictors that we consider in Section 3.1.2.

4 Results

We begin by giving an overview of our results. We present the details in Sections 4.1 and 4.2, followed by formal recommendations in Section 4.3. An even more granular discussion follows in Section 4.4, where we focus on the most influential factors and factor combinations. We present results mainly from the comprehensive study, but we also comment on the pilot study where relevant.

Based on the totality of the results from our studies, all six methods included in our comprehensive study were relatively competitive overall, with generally good performance across most FLCs and rare poor performance. As detailed below, some methods (notably, RotSF and RIST) performed especially well in certain settings, while others (e.g., RSF-log-rank) performed more consistently below average relative to the six top-performing methods.

Our pilot study revealed that RSF-log-rank-score was clearly the worst-performing method and also that CRF-GRF, CRF-QRF, RSF-Brier, and RSF-C-index were relatively uncompetitive both for point prediction and for survival function estimation (see Supplementary Information for relevant plots).

For the ensuing discussion, we separate our recommendations of methods for point prediction and those for survival function estimation based on the rationale provided in Section 3.3; that is, in some clinical settings, practitioners are likely more concerned with predicting a survival time, whereas in other settings, the entire survival function will be of interest.

4.1 Descriptive statistics

In Figure 1, we present boxplots of relative sample mean error values. The boxplots are comparisons of methods’ relative performance across all 144 FLCs in the comprehensive simulation study. Specifically, for each error metric, we computed the sample mean error across the 10 replicate datasets for each method separately within each FLC. We then centred these sample means around the overall sample mean (across methods) for that FLC. Negative values therefore represent better observed performance (smaller average error) for that FLC. Each boxplot consists of 144 relative mean error estimates.

Figure 1:

Relative sample mean absolute loss by method (left); relative sample mean IBS by method.

We emphasize that we present the boxplots to summarize and visualize differences in method performance across FLCs. We conduct formal inference in Section 4.2.

4.1.1 Point prediction

We used the estimated median for point prediction of survival time, and we quantified the accuracy of the predictions using absolute loss and the C-index (as discussed in Section 3.3.1).

In the sample from our comprehensive study, the median performances of all six methods were comparable, while the variability and distributions differed (see left panel of Figure 1). Using a few criteria, we can identify better and worse methods. Specifically, RIST had absolute loss below the mean for at least 3/4 of the FLCs (upper limit of box is below 0). All other methods – especially, CIF, L ₁, and RSF-log-rank – frequently had absolute losses above the mean, with their rather variable upper quartiles entirely above zero, although their medians occur essentially at the mean among methods. The tails at one or both ends of each boxplot are long, which may have several causes. It may relate to relative ease or difficulty that different methods have in producing accurate survival time predictions under different FLCs. It may also reflect the variability of absolute loss, and more generally the difficulty in doing point prediction accurately using any method in any situation compared to survival function estimation – see next section.

Somewhat different results emerged using E = 1 − C (see Figure 12 in Appendix B). Unlike with absolute loss, CIF and ORSF more often performed better than the mean according to E. The other methods’ relative performance was similar using both point prediction error metrics.

Within each method and each FLC, the sample correlations of the two error metrics averaged roughly −0.11 but varied from −0.53 to 0.45. Since the two metrics represent different underlying quantities, divergence in performance is not necessarily concerning or surprising. Simple scenarios involving observed and predicted survival times can be constructed such that the two metrics yield opposite conclusions regarding which method is “best”. For example, the predicted rank could be perfect, yet the absolute differences could be relatively large; or the predicted rank could be terrible, yet the absolute differences could be relatively tiny. Clearly, E cannot be used as a surrogate for absolute loss. Therefore, we caution against using E when the specific goal is minimizing the absolute difference between predicted and true survival times. It is an appropriate measure only when predicting rank or ordering of survival times is the primary interest. See the Supplementary Information for histograms of the sample correlations between the two error metrics.

4.1.2 Survival function estimation

For survival function estimation, in the sample from our comprehensive study, RotSF had estimated mean IBS that was lower than the average for all methods for at least 3/4 of all FLCs (see right panel of Figure 1). Conversely, RSF-log-rank performed worse than average for over 3/4 of all FLCs. The other three methods, CIF, ORSF, and RIST, performed roughly similarly in both median performance and variability, with RIST performing relatively badly in some outlying cases. Recall that we excluded L ₁ (as well as CRF-GRF and CRF-QRF in the pilot study) due to difficulties in computing IBS from their output (see Section 3.3.2).

4.2 Effects of sample size and censoring

In this subsection, we examine methods in more detail for specific combinations of the Censoring and SampleSize factors because they influence method performance and are known (or can be trivially estimated) prior to the analysis of a dataset. The other factors, DGM-Y, DGM-X, and ExtraVars, while possibly important, are hard to determine from the data (though subject-area knowledge may provide some information). We acknowledge that performance differences among the combinations of Censoring and SampleSize levels can be influenced by differences in their effects at the levels of other factors. But we want to present results for guidance in the common case where the analyst is blind to the levels of these factors in their study population. We defer to Section 4.4 discussion of how these more elusive factors affect method performance.

As in Section 4.1, we separate our recommendations for point prediction and survival function estimation. In each subsection, we discuss relative performance differences of the methods for each combination of Censoring and SampleSize, then summarize the findings in Table 3 (in Section 4.3).

Table 3:

The better-than-average (“+”) and worse-than-average (“−”) performing survival forest construction methods for each combination of Censoring (C) and SampleSize (n), for each error metric.

Loss	C	n	CIF	ORSF	RIST-1	RotSF	RSF-log-rank	L1
AL	10	200	−	−	+	+	+	−
	10	1500	−	−	+	−	+	+
	70	200
	70	1500					−
IBS	10	200	−		+	+
	10	1500	−	+		−	−
	70	200		−	−	+	−
	70	1500	+	−	−	+	−

We use the emmeans package on each LMM to extract the estimated marginal mean error for each method at each level of a selected factor or subset of factors. These marginal means are estimated by averaging the estimated conditional mean errors over the levels of all other factors using equal weights. We then use comparison contrasts that represent the difference between each of the marginal means and the overall marginal mean for each combination of Censoring and SampleSize. The figures in Sections 4.2.1 and 4.2.2 represent 95 % simultaneous confidence intervals for these differences. Thus, confidence intervals that are fully below or above 0 represent statistically significant differences from the mean (better and worse performances, respectively) among methods.

The confidence interval plots for the pilot study are available in the Supplementary Information.

4.2.1 Point prediction

4.2.1.1 10 % censoring

At our lowest Censoring level and at n = 200, RIST, RotSF, and RSF-log-rank had mean absolute losses significantly below the average across all methods, while CIF, ORSF, and L ₁ had mean absolute losses significantly above the average (see left panel of Figure 2). Moreover, we confirmed through pairwise comparisons that each method’s mean absolute loss was significantly different from all others (except for ORSF and CIF), so we can rank methods by performance for these factor levels, with RIST best and CIF worst.

Figure 2:

95 % simultaneous CIs for relative mean absolute loss (AL) for all methods at SampleSize = 200 (left) and SampleSize = 1500 (right) at Censoring = 10 %.

At n = 1500 (see right panel of Figure 2), RIST, RSF-log-rank, and L ₁ had mean absolute losses significantly lower than the mean, while ORSF, CIF, and RotSF had mean absolute losses significantly higher than the mean. Once again, we confirmed through pairwise comparisons that RIST had mean absolute loss significantly lower than all other methods, while L ₁ and RSF-log-rank had mean absolute losses significantly lower than the three worse-than-average methods. Interestingly, RotSF performed significantly worse than average at this higher sample size (compared to better at the lower sample size). Conversely, L ₁ performed significantly better than average at the higher sample size and worse at the lower sample size.

4.2.1.2 70 % censoring

At our high Censoring level and at our smaller SampleSize level, comparisons indicated that no method had mean absolute loss that was significantly different from the mean across methods (Figure 3). These results are presumably due to the high variability in performance that is induced primarily by the high censoring rate and, to a lesser extent, the smaller sample size. Accordingly, at our higher SampleSize, no method had mean absolute loss significantly lower than the mean, while only RSF-log-rank performed significantly worse than average. The relative performance differences we identified at the 10 % Censoring level were almost entirely absent at 70 % Censoring. Based on our results, when Censoring is high, SampleSize is large, and the levels of other factors hard to determine, we recommend any of the methods except RSF-log-rank.

Figure 3:

95 % simultaneous CIs for relative mean absolute loss (AL) for all methods at n = 200 (left) and n = 1500 (right) at 70 % censoring.

4.2.1.3 C-index

Use of E = 1 − C as the error metric results in a substantially different picture of the relative performance of the methods. In particular, some methods that performed better than average according to C-index (i.e., CIF, ORSF) performed worse-than-average according to absolute loss, and vice versa (i.e., L ₁, RSF-log-rank). See confidence interval plots in Appendix B.

4.2.2 Survival function estimation

The best- and worst-performing methods differed somewhat for survival function estimation compared to point prediction. Of the five methods compared, all but RSF-log-rank had some combinations of Censoring and SampleSize where their performance was better than average and some combinations where their performance was worse than average.

4.2.2.1 10 % censoring

At our lowest Censoring level and at n = 200, both RIST and RotSF had mean IBS significantly lower than the mean across methods, while CIF and RSF-log-rank had significantly higher mean IBS, and ORSF had an indistinguishably different mean IBS (see left panel of Figure 4). As verified through pairwise comparisons, RIST performed significantly better than RotSF, and RSF-log-rank performed significantly better than CIF.

Figure 4:

95 % simultaneous CIs for relative mean IBS for all methods at n = 200 (left) and n = 1500 (right) at 10 % censoring.

At n = 1500 (see right panel of Figure 4), RIST once again performed significantly better than average and significantly better than all other methods. Like at n = 200, CIF and RSF-log-rank had mean IBS significantly higher than the mean. Unlike at n = 200, RotSF had mean IBS significantly higher than the mean (and higher than that of CIF and RSF-log-rank, as verified through pairwise comparisons). As when absolute loss was used as the response, the lower Censoring level seems to lead to a reversal of the relative performance of RotSF at the two different sample sizes: RotSF appeared to struggle in a relative sense at the higher level of SampleSize when Censoring was at the low level.

4.2.2.2 70 % censoring

At our high Censoring level and at our smaller SampleSize, only RotSF had mean IBS significantly lower than average, while RSF-log-rank, ORSF, and RIST each had mean IBS significantly higher than the mean, and CIF performed indistinguishably differently from the mean (Figure 5). As confirmed through pairwise comparisons, RSF-log-rank performed much worse than all other methods.

At our higher SampleSize, RotSF yet again performed by far the best, achieving a mean IBS significantly lower than all other methods. Both RotSF and CIF performed better than the mean. The other three methods performed significantly worse than the mean, with RSF-log-rank performing by far the worst. Unlike with other high SampleSize FLCs, RotSF performed relatively better at these factor levels.

Figure 5:

95 % simultaneous CIs for relative mean IBS for all methods at n = 200 (left) and n = 1500 (right) at 70 % censoring.

4.4 Additional results: most influential factors and factor combinations

In this subsection, we discuss the most influential three-way interaction effects. We present these results separately because they concern at least one factor that is hard to determine and therefore are of more theoretical than practical interest. To identify these effects, we once again used the ANOVA output (see Appendix D) from each LMM, which indicates that the vast majority of effects are significant at the 5 % level. We also examined interaction plots to identify the more specific settings in which some methods are likely to perform better or worse than others, where “better performance” means a lower error.

Conveniently, the three-way interaction terms with the greatest F-values from the ANOVA ouptut align closely for all three error metrics. The magnitude of some interaction effects was often driven higher by a small subset of methods that performed much worse than the others, such that when we removed these weak methods from consideration, the estimated interaction effects were greatly reduced. In these cases, we discuss the relevant lower-order effects only.

4.4.1 Point prediction

4.4.1.1 DGM.Y:Censoring:Method

Using absolute loss as the response, the most important three-way interaction effect according to the ANOVA table is DGM-Y:Censoring:Method. Interaction plots indicate the presence of this three-way effect (see Figure 6). Most methods appear to contribute to this effect, judging from how the two-way interactions for Y-Weibull look different from the two-interactions for Y-GG. (We omit Y-AFT as the results look similar to those for Y-Weibull.) Focusing just on Y-Weibull, the two-way interaction between Censoring and Method is obvious: at the 70 % Censoring level, no method performs discernibly better or worse than any other, whereas competitive advantages and disadvantages are very clear at the 10 % Censoring level. In contrast, focusing just on Y-GG, differences in performance among methods are much more consistent across Censoring levels.

Figure 6:

Estimated mean log(log(absolute loss)) interaction plot of Method:Censoring plotted for two levels of DGM-Y: Y-Weibull (left) and Y-GG (right).

4.4.1.2 Censoring:SampleSize:Method

Censoring:SampleSize:Method was the second most influential three-way interaction effect on absolute loss. Figures 7 and 8 suggest that the magnitude of this interaction is largely due to RotSF. This finding is consistent with the CI plots in Section 4.1.1, which revealed a reversal in RotSF’s relative performance at our two levels of SampleSize when holding Censoring constant. We refit the associated LMM without this method and observed that the F-values decreased substantially. These results confirm our recommendations in Section 4.3 to avoid RotSF for point prediction when both the sample size is large and censoring is low.

Figure 7:

Estimated mean log(log(absolute loss)) interaction plot of Method:SampleSize plotted for Censoring = 10 % (left) and Censoring = 70 % (right).

Figure 8:

Estimated mean log(log(absolute loss)) interaction plot of Method:Censoring plotted for SampleSize = 200 (left) and Sample Size = 1500 (right).

4.4.2 Survival function estimation

4.4.2.1 Censoring:Method:DGM-Y

When we used IBS as the response, the most important three-way interaction term was Censoring:Method:DGM-Y. The pairs of interaction plots in Figures 9 and 10 show visual evidence of three-way interactions. In particular, the DGM-Y:Method plots at two different Censoring levels suggest a three-way interaction driven primarily by methods’ different performance with Y-GG compared to that with other DGM-Y levels. In particular, compared to their performance for other DGM-Y levels (a) RotSF and CIF have comparatively lower IBS for Y-GG than other methods; (b) RIST and RSF-log-rank have relatively higher IBS; and (c) the average difference in response to Y-GG across all methods is bigger for the higher censoring level. These results underscore that, for survival function estimation, RIST’s relative struggles at our high Censoring level occurred primarily in the Y-GG setting. Moreover, the results exemplify how RotSF’s relative strength at our high Censoring level was consistent across levels of DGM-Y.

Figure 9:

Estimated mean log(IBS) interaction plot of Censoring:Method plotted for Censoring = 10 % (left) and Censoring = 70 % (right).

Figure 10:

Estimated mean log(IBS) interaction plot of DGM-Y:Method plotted for Censoring 10 % (left) and Censoring 70 % (right).

4.4.2.2 SampleSize:Method:DGM-X

The next most important three-way interaction term according to the ANOVA table was SampleSize:Method:DGM-X. The patterns of differences in method performance change across each level of Censoring (see Figure 11). The bulk of the three-way interaction effect appears to be driven primarily by X-piecewise, the DGM-X level whose parameter values change depending on the values of X. Specifically, at the 10 % Censoring level, differences in method performance within each level of DGM-X are small and relatively consistent, whereas at the 70 % Censoring level, differences are larger, especially for X-piecewise.

Figure 11:

Estimated mean log(IBS) interaction plot of Method:DGM.X plotted for Censoring = 10 % (left) and Censoring = 70 % (right).

5 Discussion

Our results raise the obvious question: why do some methods perform relatively better and worse than other methods in certain settings? In this section, we hypothesize that certain algorithmic differences in forest construction methods may explain their empirical performances in our simulation study.

The most common splitting criterion, RSF-log-rank, attained a middling or relatively good performance for point prediction at our low Censoring level, especially at our higher SampleSize level and when hazards did not cross. However, RSF-log-rank struggled with survival function estimation for other combinations of DGM-Y, Censoring, and SampleSize (see Table 8 in Appendix E). Recall that Y-Weibull has proportional hazards; Y-AFT has non-proportional, non-crossing hazards and can have crossing survival functions; and Y-GG has both crossing survival and crossing hazard functions. The log-rank test may have low power if the hazard functions are not proportional, so RSF-log-rank may lose power to detect survival function differences outside this setting, as was suggested by Moradian et al. [10]. For survival function estimation, RSF-log-rank did perform consistently below average for Y-GG (where hazard functions may cross). However, for point prediction, RSF-log-rank sometimes performed above average for Y-AFT. Moreover, for survival function estimation, even for Y-Weibull, where hazard functions do not cross, RSF-log-rank tended to perform below average. This latter finding suggests that some other methods (notably RotSF and RIST) have competitive advantages due to reasons unrelated to the power of the log-rank test.

In light of concerns about the performance of RSF-log-rank when hazard or survival functions cross, Moradian et al. [10] proposed L ₁ as an alternative splitting criterion. But our results showed that, for point prediction, L ₁ did not out-perform RSF-log-rank at any DGM-Y level for any combination of Censoring and SampleSize, even those where the hazard or survival functions cross. An explanation may be as follows. Like all the forest methods we are considering, L ₁ determines an optimal split of each node by assuming that the observations within each daughter node are identically distributed. However, the observations are likely far from identically distributed at nodes closer to the tops of trees. In such settings, the performance of the test based on L ₁ relative to that of the log rank test is unknown. Moreover, even if the hazard or survival functions cross for certain values of the covariates, the estimated hazard or survival functions for a given pair of daughter nodes may not cross. In other words, splitting criteria that are based on tests designed to handle crossing hazard or survival functions may not be necessary for identifying optimal splits in some locations within a tree. In addition, L ₁’s overall middling or subpar performance for Y-AFT and Y-GG may reflect strengths of other methods unrelated to the splitting criterion.

Referring back to Table 3, CIF performed worse than average at our low Censoring level using both error metrics at both levels of SampleSize. A major difference in the CIF algorithm compared to RSF-log-rank is that splitting terminates earlier if none of the candidate variables is significantly associated with the survival outcome. Earlier termination leads to shallower tree depth, which may lead to greater bias in both point predictors and survival function estimators. The effects of early termination may be more prominent when censoring is low because, in this case, the other methods may be able to form substantially larger and less biased trees using the large proportion of observed survival times.

In contrast, RIST performed better than average at our low Censoring level at both levels of SampleSize for point prediction and at our lower SampleSize level for survival function estimation. Recall that RIST recursively imputes censored observations using information from their estimated conditional survival distributions. One advantage of this procedure is that deeper trees can be grown due to a greater number of available survival times by, in effect, converting some censoring times into survival times (recall that terminal nodes generally must have a minimum of one survival time). In this way, the bias of the estimated quantities of interest may be reduced. At our low Censoring level, this imputation procedure does indeed appear to extract useful additional information from the data. At our high Censoring level, RIST either had no competitive advantage (using absolute loss) or performed worse than average (using IBS), suggesting that the imputed times may increase the bias or variance of the estimated survival functions. However, a closer inspection of the FLCs where RIST performed relatively poorly reveals that all were associated with Y-GG and that, for survival function estimation, RIST generally performed better than average at our high Censoring level for Y-Weibull and Y-AFT. So crossing hazards appear to pose a greater challenge to RIST than do high censoring rates, and the likely benefit of greater tree depth does not compensate for the challenges posed by Y-GG.

RotSF was unique in that it tended to perform better than average at our low SampleSize level and worse than average at our high SampleSize level (while still performing better in an absolute sense at the latter). For survival function estimation, RotSF also performed relatively better than all but one other method at our high Censoring level. That is, when information is more plentiful – i.e., a low censoring rate and a larger sample size – RotSF appears to do relatively worse, but with information scarcity, its algorithm appears to harness that information more efficiently. One possible explanation for these findings is as follows. RotSF begins by finding principal components of random subsets of the covariates. If the number of influential variables is sufficiently large, then their relationship with the response may be complex enough that it is best modelled with splits based on linear combinations of variables (e.g., principal components) rather than the variables themselves. In these cases, the most important principal components will explain more variance in the response than do the most important covariates, thus allowing for the construction of forests with better predictive performance. In an information-scarce setting (small sample size, high censoring rate), this procedure may thus outcompete other methods that randomly select covariates. But this advantage may be eliminated in more information-rich settings.

RotSF also performs comparably better in the presence of added noisy variables. In this setting, principal component analysis spreads the important variables within each subset into all the principal components in that subset. Thus, compared to other methods, fewer splits based entirely on noisy variables may occur, resulting in more accurate predictions.

While ORSF also transforms the covariates when building trees, it does so by fitting regularized Cox proportional hazards models to the data in each non-terminal node using random subsets of covariates and three different values of one of the regularization parameters. The estimated hazard functions (which are transformations of the subset of covariates chosen at that node) for each choice of regularization parameter become the candidate splitting variables. This model-based approach to building trees is quite different from the usual, entirely non-parametric approach of standard random forests. Perhaps due to these imposed parametric constraints, ORSF did not perform as well as RotSF – including in most settings at which RotSF performed above average. Even in the Y-Weibull setting, where the proportional hazards assumption is met, ORSF did not perform relatively better than other methods. The reason may be that, at each non-terminal node, not all important variables are included in the model, i.e., the model is misspecified.

6 Conclusions

We conducted an extensive review of all recent survival forest methods and a comprehensive simulation study to investigate the performance of some of these methods under various conditions. Based on the results from this study, we summarized relative performance differences in a variety of settings and determined recommendations for the practitioner.

Most methods we retained for the comprehensive simulation study (CIF, RSF-log-rank, RIST-1, RotSF, L ₁, and ORSF) were competitive for point prediction under different settings, so we would not a priori rule out any of these methods for this purpose. For survival function estimation, RSF-log-rank tended to perform consistently below average; other methods were competitive overall, while excelling or struggling in more specific settings.

In contrast, we ruled out RSF-log-rank-score, RSF-C-index, RSF-Brier-score, CRF-GRF, and CRF-QRF, which we tested in the pilot study, due to their more widespread subpar performance. In particular, we found that RSF-log-rank-score was by far the worst-performing method for both point prediction and survival function estimation. Moreover, we removed RIST-3 from the comprehensive study because we found that it did not perform substantially better than RIST-1 in the pilot study but takes considerably more computational time.

Among other notable results, for point prediction, when censoring was low, RotSF performed better than average at our lower sample size and worse than average at our higher sample size. Moreover, when censoring was low, both RIST-1 and RSF-log-rank performed better than average, except at Y-GG (where they performed no differently from average). Interestingly, L ₁ had no advantage over RSF-log-rank for any combination of censoring and sample size levels. CIF and ORSF both performed below average when censoring was low. At our high censoring level, no method performed consistently better than average.

For survival function estimation, perhaps the most notable result was the relative performance advantage of RotSF in information-scarce settings (low sample size and/or high censoring rate) and relative performance disadvantage in information-rich settings (high sample size and low censoring). RIST-1 performed better than average when censoring was low, except at Y-GG. RSF-log-rank performed consistently below average for every combination of censoring and sample size. Whereas ORSF consistently struggled when censoring was high, CIF consistently struggled when censoring was low.

Plenty of opportunities exist for future work. Testing our hypotheses about the reasons for methods’ relative performance advantages and disadvantages is of particular importance. For example, is tree depth a key factor in determining the below-average performance of CIF and above-average performance of RIST in certain settings?

Another factor that could be investigated for its effect on method performance is the signal-to-noise ratio (SNR). We used a constant, relatively high SNR, but relative performance of the methods might differ under different SNR settings. For example, Mentch and Zhou [36] found that, to minimize error, low SNR settings are best accompanied by a smaller number of candidate variables from which to select an optimal split.

We also used the default tuning parameter values of each forest method. Our results might have differed had we tuned each method on each dataset.

Perhaps most importantly for practitioners, an implementation of these methods via packages on CRAN would be tremendously beneficial, especially if one package could amalgamate all methods. (We remind the reader that we obtained code for RIST and L ₁ from the authors who proposed them and the code for RotSF and CRF from their authors’ GitHub pages.) This synthesis would simplify the comparison of the performance of the methods by standardizing the choice of point predictions and error metrics across methods. Without this standardization, inconsistencies across methods can occur. For instance, we used the estimated median as the point prediction when computing the C-index, whereas Ishwaran et al. [6] used “predicted mortality” in their randomForestSRC package. Similarly, the range of the estimated survival function for a given set of covariates differs across packages. Like Kaplan-Meier and other non-parametric methods, most survival forest methods provide probability estimates at each observed survival time in the training set. However, the L ₁ implementation, for example, provides estimates of the survival function only at certain subsets of the observed survival times.

Another challenge worthy of investigation is the issue of point prediction (via the estimated median) when the estimated survival function does not drop below 0.5.

We caution that the datasets that we generated to evaluate method performance should not be misconstrued as a representative sample of all survival datasets. For example, one-third of all FLCs were at Y-GG, but crossing hazards may not be a property seen with such a high frequency in practice. Rather, we sought to investigate a range of possible factors, factor levels, and FLCs that a practitioner might encounter and to understand how methods might interact with them. Moreover, we did not test the methods on real datasets because the number of publicly available survival datasets is very limited, particularly those with a substantial number of covariates. We were unable to locate a repository of such datasets, but it would be a welcome addition to the resources available to researchers in machine learning.

We also recommend investigating method performance at more censoring levels to inform us when some methods – especially those that performed very differently at our 10 % and 70 % censoring levels – start to suffer or excel. Relatedly, do methods that require estimation of the censoring distribution perform relatively poorly with lower censoring rates?

In summary, we have reviewed and compared all recent survival forest methods proposed in the literature and discovered various settings in which some perform relatively better or worse. Our results provide guidance to practitioners about which methods to use for their specific application. Our work also paves the way for future research on the mechanisms underlying these methods and the development of new methods with even better performance.