Survival outcome has been one of major endpoints for clinical trials; it gives information on the probability of time-to-event of interest. Recently, there has been increasing demands for survival analysis tools, especially for high dimensional survival data. Most of the classical statistical approaches including nonparametric, semi-parametric and complete parametric models, have limitations. Typical nonparametric approaches, such as the logrank (or Cox-Mantel) test, do not have stringent models assumptions, but can only deal with a limited number of categorical predictors. Typical semi-parametric approaches, such as Cox proportional hazard (PH) model, have stringent model assumptions, such as linearity, interactions and proportionality; and they can only deal a limited number of predictors. Complete parametric models, such as accelerated failure time models, are similar to semi-parametric models except that they make further assumptions about the baseline hazard function.

In this paper, non-parametric random survival forest approaches were introduced with many advantages over typical nonparametric, semi-parametric or parametric approaches. They have very no model assumptions and can deal with many more predictors than typical survival models. Two Cox PH models were evaluated and the model performances were compared with the random survival forest models based on a simulation study. Additionally, the random survival forest approaches were further studied based on a real world case study, which included more predictors then a typical Cox PH model can handle.

Keywords: cox proportional hazard model, nonlinear, proportionality, random forest, survival, prediction errors, AUCs, time-dependent, time-varying

Survival data arise from many fields, such as medicine, biology, public health, epidemiology, etc. Although the analysis approaches discussed in this paper should be general and applicable to all of these disciplines, the research interest focused on biological, medical or clinical applications.

Most if not all survival analysis approaches are focused on survival (or time-to-event) outcomes, which are usually associated with serious disease conditions, such as death, heart failures and recurrence of cancers. Typically, survival outcome includes one binary variable for occurrence(s) of the event(s) of interest and at least one continuous variable for the time of the occurrence(s) or the censoring time (2 time variables may be needed for the intervals within which the event occurs). The outcome is referred to as time-to-event, which has been one of the major endpoints for clinical studies.

To best utilize such data, different approaches are developed to exploit the complete information about the occurrence of the event and the timing of the event. Recently, there has been increasing demand for survival analysis tools, especially for high dimensional survival data. Most of the classical statistical approaches including nonparametric, semi-parametric and complete parametric models have limitations. Typical nonparametric approaches, such as the logrank (or Cox-Mantel) test, do not have stringent model assumptions, but can only deal with a limited number of categorical predictors. Typical semi-parametric approaches, such as Cox PH models, however rely on more stringent model assumptions, such as linearity, interactions and proportionality; and they work only if the number of predictors is fewer than the number of events. Complete parametric models, such as accelerated failure time models, are similar to semi-parametric models except that they make further assumptions about the baseline hazard function.

In this paper, we employed a systematic process to assess the assumptions for typical Cox models. For linearity assumption, appropriate nonlinear functional forms were recommended if necessary; for interaction effects, appropriate interactions were suggested if detected; for proportionality assumptions, appropriate extension model was attempted to resolve non-proportionality.

In summary, two typical Cox PH models were intended, Cox PH linear model without consideration of nonlinearity and Cox PH model with RCS transformations for adjustment of nonlinearity; two nonparametric random survival forest approaches, a logrank based random survival forest¹(LR-RSF) developed by Mogensen and Ishwaran (2012) and a conditional-inference based random survival forest^2,3(CINF-RSF) suggested by Hothorn et al. (2006) were also implemented. All of these intended approaches were compared in terms of prognostic factor detections and prediction performance.

These approaches, including the Cox PH linear model and Cox PH model with RCS transformation, LR-RSF and CINF-RSF models, were assessed via a simulation study; the performance of these models was compared in terms of prediction errors and prediction AUCs. The LR-RSF and CINF-RSF models were further evaluated based on a real world case study, where typical Cox PH models did not work due to too many predictors and too little events.

For both studies, the original data were randomly partitioned into training and testing sets with 3:1 ratio; the training set was used for model selections and cross validations; the test set was used for assessment of the prediction performance

Simulation studies with time-varying treatment effect

