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.

The first advantage is structural. A survival forest does not require a linear log-hazard predictor, and it does not assume that the hazard ratio between any two covariate profiles is constant over time. A tree partitions the covariate space into rectangles and estimates the hazard separately inside each rectangle, so a covariate can matter in one region and not another, effects can be non-monotone, and interactions are discovered rather than specified. The ensemble is not restricted to a common proportional-hazards form, so two covariate profiles may have survival curves that cross, which a Cox model with constant coefficients cannot produce.

The second advantage is that censored observations contribute through risk sets. A leaf's Nelson–Aalen estimate counts each patient until their observed event or censoring time. It does not fill in an unobserved failure time. Interpretation still requires assumptions: independent censoring conditional on covariates is a standard starting point, but a leaf pools several covariate profiles. Conditional independence given the full covariates need not imply independent censoring within that pooled leaf. Finite leaves, adaptive splitting, and heterogeneity can therefore affect estimation; correct risk-set bookkeeping alone does not guarantee unbiased predictions.

The third advantage is practical. Trees compare covariate values with thresholds, so they do not need standardization to make gradient steps well conditioned. Strictly increasing transformations preserve training-value order and the partitions available to an exhaustive threshold search. Non-strict transformations can merge values, and interpolated quantile thresholds, ties, finite precision, and predictions on new values prevent a blanket claim of identical fitted predictions under every monotone map. Out-of-bag predictions provide a useful internal evaluation without holding out each observation from the entire forest fit. They still require care with preprocessing, tuning, censoring, and the evaluation grid.

One cost is interpretability. A fitted linear Cox model supplies coefficients that describe conditional hazard ratios under its assumptions. A survival forest supplies curves and, often, a single mortality score obtained by summing estimated cumulative hazard over a shared time grid. The score is not a hazard ratio or a probability, and its numerical scale depends on the grid. Permutation importance measures changes in a chosen predictive score; split-based importance measures a different quantity and can favor variables offering many candidate splits. Neither should be interpreted as a treatment effect. Partial dependence and predicted curves can help describe the fitted responses.

Another cost is the amount of information available within each leaf. Log-rank splitting uses event ordering and risk-set counts rather than the spacing between successive event times. Cox partial likelihood also relies on risk sets, so this fact alone does not establish that Cox uses more information or has a smaller mean squared error. A correctly specified proportional-hazards predictor can use its structure effectively, whereas a flexible forest must estimate many local curves. Small event counts and small late risk sets can make those estimates unstable. Event-count and leaf-size constraints help manage that trade-off, without guaranteeing precision.

Survival prediction also requires deliberate conventions. Our forest averages cumulative hazards and then exponentiates the negative average. Averaging the exponentiated tree hazards instead gives a different curve: Jensen's inequality places that average at or above the exponential of the average hazard. This comparison concerns the same tree hazard estimates, not a claim that averaged Kaplan–Meier curves equal either expression. The fitted staircases are flat after their last observed event, which supplies no evidence about later failures. The forest also does not solve informative censoring or automatically correct the consequences of pooling heterogeneous patients.

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

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:

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$:

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

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

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,

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

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,

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

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,

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

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

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.

The predictions are step functions with jumps at the union of the training event times, so they have jumps rather than a smooth time derivative and supply no learned extrapolation beyond the largest observed event time. At that point the cumulative hazard is flat, the survival function is constant, and the model is silent about what happens next.

Before fitting anything, we should see how the two ingredients of the model work. We simulate a dataset in which the truth is known, and use it to look at the split search and at the terminal-node estimate.

The data generator is a Weibull model in which the log hazard is a nonlinear function of the covariates. The Weibull generator uses a shape parameter and a rate parameter. We distinguish the generator parameters from the earlier tree-correlation symbol $\rho$: the Weibull shape parameter is $a = 1.5$ and the Weibull rate parameter is $b = 0.2$, so the conditional survival is $S(t \mid x) = \exp(-b\, t^{a} e^{\eta(x)})$. The covariates are five independent standard normal variables, the first two of which are informative and the last three of which are noise. The covariate contribution to the log hazard is

