Random Survival Forests › Formula

We organize the construction around three questions: how to choose a split, what to store at a leaf, and how to combine the trees.

1. How do we choose a split? A classification tree compares class purity in the two children, and a regression tree compares the spread of the numeric outcome. Their usual complete-data formulas require observed labels or numeric outcomes. Censored event times cannot be inserted as though their observed lower bounds were the outcomes. We need a criterion that uses only what is observed: the ordering of events in time and the sizes of the risk sets on each side of the candidate threshold.
2. What do we store at a leaf? A regression leaf stores a numeric prediction; a classification leaf stores class probabilities (a probability vector for multiple classes). A survival leaf has to store a function: the survival curve of the patients in the node, estimated from its own risk sets without assuming a parametric shape.
3. How do we combine the trees? Scikit-learn classification forests average class probabilities, while regression forests average numeric predictions; some classification forests instead combine class votes. A survival forest combines function-valued predictions. Averaging cumulative hazards and averaging survival curves are different choices; either can be meaningful, and averaging alone does not establish unbiasedness or consistency.

Bootstrap sampling, random feature subsets, and recursive axis-aligned partitioning reuse the usual forest scaffold. Event-count stopping rules and out-of-bag scoring also have to account for censored outcomes. We begin with risk sets, derive the split criterion and leaf estimator, then combine the trees.

In the formulas below, "death" denotes the event of interest, including readmission in the opening example.

Everything in the forest is built from the same three quantities observed at each event time: how many patients are still at risk, how many fail, and how the failures are distributed across the children of a split. To write those quantities down, we start with the observation model. Let $T_i$ be the observed time for patient $i$ and $\delta_i$ the event indicator, so that $\delta_i = 1$ for an event and $\delta_i = 0$ for a censored observation. The underlying event time is $\tilde T_i$ and the censoring time is $C_i$, with

$$T_i = \min(\tilde T_i, C_i), \qquad \delta_i = \mathbb{1}\{\tilde T_i \le C_i\}.$$

The quantities are:

• $T_i$: the observed time for patient $i$, equal to the smaller of the event time and the censoring time
• $\tilde T_i$: the underlying (latent) event time for patient $i$
• $C_i$: the censoring time for patient $i$
• $\delta_i$: the event indicator, equal to 1 when the event is observed at or before the censoring time and 0 otherwise
• $\mathbb{1}\{\cdot\}$: the indicator function, equal to 1 when its argument is true and 0 otherwise

The feature vector is $X_i \in \mathbb{R}^p$. The central object in everything that follows is the risk set: the set of patients still observable just before time $t$. Its size is the number at risk:

$$R(t) = \{i : T_i \ge t\}, \qquad Y(t) = |R(t)| = \sum_{i=1}^{n} \mathbb{1}\{T_i \ge t\}.$$

The quantities are:

• $R(t)$: the risk set, the set of patients still under observation just before time $t$
• $Y(t)$: the number at risk at time $t$, the size of the risk set
• $n$: the total number of patients
• $T_i$: the observed time for patient $i$
• $\mathbb{1}\{\cdot\}$: the indicator function, equal to 1 when its argument is true and 0 otherwise

A risk set contains patients whose events could still be observed at time $t$. Under the convention $T_i \ge t$, a patient censored at time 5 contributes at earlier event times and at time 5 itself, then leaves all later risk sets. Dropping that patient at the start would undercount those earlier denominators. Keeping the patient after censoring would overcount later denominators. Risk-set bookkeeping uses the observed follow-up without imputing a later failure time.

The quantity $d(t)$ is the number of deaths observed exactly at $t$:

$$d(t) = \sum_{i=1}^{n} \mathbb{1}\{T_i = t, \delta_i = 1\}.$$

The quantities are:

• $d(t)$: the number of deaths observed exactly at time $t$
• $T_i$: the observed time for patient $i$
• $\delta_i$: the event indicator for patient $i$, equal to 1 for an event and 0 for a censored observation
• $\mathbb{1}\{\cdot\}$: the indicator function, equal to 1 when its argument is true and 0 otherwise

Grouping events at their recorded times also handles tied observations; it does not require an intrinsically discrete event-time model. At any time $t$, we observe the pair $(d(t), Y(t))$: how many failed and how many were eligible to fail. For a fixed assignment of patients to candidate children or terminal leaves, these event-time counts determine the log-rank and Nelson–Aalen calculations. The search still needs covariate values and thresholds to form those assignments.

Now we can state the split criterion. A survival tree recursively partitions covariate space. At each node it compares thresholds on candidate features, forming a left child $L$ and right child $R$. Ordinary squared-error splitting cannot use unobserved failure times as if they were known. Log-rank splitting instead compares the observed event counts with those expected from the two children's risk-set sizes.

