The Cox model and its parametric cousins assume a particular shape for how covariates act on the hazard (usually linear and proportional). When you have many covariates, non-linear effects, or interactions you did not anticipate, tree-based machine-learning models can predict better without those assumptions. Greenwood ships three, all fully non-parametric and all producing per-subject survival and cumulative-hazard functions rather than a single risk score: a SurvivalTree, the bagged RandomSurvivalForest (and its extremely-randomized ExtraSurvivalTrees variant), and GradientBoostingSurvivalAnalysis. This page introduces each, then shows how to score, interpret, and visualize them.
We’ll work with lung-dataset outcomes, keeping the rows with complete covariates so the examples are easy to follow. Every model below is fit from the same formula, and covariates lists the predictor columns to pass when predicting.
import greenwood as gw
# Load the lung dataset and keep rows with complete covariates
lung = gw.load_dataset("lung", backend="polars").drop_nulls(
subset=["ph.ecog", "ph.karno", "wt.loss"]
)
# Define the model formula, the survival response, and the covariate list
formula = "Surv(time, status) ~ age + sex + ph.ecog + ph.karno + wt.loss"
death = gw.Outcome.surv(time="time", event="status")
covariates = ["age", "sex", "ph.ecog", "ph.karno", "wt.loss"]
# Bind the response to the data to display its summary
death.bind(lung)
Surv(type=right, n=213, events=151)
Survival trees
A SurvivalTree recursively splits subjects into groups whose survival curves differ as much as possible, measured by the log-rank statistic at each candidate split. Each leaf stores a Kaplan-Meier survival curve and a Nelson-Aalen cumulative hazard estimated from the subjects that land in it. A single tree is easy to reason about but high-variance; it is most useful as the interpretable building block of an ensemble.
# Grow a shallow, interpretable survival tree
tree = gw.SurvivalTree(max_depth=3, random_state=0).fit(formula, data=lung)
tree
SurvivalTree (log-rank splits, 15 nodes, 8 leaves)
n = 213, events = 151, features = 5
Calling predict with the default type="risk" returns Ishwaran’s ensemble mortality (the sum of the predicted cumulative hazard over the training event times). Larger values mean higher risk.
# Risk score for the first eight subjects (higher = higher risk)
tree.predict(lung[covariates])[:8]
array([ 95.15345288, 95.15345288, 95.15345288, 95.15345288,
95.15345288, 168.54325397, 168.54325397, 95.15345288])
Random survival forests
A RandomSurvivalForest averages many decorrelated survival trees, each grown on a bootstrap sample of the data and considering a random subset of covariates at every split. Averaging drives down the variance of the single tree and yields a well-calibrated, flexible model. Passing oob_score=True computes an out-of-bag concordance estimate (an honest, held-out measure of discrimination) using only the trees for which each subject was not in the bootstrap sample.
# Fit a forest and report the out-of-bag concordance
rsf = gw.RandomSurvivalForest(
n_estimators=200, oob_score=True, random_state=0
).fit(formula, data=lung)
rsf
RandomSurvivalForest (200 log-rank trees, max_features='sqrt')
n = 213, events = 151, features = 5
out-of-bag concordance = 0.6061
Predictions
Like the tree, the forest predicts a scalar risk score by default. It also predicts full survival or cumulative-hazard curves: pass type="survival" (or "cumulative_hazard") and optional times, and you get one column per subject.
# Predicted survival probabilities at one and two years for three subjects
rsf.predict(lung[covariates][:3], type="survival", times=[365, 730])
PolarsRows2Columns4 |
|
|
|
|
|
| 0 |
365 |
0.570687942169 |
0.664342402851 |
0.433306613867 |
| 1 |
730 |
0.0856037841238 |
0.366642713531 |
0.106499220759 |
Variable importance
variable_importance() measures how much each covariate contributes by permuting its values in the out-of-bag samples and recording the drop in concordance. Larger values indicate more important covariates. The same table is what tidy() returns for a fitted forest.
# Permutation variable importance (out-of-bag)
rsf.variable_importance()
PolarsRows5Columns2 |
|
|
|
| 0 |
ph.ecog |
0.0627294117647 |
| 1 |
sex |
0.0198352941176 |
| 2 |
age |
0.00245882352941 |
| 3 |
wt.loss |
-0.00311764705882 |
| 4 |
ph.karno |
-0.00812941176471 |
Gradient-boosted survival
GradientBoostingSurvivalAnalysis takes a different route. Instead of averaging independent trees, it builds an additive model by fitting each new (shallow) regression tree to the negative gradient of the Cox partial-likelihood loss (the martingale residuals under the current model) and adding a shrunken version of it. The result is a non-linear generalization of the Cox model that often predicts very well on structured tabular data.
# Fit a gradient-boosted survival model
gbm = gw.GradientBoostingSurvivalAnalysis(
n_estimators=300, learning_rate=0.05, max_depth=2, random_state=0
).fit(formula, data=lung)
gbm
GradientBoostingSurvivalAnalysis (300 trees, learning_rate=0.05, max_depth=2)
n = 213, events = 151, features = 5
The additive score F(x) is the boosted analogue of the Cox linear predictor. predict(type="lp") returns it, type="risk" returns exp(F(x)), and type="survival" combines it with a Breslow baseline hazard to produce survival curves.
# Survival curves at one and two years for three subjects
gbm.predict(lung[covariates][:3], type="survival", times=[365, 730])
PolarsRows2Columns4 |
|
|
|
|
|
| 0 |
365 |
0.384085259396 |
0.633711204665 |
0.331726873354 |
| 1 |
730 |
0.0432828906933 |
0.22382813027 |
0.0267584698433 |
Boosting reports impurity-based variable importance (the total reduction in squared error credited to each covariate, normalized to sum to one):
# Impurity-based variable importance
gbm.variable_importance()
PolarsRows5Columns2 |
|
|
|
| 0 |
wt.loss |
0.409465130649 |
| 1 |
ph.ecog |
0.212804757602 |
| 2 |
age |
0.205707927573 |
| 3 |
sex |
0.11564191931 |
| 4 |
ph.karno |
0.0563802648662 |
IPC-weighted ridge regression
Not every alternative to Cox needs a tree. IPCRidge fits a penalized linear model on log survival times, but instead of the Cox partial likelihood it uses inverse-probability-of-censoring (IPC) weights to correct for censoring. Only uncensored subjects contribute to the fit, and each is weighted by 1 / G(t-) from the Kaplan-Meier censoring estimate so that subjects whose events occur when censoring is heavy receive greater influence. The result is a simple linear model that does not assume proportional hazards.
# Fit an IPC-weighted ridge regression model
ridge = gw.IPCRidge(alpha=1.0).fit(formula, data=lung)
ridge
IPCRidge (IPC-weighted ridge, alpha=1.0)
coef
(Intercept) 5.315
age -0.01992
sex 0.6572
ph.ecog -0.5318
ph.karno 0.01457
wt.loss 0.001947
n = 213, events = 151 (used for fitting)
The alpha parameter controls the ridge penalty. Setting alpha=0 gives unpenalized IPC-weighted OLS. Larger values shrink the coefficients toward zero, which helps when covariates are correlated or numerous.
Predictions are on the log-time scale by default (type="lp"). Pass type="response" to get predicted survival times on the original scale.
# Predict survival times on the original scale for the first five subjects
ridge.predict(lung[covariates][:5], type="response")
array([386.72519998, 491.17815835, 280.69111386, 509.60077987,
109.31469154])
The coefficient table is available via tidy(), and model-level summaries via glance().
# Coefficient table for the ridge model
gw.tidy(ridge)
PolarsRows6Columns2 |
|
|
|
| 0 |
(Intercept) |
5.31469425325 |
| 1 |
age |
-0.0199243788921 |
| 2 |
sex |
0.657164509677 |
| 3 |
ph.ecog |
-0.531839726515 |
| 4 |
ph.karno |
0.0145723160522 |
| 5 |
wt.loss |
0.00194699468604 |
# Model-level summary statistics
gw.glance(ridge)
PolarsRows1Columns4 |
|
|
|
|
|
| 0 |
213 |
151 |
1 |
5 |
Because the linear predictor is on the log-time scale (higher means longer survival), negate it to get a risk score compatible with concordance_index().
# Negate the linear predictor so higher values mean higher risk
risk = -ridge.predict(lung[covariates], type="lp")
# Compute concordance for the ridge model
gw.concordance_index(death, risk=risk, data=lung)
Visualizing per-subject curves
Because these models predict an individual survival curve for every subject, plot_predicted_survival() draws one step curve per row of new data. It works with any of the machine-learning models (and with CoxPH), so you can compare how predicted survival varies across individuals.
# Compare predicted survival for four patients from the forest
gw.plot_predicted_survival(rsf, newdata=lung[covariates][:4])
Pass type="cumulative_hazard" for the cumulative-hazard scale, or labels= to name the curves.
Comparing models
All of these models expose a risk score, so you can judge and compare them with the same censoring-aware metrics as any other model. Here we compare the three ensembles by concordance on the training data (for an honest estimate, evaluate on held-out data or use the forest’s out-of-bag score).
# In-sample concordance for the three ensemble models
scores = {
"random forest": gw.concordance_index(death, risk=rsf.predict(lung[covariates]), data=lung),
"extra trees": gw.concordance_index(death, risk=ext.predict(lung[covariates]), data=lung),
"gradient boosting": gw.concordance_index(death, risk=gbm.predict(lung[covariates]), data=lung),
}
scores
{'random forest': 0.8101176470588235,
'extra trees': 0.726764705882353,
'gradient boosting': 0.7446764705882353}
See Prediction performance for concordance on held-out data, the Brier score, and time-dependent AUC (the right way to compare these models).
SurvivalBoost: competing-risks gradient boosting
The models above apply to single-event survival data. When you have competing risks, where subjects can experience one of several mutually exclusive event types, SurvivalBoost is the machine-learning counterpart. It estimates cause-specific cumulative incidence functions (CIFs) and the any-event survival function directly, using multi-class softmax gradient boosting with IPCW-weighted targets. The approach follows Alberge et al. (2025).
How it works
At each boosting iteration, SurvivalBoost samples random time horizons and appends them as an extra feature column alongside the original covariates. For every subject-horizon pair, it computes a multi-class target: event type k if cause k occurred before the horizon, or “event-free” if the subject was still at risk. Subjects censored before the horizon receive zero weight and are excluded. Uncensored subjects are reweighted by the inverse probability of censoring (IPCW) to correct for informative dropout.
The model then fits one regression tree per class to the softmax gradient residuals of the current ensemble, exactly as in multi-class gradient boosting for tabular data. Because the time horizon varies from iteration to iteration, the ensemble learns how the class probabilities evolve over time without being tied to a fixed grid. Predictions are probability-calibrated by construction: the softmax output sums to one at every time point.
The hard_zero_fraction parameter forces a fraction of time horizons to exactly zero each iteration. This teaches the model that no events have occurred at t = 0, preventing CIFs from starting at a nonzero value.
Fitting the model
SurvivalBoost consumes a multi-state Surv response (from gw.Surv() with a categorical event, or from gw.first_event()) and a covariate matrix. We use simulate_competing_risks() to generate synthetic data with known properties, but the model works with any competing-risks dataset.
import numpy as np
# Generate competing-risks data: 2 causes, 3 covariates
sim = gw.simulate_competing_risks(n=500, n_causes=2, n_covariates=3, seed=42)
# Fit SurvivalBoost
sb = gw.SurvivalBoost(
n_estimators=50, learning_rate=0.1, max_depth=3, random_state=0
)
sb.fit(sim.surv, covariates=sim.covariates)
sb
SurvivalBoost (50 rounds, learning_rate=0.1, max_depth=3)
n = 500, events = 296, causes = 2, features = 3
The fitted model reports the number of subjects, events, causes, and features. Internally it stores one set of trees per class per boosting round.
Predicting CIFs and survival
predict_cumulative_incidence returns a 3-D array of shape (n_subjects, n_classes, n_times). Index 0 along the class axis is the survival probability (the probability that no event has occurred by time t). Indices 1 through K are the cause-specific CIFs.
# Predict CIFs at four time points; result is (n_subjects, n_classes, n_times)
times = np.array([2.0, 5.0, 10.0, 15.0])
cif = sb.predict_cumulative_incidence(sim.covariates, times=times)
cif.shape
Because the predictions pass through a softmax, the probabilities always sum to one across classes at every time point. That means the survival probability is exactly 1 - sum(CIFs), which is the correct accounting for competing risks.
# Verify: probabilities sum to 1 at each time for every subject
cif.sum(axis=1)[:3, :]
array([[1., 1., 1., 1.],
[1., 1., 1., 1.],
[1., 1., 1., 1.]])
The survival function alone is available via predict_survival_function, which returns a 2-D array of shape (n_subjects, n_times).
# Extract the survival function alone as a 2-D array
surv_pred = sb.predict_survival_function(sim.covariates, times=times)
surv_pred[:3, :] # first three subjects
array([[0.96641157, 0.35845857, 0.05182967, 0.02070598],
[0.96287637, 0.31153838, 0.05339982, 0.02709674],
[0.97245538, 0.37588905, 0.04828712, 0.02209601]])
For a classifier-style interface at a single time horizon, predict_proba returns a 2-D array where each row sums to one.
# Class probabilities at t = 8 for the first five subjects
sb.predict_proba(sim.covariates[:5], time_horizon=8.0)
array([[0.08437444, 0.63185384, 0.28377172],
[0.08310677, 0.36864524, 0.54824798],
[0.08290381, 0.6172399 , 0.29985629],
[0.11228157, 0.47509787, 0.41262056],
[0.1260996 , 0.6473551 , 0.2265453 ]])
DataFrame output
Both prediction methods accept a format= parameter for DataFrame output. For predict_cumulative_incidence, the result is a long-format table with columns time, cause, subject, and probability.
# Long-format CIF DataFrame for three subjects at two time points
sb.predict_cumulative_incidence(
sim.covariates[:3], times=np.array([5.0, 10.0]), format="polars"
).head(6)
PolarsRows6Columns4 |
|
|
|
|
|
| 0 |
5 |
survival |
1 |
0.358458568025 |
| 1 |
5 |
survival |
2 |
0.31153837967 |
| 2 |
5 |
survival |
3 |
0.375889051292 |
| 3 |
5 |
cause1 |
1 |
0.366668870635 |
| 4 |
5 |
cause1 |
2 |
0.291875374158 |
| 5 |
5 |
cause1 |
3 |
0.377649004977 |
For predict_survival_function, the result is a wide-format table with one column per subject, matching the format used by GradientBoostingSurvivalAnalysis.
# Wide-format survival DataFrame with one column per subject
sb.predict_survival_function(
sim.covariates[:3], times=np.array([5.0, 10.0]), format="polars"
)
PolarsRows2Columns4 |
|
|
|
|
|
| 0 |
5 |
0.358458568025 |
0.31153837967 |
0.375889051292 |
| 1 |
10 |
0.051829671427 |
0.0533998235358 |
0.0482871222538 |
Evaluating predictions
Because SurvivalBoost produces CIF estimates, you can evaluate them with the competing-risks metrics from Prediction performance. The cause-specific Brier score measures calibration at each time horizon. Lower values mean the predicted probabilities are closer to the observed outcomes.
# Brier score for cause 1 at several time horizons
eval_times = np.array([3.0, 6.0, 10.0])
cif_eval = sb.predict_cumulative_incidence(sim.covariates, times=eval_times)
bs = gw.brier_score_incidence(sim.surv, incidence_prob=cif_eval[:, 1, :], times=eval_times, cause=1)
dict(zip(eval_times, bs))
{np.float64(3.0): np.float64(0.21736902139521821),
np.float64(6.0): np.float64(0.21389327807050323),
np.float64(10.0): np.float64(0.20768231209439034)}
The concordance index for incidence measures discrimination: how well the model ranks subjects by their predicted CIF at a given time horizon. Values above 0.5 indicate better-than-random ranking.
# Concordance for cause 1 at median follow-up
tau = float(np.median(sim.surv.stop))
cif_tau = sb.predict_cumulative_incidence(sim.covariates, times=np.array([tau]))
c_index = gw.concordance_index_incidence(
sim.surv, incidence_prob=cif_tau[:, 1, 0], cause=1, tau=tau
)
c_index
np.float64(0.5983861205434993)
You can also check calibration (how close the predicted CIFs are to the Aalen-Johansen reference) and accuracy in time (classification accuracy at each horizon). See Prediction performance for details on all four metrics.
Model summaries
SurvivalBoost supports the broom-style tidy() and glance() interface. glance() returns a one-row summary with the model’s key parameters and counts.
# One-row summary of the fitted SurvivalBoost model
gw.glance(sb)
PolarsRows1Columns7 |
|
|
|
|
|
|
|
|
| 0 |
500 |
296 |
2 |
50 |
0.1 |
3 |
3 |
Tuning
The main hyperparameters to tune are:
n_estimators: More boosting iterations improve fit up to a point, after which the model overfits. Typical values range from 50 to 500.learning_rate: Smaller rates (0.01 to 0.1) need more iterations but generalize better. Together, n_estimators and learning_rate control the bias-variance tradeoff.max_depth: Controls the interaction depth of each tree. Depth 2 to 4 is typical. Deeper trees capture higher-order interactions but are more prone to overfitting.subsample: Setting this below 1.0 gives stochastic gradient boosting, where each iteration uses a random fraction of the training data. Values around 0.8 often improve generalization.n_horizons_per_observation: The number of time horizons sampled per subject per iteration. Higher values (3 to 5) give more stable gradients at the cost of proportionally slower training.hard_zero_fraction: Ensures the model learns zero incidence at t = 0. The default of 0.1 works well in most settings. Increase it if predicted CIFs at t = 0 are not close to zero.
Use held-out data or cross_validate to select these. Because the model is non-parametric, it can capture complex covariate effects that Fine-Gray would miss, but it also needs enough data and regularization (depth limits, learning-rate shrinkage) to avoid overfitting.
Fine-Gray fits a single cause at a time and assumes a log-linear effect on the subdistribution hazard. SurvivalBoost estimates all causes jointly, captures non-linear effects and interactions automatically, and produces probability-calibrated predictions without post-processing. The tradeoff is interpretability: Fine-Gray gives a coefficient per covariate with a clear hazard-ratio interpretation, while SurvivalBoost is a black-box ensemble. When prediction accuracy matters more than interpretability, SurvivalBoost is the stronger choice.
More than two causes
SurvivalBoost handles any number of competing event types. The number of output classes is always K + 1 (one per cause, plus the survival class).
# Simulate data with three competing causes and fit SurvivalBoost
sim3 = gw.simulate_competing_risks(n=400, n_causes=3, n_covariates=2, seed=7)
sb3 = gw.SurvivalBoost(n_estimators=30, learning_rate=0.1, max_depth=2, random_state=0)
sb3.fit(sim3.surv, covariates=sim3.covariates)
# CIF at t = 5: shape is (400, 4, 1) because there are 4 classes (survival + 3 causes)
cif3 = sb3.predict_cumulative_incidence(sim3.covariates, times=np.array([5.0]))
cif3.shape
The model scales naturally: each additional cause adds one class to the softmax and one tree per boosting round. The probabilities still sum to one across all classes at every time point, so the survival function and CIFs remain mutually consistent regardless of how many causes are present.
For background on how competing risks are represented and why the cumulative incidence function is the right estimand, see Competing risks. For evaluating CIF predictions with discrimination, calibration, and classification metrics, see Prediction performance.