Random Survival Forests

Consider a cohort of patients discharged after a heart attack. For each patient we record a handful of covariates and the time until readmission. Some patients move to another city, and some are still event-free when the study closes; for those we know only that they were free of readmission up to their last contact. For this right-censored outcome, we record a duration and an indicator showing whether the duration is the event time or a lower bound on it. As we saw in Time-to-Event Data, Censoring, and Hazard Functions, this is right censoring. The model must use those lower bounds without treating them as observed failure times.

The Cox proportional hazards model from Cox Proportional Hazards provides a useful starting point. With a linear predictor and constant coefficients, it models the log hazard ratio as a linear function of the supplied covariates and assumes proportional hazards over time. A threshold effect or an interaction omitted from those covariates can misspecify that predictor. Transformations, splines, and interaction terms can address the problem, but their form must be selected. A survival forest offers a different way to learn nonlinear covariate relationships.

The random forest of Random Forest learns nonlinear relationships for regression and classification by combining trees. Bootstrap samples and random feature subsets vary the trees, while averaging can reduce their prediction variance. Reusing an ordinary regression tree for survival data requires more care: event times hidden by censoring are unavailable for computing the complete-data sample mean and squared-error split criterion. Those population quantities can still exist; the problem is that they are not fully observed.

The random survival forest (RSF) of Ishwaran, Kogalur, Blackstone, and Lauer (2008) adapts the splitting rule and terminal-node prediction to right censoring. We use a log-rank split statistic and a Nelson–Aalen cumulative-hazard estimate in each leaf, retaining bootstrap samples, random feature subsets, and out-of-bag predictions. Our scratch implementation obtains survival by exponentiating the negative ensemble cumulative hazard. Other implementations can use a different survival convention: scikit-survival averages leaf Kaplan–Meier curves for predict_survival_function, while averaging leaf Nelson–Aalen estimates for cumulative-hazard predictions. We will keep those outputs distinct throughout the comparison.

We will build it from scratch in this chapter. First the two components in detail, then a visualization of each, then an eight-patient worked example in which every number can be checked by hand, then a NumPy implementation that we compare with Cox regression and with the reference implementation in scikit-survival.