What we can compare is what the risk sets tell us at every event time: the fraction of deaths that occur in the left child versus the fraction of patients who were at risk there. Under equal hazards and suitable censoring assumptions, their conditional expectations agree. Observed fractions can differ by chance even for an uninformative split, and crossing hazards can make signed differences cancel. The log-rank statistic sums those observed-minus-expected differences and divides by an aggregate standard deviation. That is the intuition. The formula makes it precise.

For a candidate split, let

$$Y_L(t) = |R(t) \cap L|, \qquad Y_R(t) = |R(t) \cap R| = Y(t) - Y_L(t),$$

be the number at risk in each child, and let $d_L(t)$ be the number of deaths at time $t$ that fall in the left child, where:

• $R(t)$: the risk set at time $t$, the set of patients still under observation just before $t$
• $L$: the set of patients in the left child of the split
• $R$: the set of patients in the right child of the split (this symbol is overloaded with the risk set $R(t)$; the argument $t$ disambiguates the risk set here)
• $Y_L(t)$: the number at risk at time $t$ in the left child
• $Y_R(t)$: the number at risk at time $t$ in the right child
• $Y(t)$: the total number at risk at time $t$
• $d_L(t)$: the number of deaths at time $t$ that fall in the left child

If the two children have the same hazard, then the $d(t)$ deaths at time $t$ are allocated between the children in proportion to how many people were at risk there, so the expected number of deaths in the left child is

$$\mathbb{E}[d_L(t)] = d(t)\,\frac{Y_L(t)}{Y(t)}.$$

The quantities are:

• $\mathbb{E}[d_L(t)]$: the expected number of deaths in the left child at time $t$ under the null of equal hazards
• $d(t)$: the total number of deaths at time $t$
• $Y_L(t)$: the number at risk at time $t$ in the left child
• $Y(t)$: the total number at risk at time $t$

The observed number minus the expected number,

$$d_L(t) - d(t)\,\frac{Y_L(t)}{Y(t)},$$

is the observed-minus-expected event contribution at that time. The log-rank numerator sums these contributions, and its denominator is the square root of their summed null variances. It does not standardize each event-time contribution separately. The variance calculation below accounts for the number of deaths and the allocation of patients between the two risk sets.

The probability of exactly $k$ deaths in $L$ when $d(t)$ deaths are drawn without replacement from a population of $Y(t)$ with $Y_L(t)$ members is hypergeometric, and its variance is

$$V(t) = \frac{d(t)\,\bigl(Y(t) - d(t)\bigr)\,Y_L(t)\,Y_R(t)}{Y(t)^2\,\bigl(Y(t) - 1\bigr)}.$$

The quantities are:

• $V(t)$: the hypergeometric variance of the number of deaths in the left child at time $t$
• $d(t)$: the total number of deaths at time $t$
• $Y(t)$: the total number at risk at time $t$
• $Y_L(t)$: the number at risk at time $t$ in the left child
• $Y_R(t)$: the number at risk at time $t$ in the right child

Write $p_L=Y_L(t)/Y(t)$. For $Y(t)>1$, the same hypergeometric variance is $d(t)p_L(1-p_L)[Y(t)-d(t)]/[Y(t)-1]$. The first factor is the binomial variance for $d(t)$ draws, and the last factor corrects for sampling without replacement. An empty child or an event count equal to the whole risk set gives zero variance. When $Y(t)=1$, the displayed fraction is undefined; the implementation assigns zero because no between-child comparison remains.

Summing the surprises and the variances over all event times gives a single standardized statistic. Writing $t_1 < t_2 < \dots$ for the distinct event times,

$$L(s) = \frac{\sum_{j} \left( d_L(t_j) - d(t_j)\,\dfrac{Y_L(t_j)}{Y(t_j)} \right)}{\sqrt{\sum_{j} V(t_j)}},$$

The quantities are:

• $L(s)$: the standardized log-rank statistic for candidate split $s$
• $s$: the index of the candidate split
• $t_j$: the $j$-th distinct event time
• $d_L(t_j)$: the observed number of deaths in the left child at time $t_j$
• $d(t_j)$: the total number of deaths at time $t_j$
• $Y_L(t_j)$: the number at risk at time $t_j$ in the left child
• $Y(t_j)$: the total number at risk at time $t_j$
• $V(t_j)$: the hypergeometric variance of the left-child death count at time $t_j$