A total of 2000 subjects were simulated using R⁴ Software; 7 factors included age, sex, race, systolic blood pressure (SBP), diastolic blood pressure (DBP), body mass index (BMI) and treatment were generated from computer programs. SBP and DBP were generated from two normal distribution with correlated coefficients of 0.9; and they were weakly correlated with BMI (SBP and BMI: $\rho$ =0.3; SBP and BMI: $\rho$ = 0.2). Once all factors were generated, time-to-events were simulated from the factors following a pre-specified exponential model. A time-varying treatment was also generated: subjects were randomized to active and placebo treatment at baseline. However, subject randomized to placebo at baseline were allowed to switch to active treatment after the 1 year follow-up. The actual reasons for treatment switch could be very subjective, but for this study, the treatment switch was randomly generated for the placebo randomized subjects. Further simulation details are not presented in this paper.

Real world case study

A breast cancer data with 5 clinical factors and 70 gene signatures obtained from the Netherlands Cancer Institute was used to evaluate metastasis-free survival, which will be called NKI70⁵ data onward. The data included a total of 144 independent lymph-node-positive breast cancer subjects, followed for 17 months; 48 subjects experienced a metastasis event.

For this study, the number of factors were more than the number of event $\left(p>>N\right)$ ; typical Cox models did not work, only LR-RSF and CINF-RSF models were implemented.

Analysis method

In general, the following 5 functions are used to describe survival data: (1) survival function, (2) density probability function, (3) cumulative probability function, (4) hazard function and (5) cumulative hazard function. A basic problem in survival analysis is to estimate one or more of these five functions from the collected survival data and to draw inferences about the survival patterns in the population.⁴ The formulations of the 5 functions do not involve any assumptions; hence they are still considered as nonparametric and can be applied across all survival models.

Cox PH linear model: Typical Cox PH linear model are commonly expressed as,

$\lambda \left(t|X\right) =\lambda \left(t\right)\exp \left(X\beta \right)$  (1)

where all factors or composite factors are included in their 1st order linear form and potential 2nd order interaction terms between the original factors (or composite factors) are also considered; this model is referred to as Cox PH linear model. For this model, if composite factors are derived, they should be included to replace the original factors they are clustered from. The best model are achieved following a backward step-down selection procedure. This model was implemented with R/rms⁷ package from Harrell (2014).

Cox PH model with RCS transformation: A Cox PH model with RCs transformation is implemented to account for nonlinearity; potential interactions between the actual functional forms of the factors are also considered. For this model, nonlinear continuous factors are transformed using restricted cubic spline (RCS) function. If interactions are detected, the terms should include the RCS transformations of the factors involved in the interactions. For time varying treatment effect, the subject who switched treatment were considered as having two different observations, each with a different treatment and a corresponding time interval; Andersen-Gill extension of the Cox PH model are implemented for the adjustment. The Cox PH models with RCS transformations can be formulated similar to Equation 1, except that some covariates may have nonlinear forms of the original factors. This model were implemented with R/rms⁷ package.

Random survival forest (RSF) models: Two forms of random survival forest (RSF) were implemented. One was based on the logrank test; this approach is referred to as logrank-based RSF model (LR-RSF). Another form was built on conditional inference based logrank test; where multiple significance tests are computed at each start of the algorithm based on permutation test; this approach is referred to as conditional-inference-based (CINF-RSF) model. Both of the RSF mode are based on logrank scores; of which, the ${\chi }^2$ test statistics can be obtained as

${\chi }^2={\frac{\left[{\displaystyle {\sum }_{i=1}^k\left(d_i^{\left(a\right)}-r_i^{\left(a\right)} d_i/r_i\right)}\right]}{{\displaystyle {\sum }_i^k\frac{r_i^{\left(a\right)}{r}_i^{\left(b\right)}d\left(r_i-d_i\right)}{r_i^2\left(r_i-1\right)}}}}^2$  (2)

where the superscripts $\left(a\right)$ and $\left(b\right)$ refer to the two offsprings of the parent node from the tree splitting; $k$ is the number of distinct event times and $i$ refers to the time $i$ .

The LR-RSF model was implemented with R/random Survival SRC^8-10 package developed by Ishwaran and Kogalur (2013) and CINF-RSF model was implemented with R/party^2,11,12 package developed by Hothorn et al. (2006).

Evaluation of model performance: The model performance is evaluated using Brier scores¹³ as proposed by Gerds et al. (2006) and time-dependent AUCs. The Brier Score at time, $t$ , is defined as:

Brier Score at $t$ $L(t)=1/N{\displaystyle {\sum }_{j=1}^N\left(Y\left(t\right)\right)}-{\left(\stackrel{⌢}{Y}\left(t\right)\right)}^2$ (3)

