Random Survival Forests › Implementation in scikit-learn and scikit-survival
Machine Learning from ScratchAll chapters Random Survival Forests › Implementation in scikit-learn and scikit-survivalSection view · Sections in this chapter7 / 10 ↓
Reading tools Definitions, links & notes
Choose your expertise level to adjust how many terms are explained. Beginners see more tooltips, experts see fewer to maintain reading flow. Hover over underlined terms for instant definitions.
Article links
Make inline references clickable
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.
From a log-rank split to a survival tree
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.
PythonIn [22]Expand
class SurvivalNode:
"""A node in a survival tree; leaves hold a Nelson-Aalen step function."""
__slots__ = ("times", "cumhaz", "feature", "threshold", "left", "right")
def __init__(self):
self.times = np.empty(0)
self.cumhaz = np.empty(0)
self.feature = None
self.threshold = None
self.left = None
self.right = None
@property
def is_leaf(self):
return self.left is None
class SurvivalTree:
"""A single survival tree grown by maximising the absolute log-rank statistic."""
def __init__(
self,
max_depth=7,
min_leaf=15,
min_events=8,
n_cuts=10,
mtry=None,
random_state=0,
):
self.max_depth = max_depth
self.min_leaf = min_leaf
self.min_events = min_events
self.n_cuts = n_cuts
self.mtry = mtry
self.random_state = random_state
def fit(self, X, time, event):
X = np.asarray(X, dtype=float)
time = np.asarray(time, dtype=float)
event = np.asarray(event, dtype=bool)
self.n_features_ = X.shape[1]
self.mtry_ = self.mtry or max(1, int(np.sqrt(self.n_features_)))
self.rng_ = np.random.default_rng(self.random_state)
self.root_ = self._grow(X, time, event, np.arange(X.shape[0]), 0)
return self
def _grow(self, X, time, event, idx, depth):
node = SurvivalNode()
node.times, node.cumhaz = nelson_aalen(time[idx], event[idx])
if (
depth >= self.max_depth
or idx.size < 2 * self.min_leaf
or event[idx].sum() < 2 * self.min_events
):
return node
features = self.rng_.choice(
self.n_features_, size=self.mtry_, replace=False
)
split = best_split(
X[idx],
time[idx],
event[idx],
features,
n_cuts=self.n_cuts,
min_leaf=self.min_leaf,
min_events=self.min_events,
)
if split is None:
return node
node.feature, node.threshold, _ = split
left = X[idx, node.feature] <= node.threshold
node.left = self._grow(X, time, event, idx[left], depth + 1)
node.right = self._grow(X, time, event, idx[~left], depth + 1)
return nodeThree direct stopping checks are in play: the depth limit, fewer than twice the minimum leaf size, and fewer than twice the minimum event count. A node also stops when no valid candidate split is found. Each accepted split enforces both child constraints, but the root can still be terminal with fewer events than min_events. In a bootstrap tree, repeated copies of an observation count repeatedly toward these constraints, so the counts are not numbers of distinct original patients. The feature subset is drawn freshly at each split. These controls limit the fitted partition; they do not guarantee low variance or calibrated leaf estimates.
Predictions use the right-continuous step_function from the worked example and a function that walks a tree for every row. The earlier helper returns the value after the jump at an event time, matching the Nelson–Aalen convention.
PythonIn [23]Expand
def tree_cumhaz(node, X, grid):
"""Cumulative hazard on a grid for every row of X, from a single tree."""
out = np.empty((X.shape[0], grid.size))
for i in range(X.shape[0]):
current = node
while not current.is_leaf:
if X[i, current.feature] <= current.threshold:
current = current.left
else:
current = current.right
out[i] = step_function(current.times, current.cumhaz, grid)
return outThe first tree is a useful object to inspect before the aggregation blurs it. The figure below shows the predicted cumulative hazard at a fixed time over a slice of the covariate space, in which the three remaining features are held at zero. The blocks are the terminal nodes.
PythonIn [24]Expand
first_tree = SurvivalTree(
max_depth=4, min_leaf=25, min_events=12, n_cuts=8, mtry=2, random_state=17
).fit(X_train, t_train, e_train)
grid_1 = np.linspace(
min(np.quantile(X_train[:, 0], 0.01), X_train[:150, 0].min()),
max(np.quantile(X_train[:, 0], 0.99), X_train[:150, 0].max()),
160,
)
grid_2 = np.linspace(
min(np.quantile(X_train[:, 1], 0.01), X_train[:150, 1].min()),
max(np.quantile(X_train[:, 1], 0.99), X_train[:150, 1].max()),
160,
)
grid_1_pad = 0.03 * np.ptp(grid_1)
grid_2_pad = 0.03 * np.ptp(grid_2)
grid_1 = np.linspace(grid_1[0] - grid_1_pad, grid_1[-1] + grid_1_pad, 160)
grid_2 = np.linspace(grid_2[0] - grid_2_pad, grid_2[-1] + grid_2_pad, 160)
G1, G2 = np.meshgrid(grid_1, grid_2)
slice_rows = np.zeros((G1.size, 5))
slice_rows[:, 0] = G1.ravel()
slice_rows[:, 1] = G2.ravel()
hazard_slice = tree_cumhaz(first_tree.root_, slice_rows, np.array([2.0]))[:, 0]PythonIn [25]Expand
import matplotlib.pyplot as plt
from book_plot_style import PALETTE, polish_axes, use_book_style
use_book_style()
plt.figure(figsize=(6.0, 4.0))
ax = plt.gca()
mesh = ax.pcolormesh(
G1, G2, hazard_slice.reshape(G1.shape), shading="nearest", cmap="viridis"
)
ax.scatter(
X_train[:150, 0],
X_train[:150, 1],
s=6,
color="white",
alpha=0.5,
linewidths=0,
)
ax.set_xlabel("$x_1$")
ax.set_ylabel("$x_2$")
ax.set_title("A shallow tree carves the covariate space into rectangles")
plt.colorbar(mesh, ax=ax, label="estimated cumulative hazard at $t = 2$")
ax.set_xlim(grid_1[0], grid_1[-1])
ax.set_ylim(grid_2[0], grid_2[-1])
ax.grid(False)
plt.show()VisualizationOut [25]Expand
The blocks are axis-aligned, and the prediction is constant inside each block. That is the structural reason the tree can represent interactions and thresholds without being told about them. This threshold sensitivity is one reason a single tree can have high sampling variance: a small change in the data can move a threshold, and moving a threshold can change a whole block of predictions.
From trees to a forest
The forest is a bootstrap loop, the aggregation, and the out-of-bag machinery. The staged_* methods are there so that we can look at the effect of the number of trees: they return one row per tree, which lets us compute the ensemble and the out-of-bag estimate for any prefix of the forest.
PythonIn [26]Expand
class RandomSurvivalForest:
"""Random survival forest with Nelson-Aalen terminal nodes."""
def __init__(
self,
n_trees=100,
max_depth=7,
min_leaf=15,
min_events=8,
n_cuts=10,
mtry=None,
random_state=0,
):
self.n_trees = n_trees
self.max_depth = max_depth
self.min_leaf = min_leaf
self.min_events = min_events
self.n_cuts = n_cuts
self.mtry = mtry
self.random_state = random_state
def fit(self, X, time, event):
X = np.asarray(X, dtype=float)
time = np.asarray(time, dtype=float)
event = np.asarray(event, dtype=bool)
n = X.shape[0]
self.grid_ = np.unique(time[event])
rng = np.random.default_rng(self.random_state)
self.trees_ = []
self.in_bag_ = np.zeros((self.n_trees, n), dtype=bool)
for b in range(self.n_trees):
idx = rng.integers(0, n, size=n)
self.in_bag_[b, idx] = True
tree = SurvivalTree(
max_depth=self.max_depth,
min_leaf=self.min_leaf,
min_events=self.min_events,
n_cuts=self.n_cuts,
mtry=self.mtry,
random_state=int(rng.integers(2**31)),
)
tree.fit(X[idx], time[idx], event[idx])
self.trees_.append(tree)
return self
def predict_cumhaz(self, X):
X = np.asarray(X, dtype=float)
total = np.zeros((X.shape[0], self.grid_.size))
for tree in self.trees_:
total += tree_cumhaz(tree.root_, X, self.grid_)
return total / len(self.trees_)
def predict(self, X):
"""Risk score: the average cumulative hazard summed over the time grid."""
return self.predict_cumhaz(X).sum(axis=1)
def predict_survival(self, X, times):
cumhaz = self.predict_cumhaz(X)
index = np.searchsorted(self.grid_, times, side="right") - 1
evaluated = np.zeros((cumhaz.shape[0], np.asarray(times).size))
observed = index >= 0
evaluated[:, observed] = cumhaz[:, index[observed]]
return np.exp(-evaluated)
def staged_risk_scores(self, X):
"""One row per tree, holding that tree's risk score for every row of X."""
X = np.asarray(X, dtype=float)
return np.array(
[
tree_cumhaz(tree.root_, X, self.grid_).sum(axis=1)
for tree in self.trees_
]
)
def staged_oob_scores(self, X):
"""One row per tree, with NaN wherever the row was in that tree's bootstrap sample."""
X = np.asarray(X, dtype=float)
scores = np.full((self.n_trees, X.shape[0]), np.nan)
for b, tree in enumerate(self.trees_):
out_of_bag = ~self.in_bag_[b]
if out_of_bag.any():
scores[b, out_of_bag] = tree_cumhaz(
tree.root_, X[out_of_bag], self.grid_
).sum(axis=1)
return scoresThe predict_cumhaz method is the average over trees, and the predict method is its sum over the grid, so the risk score inherits the averaging convention: it is the sum of the averaged cumulative hazards. The survival function is the exponential of the negative averaged cumulative hazard, evaluated at the requested times. The staged_oob_scores method stores a NaN in the in-bag positions, which is the matrix form of the condition that observation only receives contributions from trees that did not see it.
Each tree is grown on a bootstrap sample, so its cumulative hazard is defined on its own event times. To combine the trees we evaluate every staircase on the shared grid self.grid_, the distinct event times of the full training sample, before averaging. The step_function used by tree_cumhaz is right-continuous, so the value at a tree's event time is the value just after the jump; between the tree's event times its step function is flat, which is the Nelson–Aalen convention extended off the event times. Averaging these piecewise-constant functions gives the ensemble cumulative hazard, and exponentiating the negative gives the ensemble survival function.
We now fit the forest and compare it with the reference implementation and with Cox regression. We use Harrell's concordance index: the fraction of comparable observed pairs whose risk ordering agrees with event-time ordering, with half credit for tied risks. Its target can depend on censoring, and it is not an unconditional concordance probability without additional assumptions. We implement it explicitly and check it against the reference implementation, because a metric implemented twice is a metric that is less likely to be silently wrong.
PythonIn [27]Expand
from sksurv.ensemble import RandomSurvivalForest as SKRSF
from sksurv.linear_model import CoxPHSurvivalAnalysis
from sksurv.metrics import concordance_index_censored, integrated_brier_score
from sksurv.util import Surv
def harrell_c(time, event, risk):
"""Harrell's concordance index: higher risk must mean shorter survival."""
time = np.asarray(time, dtype=float)
event = np.asarray(event, dtype=bool)
risk = np.asarray(risk, dtype=float)
comparable = event[:, None] & (
(time[:, None] < time[None, :])
| ((time[:, None] == time[None, :]) & ~event[None, :])
)
n_comparable = int(comparable.sum())
if n_comparable == 0:
raise ValueError("No comparable pairs for concordance")
difference = risk[:, None] - risk[None, :]
tied = (np.abs(difference) <= 1e-8) & comparable
ordered = (difference > 1e-8) & comparable
return float((ordered.sum() + 0.5 * tied.sum()) / n_comparable)PythonIn [28]Expand
forest = RandomSurvivalForest(
n_trees=60,
max_depth=6,
min_leaf=15,
min_events=8,
n_cuts=8,
mtry=2,
random_state=0,
).fit(X_train, t_train, e_train)
y_train = Surv.from_arrays(event=e_train, time=t_train)
y_test = Surv.from_arrays(event=e_test, time=t_test)
sk_forest = SKRSF(
n_estimators=60,
min_samples_split=20,
min_samples_leaf=10,
max_features="sqrt",
random_state=0,
n_jobs=-1,
).fit(X_train, y_train)
cox = CoxPHSurvivalAnalysis().fit(X_train, y_train)
risk_forest = forest.predict(X_test)
risk_sk = sk_forest.predict(X_test)
risk_cox = cox.predict(X_test)The next diagnostic varies the number of trees. We plot both OOB and held-out concordance to illustrate their sampling variation; neither curve is required to rise monotonically or to coincide. These curves are not used to tune the settings in this example. Choosing a forest size or other settings from the held-out test curve would consume its independence as a final assessment.
PythonIn [29]Expand
oob_scores = forest.staged_oob_scores(X_train)
test_scores = forest.staged_risk_scores(X_test)
def oob_risk(scores, n_trees):
"""Out-of-bag risk score using the first n_trees trees; NaN where no tree is out of bag."""
sub = scores[:n_trees]
usable = np.isfinite(sub).any(axis=0)
risk = np.full(scores.shape[1], np.nan)
risk[usable] = np.nanmean(sub[:, usable], axis=0)
return risk
tree_counts = np.arange(2, forest.n_trees + 1)
oob_curve = []
test_curve = []
for b in tree_counts:
risk_oob = oob_risk(oob_scores, b)
ok = np.isfinite(risk_oob)
oob_curve.append(harrell_c(t_train[ok], e_train[ok], risk_oob[ok]))
test_curve.append(harrell_c(t_test, e_test, test_scores[:b].mean(axis=0)))
oob_curve = np.array(oob_curve)
test_curve = np.array(test_curve)
oracle_c = harrell_c(t_test, e_test, eta_test)PythonIn [30]Expand
import matplotlib.pyplot as plt
from book_plot_style import PALETTE, polish_axes, use_book_style
use_book_style()
plt.figure(figsize=(6.0, 4.0))
ax = plt.gca()
ax.plot(
tree_counts,
oob_curve,
color=PALETTE["primary"],
linewidth=1.1,
label="out-of-bag",
)
ax.plot(
tree_counts,
test_curve,
color=PALETTE["amber"],
linewidth=1.1,
label="held-out test",
)
ax.axhline(
oracle_c,
color=PALETTE["muted"],
linewidth=1.0,
linestyle="--",
label="oracle (true log hazard)",
)
ax.set_xlabel("number of trees")
ax.set_ylabel("concordance index")
ax.set_title("The out-of-bag curve tracks the test curve")
ax.legend(loc="lower right")
polish_axes(ax)
plt.show()VisualizationOut [30]Expand
The two curves fluctuate as trees are added, and their later values summarize this particular fit and sample. A gap to the true-signal reference can reflect finite data, the chosen split grid and leaf constraints, fitting randomness, and sampling variation. It does not show that a forest could never recover the underlying smooth relationship.
The comparison table also reports integrated Brier score. At each time, it compares the predicted survival probability with the survival indicator , estimating squared error with inverse-probability-of-censoring weights. The reference implementation estimates a marginal censoring distribution from training data. Its validity needs independence appropriate to those marginal weights and positive censoring-survival probabilities over the evaluation interval. If censoring depends on covariates, conditional weights may be needed even when event and censoring times are independent given those covariates. Under those conditions the Brier score assesses probability accuracy, including calibration, rather than ranking alone. We select the time grid from training follow-up before scoring test outcomes.
PythonIn [31]Expand
score_times = np.quantile(t_train, np.linspace(0.15, 0.85, 25))
curve_forest = forest.predict_survival(X_test, score_times)
curve_sk = np.vstack(
[fn(score_times) for fn in sk_forest.predict_survival_function(X_test)]
)
curve_cox = np.vstack(
[fn(score_times) for fn in cox.predict_survival_function(X_test)]
)
curve_oracle = np.exp(
-0.20 * score_times[None, :] ** 1.5 * np.exp(eta_test)[:, None]
)| test C-index | test integrated Brier | |
|---|---|---|
| model | ||
| oracle (true log hazard) | 0.863 | 0.097 |
| Cox PH, linear | 0.722 | 0.161 |
| random survival forest, ours | 0.832 | 0.119 |
| random survival forest, scikit-survival | 0.840 | 0.111 |
ConsoleOut [33]Expand
our C-index: 0.8324 scikit-survival C-index for the same risk scores: 0.8324 difference: 0.00e+00
The two concordance calculations agree on this continuous-time sample to the reported tolerance. The helper also specifies how event–censor ties, near-tied risks, and no-comparable-pair inputs are handled; agreement on one sample is not exhaustive validation. This is one seeded train/test comparison with fixed settings, not a population ranking. Both forests use 60 trees, but their split searches and leaf constraints differ, and their survival-aggregation conventions differ. The simulated log-hazard predictor is additive in transformed features but nonlinear in the original and . A Cox model supplied with the sine and quadratic transformations could represent it; the raw-feature linear Cox model cannot. The oracle is a known-signal reference, not a guaranteed finite-sample winner.
Risk-score profiles show how the fitted responses change along one feature. For each grid value, we replace that feature in every background row and average the predictions. The plots include the known log-hazard component as a reference and standardize each curve separately to reveal turning points. Forest mortality scores and Cox log-risk scores use different scales, so these profiles do not measure comparable effect magnitudes, causal effects, or prediction errors.
PythonIn [34]Expand
def partial_dependence(model, background, feature, grid):
"""Average model prediction while a single feature is set to each value of the grid."""
values = []
for value in grid:
rows = background.copy()
rows[:, feature] = value
values.append(model.predict(rows).mean())
return np.array(values)
def zscore(values):
values = np.asarray(values, dtype=float)
scale = values.std()
return (
(values - values.mean()) / scale if scale > 0 else np.zeros_like(values)
)
x1_grid = np.linspace(
np.quantile(X_train[:, 0], 0.02), np.quantile(X_train[:, 0], 0.98), 25
)
x2_grid = np.linspace(
np.quantile(X_train[:, 1], 0.02), np.quantile(X_train[:, 1], 0.98), 25
)
background = X_train[:60]
curves = {
"forest": (
partial_dependence(forest, background, 0, x1_grid),
partial_dependence(forest, background, 1, x2_grid),
),
"Cox PH": (
partial_dependence(cox, background, 0, x1_grid),
partial_dependence(cox, background, 1, x2_grid),
),
}
truth = {
"x1": 1.2 * x1_grid + 1.2 * np.sin(2.0 * x1_grid),
"x2": 1.2 * (x2_grid**2 - 1.0),
}
pd_plot_grids = (x1_grid, x2_grid)
pd_plot_curves = {
"truth": (zscore(truth["x1"]), zscore(truth["x2"])),
**{
name: tuple(zscore(value) for value in pair)
for name, pair in curves.items()
},
}VisualizationOut [35]Expand
VisualizationOut [36]Expand
The first profile shows bends in the forest response that the Cox model's linear log-risk response cannot express. The second shows a U-shaped forest response along the quadratic feature. The fitted Cox coefficient for that feature is small but nonzero in this finite sample, and separate standardization stretches it into a visible straight line. A standardized profile therefore reveals shape while concealing relative magnitude. Use the held-out concordance and integrated Brier scores for quantitative prediction comparisons.
The final diagnostic is calibration. The forest predicts a survival curve per patient, so we can group patients by their predicted risk and compare the average predicted curve with the Kaplan–Meier curve of the observed data in that group.
This compares each group's mean predicted survival with its Kaplan–Meier estimate. In this simulation censoring is independent of the event time and covariates, which makes Kaplan–Meier an appropriate comparator within the predicted-risk groups. In an application, this requires independent censoring within each group. Agreement of group averages is one aggregate check and does not establish calibration for every covariate profile. Sampling uncertainty also matters, especially late in follow-up.
PythonIn [37]Expand
def kaplan_meier(time, event):
"""Kaplan-Meier product-limit estimate: distinct event times and survival values."""
time = np.asarray(time, dtype=float)
event = np.asarray(event, dtype=bool)
n = time.size
order = np.argsort(time, kind="mergesort")
t_sorted, e_sorted = time[order], event[order]
unique_times, first = np.unique(t_sorted, return_index=True)
at_risk = n - first
deaths = np.add.reduceat(e_sorted.astype(float), first)
keep = deaths > 0
return unique_times[keep], np.cumprod(1.0 - deaths[keep] / at_risk[keep])
cumhaz_test = forest.predict_cumhaz(X_test)
risk_test = cumhaz_test.sum(axis=1)
low_risk, high_risk = (
risk_test <= np.quantile(risk_test, 0.3),
risk_test >= np.quantile(risk_test, 0.7),
)
calibration_curves = []
for name, mask in [("lowest 30%", low_risk), ("highest 30%", high_risk)]:
km_time, km_surv = kaplan_meier(t_test[mask], e_test[mask])
calibration_curves.append(
(
name,
np.r_[0.0, forest.grid_],
np.r_[1.0, np.exp(-cumhaz_test[mask]).mean(axis=0)],
np.r_[0.0, km_time],
np.r_[1.0, km_surv],
)
)PythonIn [38]Expand
import matplotlib.pyplot as plt
from book_plot_style import PALETTE, polish_axes, use_book_style
use_book_style()
plt.figure(figsize=(6.0, 4.8))
ax = plt.gca()
for (name, pred_times, pred_surv, km_time, km_surv), colour in zip(
calibration_curves, [PALETTE["primary"], PALETTE["coral"]]
):
ax.step(
pred_times,
pred_surv,
where="post",
color=colour,
linewidth=1.1,
label=f"{name}: forest",
)
ax.step(
km_time,
km_surv,
where="post",
color=colour,
linewidth=1.1,
linestyle="--",
label=f"{name}: Kaplan–Meier",
)
ax.set_xlabel("time")
ax.set_ylabel("survival probability")
ax.set_xlim(0.0, 6.0)
ax.set_ylim(0.0, 1.02)
ax.set_title("Predicted and observed survival by risk group")
ax.legend(loc="upper center", bbox_to_anchor=(0.5, -0.18), ncol=2)
polish_axes(ax)
plt.show()VisualizationOut [38]Expand
The low-risk group has higher survival than the high-risk group, but the predicted and observed curves do not coincide. The low-risk mean prediction lies below its Kaplan–Meier curve through much of follow-up, and late high-risk estimates also differ. These gaps motivate further calibration assessment; this single grouped plot is not a calibration certificate. Sparse late risk sets add sampling uncertainty, so a gap alone does not identify its cause.
We repeat the comparison with a linear log-hazard predictor. This supplies a case in which the raw-feature Cox model is correctly specified. Its structural restriction can then help estimation, but it does not guarantee a win or tie in every finite sample, metric, or forest configuration. Read the actual result below as another seeded illustration.
PythonIn [39]Expand
rng_linear = np.random.default_rng(4242)
Xl_train, tl_train, el_train, etal_train = make_survival_data(
500, rng_linear, nonlinear=False
)
Xl_test, tl_test, el_test, etal_test = make_survival_data(
400, rng_linear, nonlinear=False
)
yl_train = Surv.from_arrays(event=el_train, time=tl_train)
forest_linear = RandomSurvivalForest(
n_trees=60,
max_depth=6,
min_leaf=15,
min_events=8,
n_cuts=8,
mtry=2,
random_state=0,
).fit(Xl_train, tl_train, el_train)
cox_linear = CoxPHSurvivalAnalysis().fit(Xl_train, yl_train)| test C-index | |
|---|---|
| model | |
| oracle (true log hazard) | 0.804 |
| Cox PH, linear | 0.799 |
| random survival forest, ours | 0.773 |
On this linear-data draw, compare each fitted model's concordance with the known predictor's value. Together the two examples illustrate why the predictor specification matters. They do not isolate a universal efficiency penalty, establish an ordering for every dataset, or validate the full survival curves from concordance alone.
Key Parameters
The number of trees (n_estimators in scikit-survival, n_trees here) controls averaging and compute cost. Additional trees can reduce randomization variation with diminishing returns. The examples use 60 trees to keep the scratch calculation manageable; that count is not a validated default for another dataset. Assess stabilization with OOB or validation predictions rather than repeatedly selecting from a final test set.
Our min_events constrains event counts in the children of an accepted split, while min_leaf constrains their observation counts. These controls affect the amount of information available to estimate each staircase. A small late risk set can still create a large jump. Scikit-survival's min_samples_leaf and min_samples_split are observation-count constraints and do not translate into an exact event-count minimum. The settings here are illustrative; select them with the event count, follow-up support, and validation behavior in mind.
The number of features considered at each split (max_features in scikit-survival, mtry here) controls the correlation between the trees. "sqrt" is the reference default and "log2" is another option. Their rounded feature counts depend on the number of features; "log2" is not always the smaller subset, and a value equal to the number of features removes feature-subset randomness. Bootstrap randomness remains, and implementations may also randomize ties or feature order.
The depth (max_depth) and the minimum leaf size (min_samples_leaf) control the bias-variance trade-off. Allowing deeper trees enlarges the available partition and can reduce approximation bias while increasing sampling sensitivity, and the leaf-size rule acts as a smoother. An unrestricted tree is rarely what we want, because the terminal nodes end up with too few events.
Key Methods
In scikit-survival, fit(X, y) takes a feature matrix and a structured outcome array whose first field is the boolean event indicator and whose second field is the numeric observed time. Surv.from_arrays(event=..., time=...) constructs that array; keyword arguments make the intended roles explicit. The scratch class instead takes fit(X, time, event) and assumes finite numeric features and valid nonnegative observed times.
predict(X) returns the risk score, which is the sum of the cumulative hazard over the ensemble's time grid. It is a rank-order summary, not a probability and not a hazard ratio. Larger values mean higher risk.
predict_cumulative_hazard_function(X) and predict_survival_function(X) return arrays of step functions, one per row. They are called with a time vector, and they return the value of the function at that time. For a single prediction, predict_survival_function(X)[0](times) is the survival function evaluated at times. Our predict_cumhaz returns an array on grid_, and predict_survival evaluates the negative exponential of that averaged hazard at supplied times. The reference survival method uses averaged Kaplan–Meier estimates, so it need not return the same survival probabilities.
predict is not the same as predict_survival_function evaluated at a single time. The former is a summary of the whole curve, and the latter is the curve. If two patients have the same risk score but different curves, the score cannot distinguish them, and the curves are what the clinician needs.
Comments
Add to the conversation
A question, a correction, or an example from your own work.
No comments yet. Be the first to share your thoughts!