which contains a linear term, a periodic term, and a quadratic term, and therefore a linear model in $x_1$ and $x_2$ cannot represent the truth. Censoring is exponential with rate $0.30$, which produces roughly half of the observations censored at the generated event-time scale.

The quantities are:

• $x_1, x_2, \ldots, x_5$: the five standard-normal covariates, of which only $x_1$ and $x_2$ enter $\eta$
• $\eta(x)$: the covariate contribution to the log hazard; the Weibull baseline also depends on time
• $a$: the Weibull shape parameter, $1.5$ here
• $b$: the Weibull rate parameter, $0.2$ here
• $S(t \mid x)$: the conditional survival function of the event time given the covariates

The simulated censoring times are independent of both the covariates and the latent event times. This gives independent censoring in this simulation.

The event counts, censoring percentages, and median follow-up above are computed from the simulated arrays and can vary with the seed. The requested row counts remain 500 for training and 400 for testing. In this simulated example, a few hundred events allow several leaves under the minimum-events rule, while censoring makes the risk-set bookkeeping visible.

The first ingredient is the log-rank statistic. The function below implements the formula we derived, with the sorted arrays used to compute the risk sets at each distinct event time. The grouping vector is true for the left child and false for the right child; a negative returned statistic indicates fewer left-child events than expected in the summed comparison, and its absolute value is what the split search maximizes.

The $x_1$ curve has a broad interior trough rather than a monotone shape. The selected threshold near $x_1=0.53$ separates a lower-risk left child from a higher-risk right child in this sample. The $x_2$ curve shows a threshold near $x_2=1.56$ that isolates the upper high-risk tail. A single threshold cannot put both tails of a quadratic effect into one child while leaving the centre in the other. Later splits can isolate the two tails in separate leaves, allowing their predictions to differ from the centre. Bootstrap samples and random feature subsets can change the selected thresholds across trees.

The second ingredient is the Nelson–Aalen estimator. The function below sorts the data, groups the times, and accumulates the increments. It returns the distinct event times and the cumulative hazard at each of them, so the result is a step function that can be evaluated anywhere.

The fitted child staircases separate in this sample; a selected split can also separate curves by chance. The jumps are visible: each event time in a node produces a vertical increment of size $d/Y$, and their sizes depend on the local event count and number still at risk. A minimum event count does not itself bound the last jump. If a leaf had only one event and one patient, the jump would be a full unit of hazard, and the survival function would collapse to $e^{-1}$ in one step. Our implementation constrains both observations and events per child. Other implementations, including the reference used below, express their leaf-size constraints differently.

The formulas are easier to trust once every number has been computed by hand. We use eight patients and a single binary covariate, so that the split is given and the arithmetic is the only thing left to do. Every quantity we compute below maps to one of the equations above: the risk sets and counts feed the split criterion, and the same counts feed the terminal-node estimate.

Group A is the low-risk group: four patients, of whom one has an event at time 8 and three are censored at times 1, 2 and 3. Group B is the high-risk group: four patients with events at times 1, 2, 3 and 4. The table below is the entire dataset.

The printed values reproduce the hand calculation. The left child (group A) is the lower-risk child, so the sign of the log-rank statistic is negative. If any printed value disagrees with the arithmetic above, the implementation has a bug, and the worked example is what makes that comparison possible. The split criterion and the terminal estimate are built from the same three numbers per event time: the risk-set size, the number of deaths, and the group membership. Those counts determine the two survival calculations once a grouping is fixed. The tree also needs the covariates and candidate thresholds to decide that grouping.

We now assemble the pieces into working code. The tree is a recursive function whose only job is to choose the best split, recurse, and store a Nelson–Aalen staircase at the leaves. The forest is a bootstrap loop around the tree, plus the out-of-bag bookkeeping. We use NumPy for the mechanism and then compare with the reference implementation in scikit-survival (sksurv), which supplies the reference estimators and metrics used here.

A node stores a leaf step function, a split feature and threshold, and two child references. The class below is deliberately thin. The __slots__ attribute avoids the overhead of a dictionary per node, which matters when a forest contains hundreds of thousands of nodes.