Where $Y\left(t\right)$ is the observed survival status at time $t$ , which is $I\left(T_j\ge t\right), \stackrel{⌢}{Y}\left(t\right)$ is the predicted survival probability at time $t$ , $N$ is the total number of subjects. The prediction error, $t$ , is used to estimate the prediction performance at time $t$ . The integrated Brier score is obtained by averaging the total Brier scores for interval $(0,t]$ , which can be expressed as $1/k{\displaystyle {\sum }_{i}Br\left(t\right)}$ , where $i=\mathrm{1....}k$ different event times.

AUC is another assessment for model performance, which is estimated using c-index, as recommended by Harrell (2012); c-index is the probability of concordant pairs of predicted and observed cases among all pairs of responses. Specifically, if the predicted survival time (or probability) is larger for the subject who (actually) lived longer, the predictions for that pair are said to be concordant with the (actual) outcomes. Similar to prediction errors, time-dependent AUCs can be estimated.

The performance statistics and corresponding 95 percentile intervals (95% PCIs) can be estimated with the mean and the 95 percentiles from bootstrapped samples of the test set.

Results of the simulation study

Summary statistics of the simulation study: Summary of demographics for the simulated study is not presented. Summary of subject survival status at 1-year post baseline and treatment switching after 1-year post baseline is presented in Table 1. For this simulation study, 106 (10.6%) subjects from active treatment arm had events prior to 1 year and 163 (16.3%) subjects from the placebo (SOC) arm had events prior to 1 year. Of the 891 subjects randomized to active treatment who survived 1 year, none of the subjects switched to placebo; of the 840 subjects randomized to placebo and survived 1 year, 201 (23.9%) switched to the active treatment, and the rest of the 639 subjects stayed in the placebo arm until censored or failed.

Preparations for analysis: As mentioned earlier, subjects who were randomized to placebo at baseline were allowed to switched to active treatment after 1-year post baseline; according to the intent-to-treat (ITT) principal, these subjects were assumed to be treated with placebo as intended and assumed to follow the hazard from the placebo arm even after the switching; neglecting the benefit from the active treatment due to treatment switching led to violation of proportionality, the scaled Schoenfeld residual was significantly correlated with time $\left(p=-0.182and << 0.0001\right)$ ; the downwarding slope of the scaled Schoenfeld residual plot as shown in Figure 1 also confirmed non-proportionality.

Figure 1 Scaled Schoenfeld Residuals for Treatment (ITT Principal): Subjects were assumed in their Original Treatment Group as Randomized – Simulation Study.

Apparently, treatment was changing over time for these subjects; to adjust for time-varying treatment effect, each of each subjects were assumed to have multiple observations corresponding to the different treatments overt time, Andersen-Gill extension of Cox PH model was implemented. The scaled Schoenfeld residuals from the adjusted model are plotted in Figure 2. The scaled Schoenfeld residual plot is almost parallel to the x-axis; the proportionality test was no longer significant ( $p-value\approx 0.41$ ).

Figure 2 Scaled Schoenfeld Residuals for Treatment vs. Time with adjustment of Time-Varying Treatment Effect – Simulation Study.

Thus, non-proportionality was resolved for treatment.

Additionally, SBP and DBP were highly correlated, from which a composite factor MAP was derived; thus the factor MAP would be used to replace SBP and DBP for survival models.

Linearity, interactions and proportionality: For Cox PH linear model, linearity was assumed for all factors; for Cox model with RCS transformations, the functional forms were visually examined with the Martingale residual plot, which was obtained from a tentative model with the original factors. Figure 3 displays the Martingale residuals vs. each of the continuous factors stratified by the actual Treatment. The observational-level residuals are presented in the top half of the figure and the subject-level residuals are presented in the bottom. From the subject-level residual plot, 3 inflection points were observed for Age and MAP, 1 inflection point was observed for BMI; to be more conservative, 5-knot RCS transformation for Age and MAP, 3-knot RCS transformation for BMI was considered in the initial model with RCS transformations.

Figure 3 Martingale Residuals against Continuous Factors Stratified by Treatment – Simulation Study.

Additionally, the smoothed lowess curves of Age were parallel between treatments, which suggested no interaction between Age and treatment. For MAP, the lowess curves for different treatments had different slopes, thus interaction was suggested. For BMI, the placebo group had shown slight quadratic pattern, but pretty much straight for the Active treatment; however the overall pattern did not suggest interaction. For further investigation, a model with all interactions involving treatment was intended. Results are presented in Table 2; significant P-value was a strong evidence of interaction. Overall, Age, BMI and Sex had no interactions with Treatment, but MAP and Race might interact with Treatment.