The tree picks the split that maximizes $|L(s)|$. A positive sign indicates more observed left-child events than expected in the summed comparison; it need not mean a larger hazard at every time. For a prespecified split, the standardized statistic has an approximately standard-normal null distribution under the usual log-rank central-limit conditions, including independent censoring, growing nondegenerate event information, and no dominating event-time contribution. With untied failures it is the standardized Cox score at $\beta=0$ for the binary group indicator. Tied-failure conventions can change the Cox information calculation, so that exact equivalence needs its convention stated.

This is a greedy search, and it is biased. The maximum of a large collection of statistics is not itself a standard normal, and features with many distinct values offer more candidate thresholds and therefore more chances to find a large maximum. Because we maximize over thresholds, we should not turn the split value into a $p$-value. The random survival forest mitigates the problem by considering only a random subset of features at each split and by requiring a minimum number of events in each child, but the mitigation is not a cure.

The split criterion gives us a partition, but a partition is only useful if we can describe the survival behavior of each cell. This is the second piece of machinery. Once a node is declared terminal, we need to turn its patients into a survival curve, and we want the curve to assume no parametric shape. The natural starting point is the nonparametric estimator from the counting-process chapter, the Nelson–Aalen estimator. For a terminal node $h$ with event times $t_{h,1} < \dots < t_{h,m_h}$, the number of deaths at $t_{h,j}$ is $d_{h,j}$ and the number at risk just before it is $Y_{h,j}$. The estimated cumulative hazard is

$$\hat{H}_h(t) = \sum_{j\,:\,t_{h,j} \le t} \frac{d_{h,j}}{Y_{h,j}},$$

The quantities are:

• $\hat{H}_h(t)$: the Nelson–Aalen estimate of the cumulative hazard in node $h$ at time $t$
• $h$: the terminal node
• $t_{h,j}$: the $j$-th distinct event time in node $h$
• $d_{h,j}$: the number of deaths at $t_{h,j}$ in node $h$
• $Y_{h,j}$: the number at risk just before $t_{h,j}$ in node $h$

and the implied survival function is $\hat S_h(t) = \exp(-\hat H_h(t))$.

Compare this estimate with the Kaplan–Meier product $\hat S_{\mathrm{KM},h}(t)=\prod_{j:t_{h,j}\le t}(1-d_{h,j}/Y_{h,j})$. Each factor is the empirical conditional survival fraction at an event time. Where every increment $u=d/Y$ is small, the expansion $-\log(1-u)=u+O(u^2)$ explains why the negative log of Kaplan–Meier is close to Nelson–Aalen. They need not be equal in a finite leaf. If the last remaining patient fails, Kaplan–Meier reaches zero while $\exp(-\hat H_h)$ remains positive. Both are valid finite-follow-up survival estimates; neither flat tail identifies the unobserved lifetime distribution beyond follow-up.

The terminal node therefore stores a step function, not a number. Its jumps are at the event times of the node, and its height is determined by the node's own risk sets. Jump sizes depend on $d/Y$ at each event time. A small late risk set can create a large jump even when earlier risk sets were large. This makes event counts and risk-set sizes useful tuning considerations, and it is why the forest's terminal nodes are not "pure" in the sense of classification trees: they are small samples from which a cumulative hazard can be estimated.

A second estimate, the "mortality" or risk score, is used to rank patients and is a summary of the whole step function. It exists because ranking and reporting often demand a single number per patient, even though the primary output is a curve. With a common grid $t_1 < \dots < t_K$ of times shared by all trees,

$$\tau(x) = \sum_{k=1}^{K} \hat H(t_k \mid x),$$

The quantities are:

• $\tau(x)$: the risk score for a patient with covariates $x$, the sum of the ensemble cumulative hazard over the grid
• $x$: the covariate vector of the patient
• $t_1 < \dots < t_K$: the common time grid shared by all trees
• $K$: the number of points on the common grid
• $\hat H(t_k \mid x)$: the ensemble cumulative hazard at time $t_k$ for covariates $x$

which is the sum of the ensemble cumulative hazard over the grid, i.e. the average across trees of the total hazard accumulated over the grid. It is a rank-order score and it is not a probability. A larger $\tau$ means a larger sum on this grid, but does not imply larger cumulative hazard or lower survival at every time when curves cross; two patients whose scores differ by a factor of two are simply in a different order, not "twice as sick".

The third piece is the ensemble. A single tree can have high sampling variance: its splits depend on the training sample, and a small change in the data can move a threshold and change predictions below it. The ensemble can reduce that variance by averaging trees with imperfectly correlated predictions. Bootstrap samples provide diversity and omitted observations for OOB assessment, subject to the evaluation caveats below.

The forest is a collection of trees, each grown on a bootstrap sample of the observations. Let $\hat H_b(t \mid x)$ be the cumulative hazard for the $b$-th tree, evaluated at a point $x$ by pushing $x$ down the tree and returning the step function of the terminal node it reaches. The ensemble estimate is