Similarly, interactions involving Sex and interactions involving Race were also examined separately. Interaction between Sex and Race was suggested; interactions between Sex and MAP, between Race and BMI were questionable. For the initial survival model, all potential interactions were included; the initial Cox PH model with RCS transformations is formulated as

$\begin{array}{l}\lambda \left(t|X\right) {\beta }_1^{T}RCS (Age,5) +{\beta }_2^{T}RCS (MAP,5) +{\beta }_3^{T}RCS (BMI,3) +\\ {\beta }_4^{T}I (Treatment) +{\beta }_5^T\left(Sex\right)+{\beta }_6^{T}I\left(Rice\right)+\\ {\beta }_9^{T}RCS (MAP,5): I(Treatment)+\\ {\beta }_9^{T}RCS (MAP,5):I\left(Sex\right) +{\beta }_3^{T}RCS (BMI,3): I(Race)+\\ {\beta }_{10}^{T}I (Treatment) :I(Race) +{\beta }_{11}^T\left(Sex\right): I(Race)\end{array}$

Proportionality assumptions were further examined using a Cox PH model with RCS transformation as presented in Equation 4; the Scaled Schoenfeld residual plots are presented in Figure 4 and the results of the Wald tests for proportionality are presented in Table 3. Non-proportionality was not found for any of the predictors.

Figure 4 Scaled Schoenfeld Residual Plot for Each Covariate Based on the Cox PH Model with RCS Transformation – Simulation Study.

Analysis (Model selection)

Cox PH linear regression: The Cox PH linear model was implemented including all factors in their original scale and the 2nd order interactions. Models were selected using backward step-down procedure with AIC as the selection rule. Once the best model was achieved, the regression coefficients, the corresponding standard errors and P-values were estimated from the final Cox PH linear model; results are presented in Table 4.

The log hazard and hazard ratios from the final model are summarized in Table 5. For continuous factor, the low and high are the minimum and maximum of the factor; the difference was the difference between the maximum and the minimum; the effect, the lower and upper 95% CI were the log hazard (/ hazard ratios), the lower and upper bounds of the 95% CIs with respect to unit change in the factor.

In this table, the interaction effects were not presented separately from the individual factors, instead they were summarized by different level of the factors involved in the interactions: the hazard ratio between Active Treatment and Placebo was obtained based on the median MAP.

The performance for the final Cox PH linear model were measured with prediction errors and time-dependent AUCs, which were assessed based on the test set; results are summarized in Table 6. The prediction errors were calculated based on Equation 3; the AUC was the area under the ROC curve; a value of 0.5 corresponded to the probability of a random guess and a value of 1 corresponded to 100% accuracy; smaller prediction error or large AUCS suggested better predictions. The prediction errors and time-dependent AUCs and the corresponding 95% PCIs are displayed in Figure 5. The survival probability of future events and median survival time for each subject can be predicted; results are not presented.

Figure 5 Prediction Performance of the Selected Cox PH linear model – Simulation Study

Cox PH model with RCS transformation: The initial Cox PH model with RCS transformation was intended with Equation 4. Similar to what was done for Cox PH linear model, the above model was backward selected following a step-down procedure using AIC as the selection rule. The final Cox PH model with RCS transformation obtained from the procedure is summarized in Table 7. The final reduced model can be formulated as

$$\mathrm{Pr}ob\left\{T\ge t\right\}=S_0{\left(t\right)}^{E^{X\beta }},WhereX\stackrel{⌢}{\beta }=$$

$\begin{array}{l}-5.5101+6.24306\left\{Male\right\}+0.3262\left\{Black\right\}+0.0454\left\{Hispanic\right\}+0.0282\left\{White\right\}-0.4751\left\{Treatment\right\}\\ +0.0260BMI-0.0004{\left(BMI-23.7275\right)}^3+0.0009{\left(BMI-28.6586\right)}_{{}^{+}}^3-0.0004{\left(BMI-33.6296\right)}^3+\\ +0.0625MAP-0.0004{\left(MAP-76.7173\right)}_{{}^{+}}^3+0.0025{\left(MAP-83.3120\right)}_{{}^{+}}^3-0.0046{\left(MAP-86.7914\right)}_{{}^{+}}^3\\ +0.002977986{\left(MAP-90.1678\right)}_{{}^{+}}^3-0.0005{\left(MAP-96.5427\right)}_{{}^{+}}^3+\left\{Male\right\}[-0.0837MAP+0.00138{\left(MAP-76.7173\right)}_{{}^{+}}^3\\ -0.0106{\left(MAP-83.3120\right)}_{{}^{+}}^3+0.0202{\left(MAP-86.7914\right)}_{{}^{+}}^3-0.0134{\left(MAP-90.1678\right)}_{{}^{+}}^3+0.0023{\left(MAP-96.5421\right)}^3+]\end{array}$

The log hazards and hazard ratios from the selected Cox model with RCS transformations are presented in Table 8. For continuous factors, the effect, the corresponding lower and upper bound of the 95% CI with respect to unit change are presented.

The prediction performance of the final Cox model with RCS transformations was measured with prediction errors and time-dependent AUCS again based on the test set; the prediction errors, the time-dependent AUC and the corresponding 95% PCIs are summarized in Table 9; the corresponding plots are displayed in Figure 6. With the selected model, the predicted log-hazards, survival and median survival time were obtained for each subject in the test; the results are not presented in this paper.

Figure 6 Prediction Errors and Time-Dependent AUCs of the Reduced Cox PH Model with RCS Transformations – Simulation Study.

Logrank based random survival forest (LR-RSF): As mentioned earlier, LR-RSF could not handle multiple observations per subject from the time- varying treatment effect; thus subjects who switched from placebo to active treatment were considered as 2 independent subjects, each with a different treatment for different durations. Nonetheless, LR-RSF is nonparametric without any assumptions, thus nonlinearity or non-proportionality was not concerned.

As previously mentioned, MAP was a composite factor derived from SBP and DBP, they were highly correlated. To keep consistent with the Cox PH models, the initial RSF models were intended with Age, MAP, BMI, Sex, Race and Treatment. Table 10 presents the VIMP of all factors excluding SBP and DBP. The plots of out-of-Bag (OOB) error rate and VIMP with SBP and DBP excluded are presented in Figure 7; significant change was not observed.

Figure 7 Out-of-Bag Error Rate and VIMP of LR-RSF with SBP and DBP Removed – Simulation Study.

The cross validated LR- RSF model was then used for predictions of future events and assessment of prediction performance. The OOB survival, cumulative hazard and hazard for a subset of 3 subjects are presented in Figure 8 and the overall OOB survival probability, Brier scores and mortality are presented in Figure 9. The predicted mortality and survival probability for each subject in the test set were obtained; results are not presented in this paper. The model performance were assessed again with the test set, the results will be presented later in Table 11 and Error! Reference source not found. for comparison with CINF-RSF model.

Figure 8 CV Survival, Cumulative Hazard and Hazard for LR-RSF (Subset of 3 Subjects) – Simulation Study.

Figure 9 OOB Survival, Brier Scores and Mortality for LR-RSF (All Subjects) – Simulation Study.

Conditional inference based random survival forest (CINF-RSF): Another nonparametric CINF-RSF model was also cross validated with the training set. A CINF-RSF tree is presented in Figure 10; survival plots are displayed in the terminal node.

Figure 10 A Sample Forest Tree from CINF-RSF – Simulation Study.

LR-RSF vs. CINF-RSF: For detecting prognostic factors, both managed to pick up the most important factors, although slight difference was observed between the two RSF models. The prediction performance of the two approaches were further assessed and compared in terms of prediction errors and time-dependent AUCs.

The prediction errors and time-dependent AUCS for LR-RSF and CINF-RSF and the corresponding 95% PCIs are summarized in Table 11 and the prediction errors and AUCs from the reduced Cox PH linear model are also presented as a reference; the corresponding plots are displayed in Figure 12. The black solid curve and the gray shaded area is the performance and the corresponding 95% PCIs for CINF- RSF; the red dashed curve with light pink shaded area is the performance and the 95% PCIs for LR-RSF; the blue dotted curves is the performance from the Cox PH linear model.

Figure 11 Prediced Survival from CINF-RSF and Kaplan Meier Curve – Simulation Study.

Figure 12 Prediction Errors and Time-Dependent AUCs for LR-RSF, CINF-RSF and Conventional Cox Linear Model– Simulation Study.