$$\hat H(t \mid x) = \frac{1}{B} \sum_{b=1}^{B} \hat H_b(t \mid x), \qquad \hat S(t \mid x) = \exp\!\left(-\hat H(t \mid x)\right).$$

The quantities are:

• $\hat H(t \mid x)$: the ensemble cumulative hazard at time $t$ for covariates $x$
• $B$: the number of trees in the forest
• $\hat H_b(t \mid x)$: the cumulative hazard of the $b$-th tree at time $t$ for covariates $x$
• $\hat S(t \mid x)$: the ensemble survival function at time $t$ for covariates $x$

Our scratch forest uses the cumulative-hazard average shown above and defines survival as its negative exponential. For any fixed collection of tree hazards, convexity gives $B^{-1}\sum_b\exp(-\hat H_b)\ge\exp(-B^{-1}\sum_b\hat H_b)$. Both sides are nonincreasing curves in $[0,1]$, but they need not agree. This finite-ensemble inequality supplies no automatic unbiasedness or consistency result. Scikit-survival's cumulative-hazard method also averages Nelson–Aalen estimates; its survival method instead averages leaf Kaplan–Meier estimates. Those survival curves are not the left-hand expression unless the leaf estimators happen to coincide. Check the actual prediction method before comparing implementations.

The following deterministic toy ensemble isolates that gap. Five nondecreasing cumulative-hazard curves share one baseline shape but differ in scale, so the difference between the two aggregation conventions is not confounded by sampling noise.

The bootstrap also supplies the out-of-bag estimate. Observation $i$ is out of bag for tree $b$ if the $b$-th bootstrap sample did not contain it. Let $B_i = \{b : i \notin \text{bootstrap}(b)\}$ be the set of trees for which $i$ is out of bag. The probability that a given observation is in a given bootstrap sample is $1 - (1 - 1/n)^n \to 1 - e^{-1} \approx 0.632$, so its exact expected size is $B(1-1/n)^n$, approximately $0.368B$ for large $n$. The out-of-bag prediction is

$$\hat H_{\mathrm{OOB}}(t \mid x_i) = \frac{1}{|B_i|} \sum_{b \in B_i} \hat H_b(t \mid x_i),$$

The quantities are:

• $\hat H_{\mathrm{OOB}}(t \mid x_i)$: the out-of-bag cumulative hazard at time $t$ for observation $i$
• $x_i$: the covariate vector of observation $i$
• $B_i$: the set of trees for which observation $i$ is out of bag
• $|B_i|$: the number of such trees
• $\hat H_b(t \mid x_i)$: the cumulative hazard of tree $b$ at time $t$ for observation $i$

Here only trees whose bootstrap samples exclude observation $i$ contribute to its prediction. If $|B_i|=0$, that prediction is unavailable. This avoids training those trees' splits and leaves on the held-out observation, but it does not make the complete OOB procedure independent of that observation: this implementation uses the full training event-time grid when forming mortality scores. Shared preprocessing or repeated OOB tuning can introduce further dependencies. Treat OOB concordance as an internal estimate, and reserve an independent test set for a final comparison.

Three properties summarize what we have built, and each points to a different reason to be careful when the forest is applied.

For a prespecified binary grouping with untied failures, the log-rank comparison agrees with the standardized Cox score at zero coefficient. This connects the split criterion to a local proportional-hazards comparison without imposing proportional hazards on the fitted forest. Choosing the largest statistic over many thresholds changes its null distribution, and crossing hazards can reduce the power of the signed comparison.

The Nelson–Aalen estimator adds the non-negative increments $d_{h,j}/Y_{h,j}$. At a fixed time $t$, a commonly used variance estimate is $\sum_{j:t_{h,j}\le t}d_{h,j}/Y_{h,j}^2$. This estimates uncertainty in the cumulative-hazard staircase; it is distinct from the hypergeometric variance used in the log-rank split statistic.

Averaging can reduce variance; it does not in itself remove bias. If random tree predictions at a fixed $(t,x)$ have common variance $\sigma^2$ and common pairwise correlation $\rho$, the covariance sum gives $\operatorname{Var}(B^{-1}\sum_b\hat H_b)=\rho\sigma^2+(1-\rho)\sigma^2/B$. This is an exchangeable-prediction model, not a claim that trees are unbiased. Independently randomized trees can be conditionally independent given the training data while remaining dependent across repeated training datasets. More trees reduce the remaining randomization variation, with diminishing gains; feature subsampling can change both correlation and individual-tree strength.

The variance formula can be read directly from a small curve. For a fixed unit individual-prediction variance, the floor is the tree correlation, and more trees cannot remove variance that the trees share.