In terms of prediction errors, the Cox PH linear model performed slightly better than the two RSF models; the two RSF models were similar to each other within the first 5 years; LR-RSF was worse than the CINF-RSF model between 5 and 10 years and at the tail (beyond year-10), LR-RSF performed slightly better than CINF-RSF. In terms of time-dependent AUCs, the Cox PH linear model still performed the best; the two RSF models had the maximum difference within the first 5 years, after which the two models were similar.

Comparisons of all intended models: In summary, the Cox PH linear model had the best performance of all intended models; while the 2 RSF models were much easier to implement and the prediction performance was satisfactory. Both RSF models can deal with more predictors than typical Cox PH models can handle; they are nonparametric without stringent assumptions, so they were much easier to implement. Between the 2 RSF models, LR-RSF is more favorable if covariates are not highly correlated; while CINF-RSF performs better for highly correlated data. There is one pitfall for LR-RSF model, it could not account for multiple observations per subject, which could lead to the slightly worse performance.

For cross comparisons, the prediction errors and time-dependent AUCs for the 4 intended models are displayed again in Figure 13 and Figure 14, respectively. Of the 4 models, the Cox PH linear model and CINF-RSF model had almost equivalent prediction errors and LR-RSF model had the worst prediction errors in general. Surprisingly, the Cox model with RCS transformation had slightly worse prediction errors than the Cox PH linear model, probably because the former picked up BMI instead of Age. In terms of time-dependent AUCs, the Cox PH linear model performed the best overall; the AUC of the LR-RSF inclined a little between year 2-5, otherwise, it had similar performance to the CINF-RSF; the AUC of the Cox PH model with RCS transformation was similar to the CINF-RSF within the first half of the study, then it became the second best at the tail (slightly worse than the Cox PH linear model). Additionally, both Cox PH models provided unbiased parameter estimates which were easy to interpret; while the RSF models did not provide parameter estimates and were relative more difficult to interpret.

Figure 13 Model Comparisons Based on Prediction Errors – Simulation Study.

Figure 14 Model Comparisons Based on Time-Dependent AUC(t) – Simulation Study.

However, considering the nonparametric nature, the prediction performance of the two RSF models was satisfactory and both models picked up the right set of factors without making stringent assumptions.

Results of real world case study

Summary of the NKI70 data: For NKI70 data, only the 5 clinical factors and 2 of the 70 gene signatures are summarized in Table 12. Contingency tables for categorical factors are presented in the top half of the table and continuous factors, including Age and 2 of the 70 gene signatures are summarized in the bottom. Full summary of all gene signatures will not be presented but can be provided upon request. One clinical factor, grade of tumor, had 3 categories; 2 dummy variables were created for the relative difference between categories. Thus, a total of 76 factors should be considered.

Data preparation: The Kaplan-Meier (KM) estimates of the survival probability and the Fleming-Harrington estimates of the log cumulative hazard are displayed in Figure 15; the green dotted line superimposed on top of the cumulative hazard was the hazard assuming exponential survival models; where little change was observed in the slope of cumulative hazard over time. Thus, it was reasonable to consider Cox model. Figure 16 displays the correlation map of all factors; blue color indicated a positive correlation and red color indicated a negative correlation. The intensity of the color indicated the severity of collinearity. Apparently, several of the gene-signatures were highly correlated.

Figure 15 Kaplan Meier Survival Probability and Fleming Harrington Cumulative Hazard for NKI70 Data.

For this case study, typical Cox models did not work for too many predictors with too few events; only the LR-RSF and CINF-RSF models were intended. Again, the two RSF models were nonparametric without stringent assumptions.

Logrank Based Random Survival Forest (RSF): For this case study, LR-RSF was first implemented. Table 13 presents the variable importance (VIMP) of a subset of factors; the gene signature, GPR180, was the least important with VIMP of -0.0027; 3 genes with the least absolute VIMP (VIMP = 0) were Grade. Intermediate, Contig40831.RC and PECI.1; the factors with VIMP ≥ 0.0027, which is the absolute value of the least VIMP, should be potentially important. Figure 17 displays the CV error and VIMP for all factors; the LR-RSF model was stabilized at 25 trees, although the model still exhibited some variations which was common for bootstrap aggregations. For VIMP plot, factors were ordered from the least to the most VIMP; the factors with negative VIMP were colored in red and positive VIMP were colored in blue.

The OOB survival, Brier scores and mortality for all subjects in the training set were obtained from cross validation (CV); results are plotted in Figure 18. The cross validated LR-RSF model was used for predictions based on the test set; the predicted mortality and survival were obtained, but the results are not presented in this paper. The prediction errors and time-dependent AUCs were assessed again based on the test set; results will be presented later in Table 14 for comparison with CINF-RSF model.

Figure 16 Factor Correlation Map for NKI70 Data

Figure 17 CV OOB Error Rate and VIMP from LR-RSF – NKI70 Data

Figure 18 Predicted Survival, OOB Errors and Mortality for LR-RSF– NKI70 Data

Conditional inference based random survival forest (RSF): CINF-RSF model was evaluated with the NKI70 data; a CINF-RSF tree is presented in Figure 19; negative survival probabilities are observed in some terminal leaves, which is possibly due to the over-optimistic predictions. Besides the forest tree, the cross validation also selected ZINF533, EGLN1, UCHL5, PRC1, MTDH, C20orf46, LGP2 and RTN4R1 as important factors.

Figure 19 Forest Tree from CINF-RSF – NKI70 Data.

The predicted survival estimated from the cross validated CINF-RSF model based on the test set and the corresponding KM estimates are presented in Figure 20. The difference in the curves suggests that the CINF-RSF model were probably over-optimistic. The actual reasons for this optimism may be partly due to the extra loss of information during the split of the forest tree.

Figure 20 Prediced Survival Probability from CINF-RSF vs. KM Curve – NKI70 Data.

The performance of the 2 RSF models was assessed again based on the test set. The prediction errors, time-dependent AUCs and the corresponding 95% PCIs of the 2 RSF models were obtained based on 200 bootstrap samples from the test set; results are presented in Table 14; the corresponding plots are displayed in Figure 21, respectively. In terms of prediction errors, the two RSF models were almost equivalent most of the time; but at the tail, the CINF-RSF model was slightly better. In terms of the time-dependent AUCs, the CINF-RSF model was consistently better than the LR-RSF model, and the difference was quite substantial.

Figure 21 Prediction Errors and AUCs for LR-RSF and CINF-RSF – NKI70 Data.

The project was completed for research and publication only and did not receive any financial support.

Special thanks go to Netherlands Cancer Institute for publishing their 70-gene-signature breast cancer data, which was utilized and analyzed in this project. I would also like to express my thanking to Professor Dinesh P. Mital for his valuable guidance and advice.

None.

1. Mogensen UB, Ishwaran H, Gerds TA. Evaluating Random Forests for Survival Analysis Using Prediction Error Curves. J Stat Softw. 2012;50(11):1–23.
2. Hothorn T, Bühlmann P, Dudoit S, et al. Survival Ensembles. Biostatistics. 2006a;7(3):355–373.
3. Hothorn T, Hornik K, Zeileis A. Unbiased Recursive Partitioning: A Conditional Inference Framework. Journal of Computational and Graphical Statistics. 2006b;15(3):651–674.
4. R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria; 2014.
5. van de Vijver MJ, He YD, van't Veer LJ, et al. A gene-expression signature as a predictor of survival in breast cancer. N Engl J Med. 2002;347(25):1999–2009.
6. Balakrishnan N, Rao CR. Handbook of Statistics. Advances in Survival Analysis. 1st ed., HandbNorth Holland; Handbook of Statistics (Book 23). 2004.
7. Harrell Jr FE. rms: Regression Modeling Strategies. R package version 4.2-0; 2014.
8. Ishwaran H, Kogalur UB. Random Forests for Survival, Regression and Classification (RF-SRC), R package version 1.4. 2013.
9. Ishwaran H, Kogalur UB. Random survival forests for R. R News. 2007;7(2):25–31.
10. Ishwaran H, Kogalur UB, Blackstone EH, et al. Random survival forests. Ann Appl Statist. 2008;2(3):841–860.
11. Strobl C, Boulesteix AL, Zeileis A, et al. Bias in Random Forest Variable Importance Measures: Illustrations, Sources and a Solution. BMC Bioinformatics. 2007;8(25).
12. Strobl C, Boulesteix AL, Kneib T, et al. Conditional Variable Importance for Random Forests. BMC Bioinformatics. 2008;9(307).
13. Gerds TA, Schumacher M. Consistent estimation of the expected Brier score in general survival models with right-censored event times. Biom J. 2006;48(6):1029–1040.