Machine Learning (beta)

Tree-based survival models: a survival tree, and a random survival forest built from them. Both accept the full surpyval data model – arbitrary censoring and truncation – and return a fitted model at each leaf (parametric, or a Nelson-Aalen or Turnbull estimate for kind="non-parametric"), so a prediction is a distribution rather than a point.

The tree kind couples the split with the leaf model, and every kind takes every kind of censoring and truncation:

kind

Observed and right-censored data, with or without left truncation

Left- or interval-censored data, or right truncation (any truncation)

"non-parametric"

Risk-set log-rank split; Nelson-Aalen leaves

Turnbull-score split (log-rank scores of the pooled Turnbull estimate, each less the score of its truncation window); Turnbull leaves

"exponential"

Exponential deviance split; Exponential leaves

The same, on the full likelihood

"weibull" (default)

Weibull deviance split; Weibull leaves

The same, on the full likelihood

For "non-parametric" the column is chosen at each node (and leaf) from the rows that reach it. Either selection ("greedy" or "ctree") works with every kind and data type.

Warning

These live under surpyval.beta because their API is not yet settled and may change between minor versions without the deprecation cycle the rest of the package follows. They are tested and usable; they are not covered by the same stability promise.

For a narrative introduction see the survival-forest section of Regression Modelling with SurPyval.

A few things worth knowing:

  • Fitted from a DataFrame (fit_from_df with Z_cols or a formula, or fit with a DataFrame Z), a tree or forest keeps the covariate names as feature_names; print(tree) shows its splits by name (temp <= 42), feature_importances is a pandas.Series keyed by name, and predictions read a DataFrame by those names. Fitted from arrays the covariates are shown as Z0, Z1, …

  • A "weibull" or "exponential" tree grows until min_leaf_samples or min_leaf_failures stops it, which suits a forest. For a tree used on its own, set min_split_gain="aic" (or "bic", or a log-likelihood gain), or use selection="ctree".

  • On observed and right-censored data the likelihood splits are found directly (no optimiser per candidate), at a cost of the same order as the log-rank split’s; with left or interval censoring or truncation each candidate needs an optimiser and growing a forest takes much longer. A parametric leaf is fitted when the tree first predicts, so the first prediction of a "weibull" forest can take longer than growing it.

  • A forest grows its trees one after another by default. n_jobs (as in joblib and scikit-learn; -1 for every core) grows them in worker processes; given a random_state the forest is the same whatever n_jobs is.

Random Survival Forest

class surpyval.beta.ml.forest.forest.RandomSurvivalForest(data: SurpyvalData, Z: ArrayLike | NDArray, n_trees: int = 100, max_depth: int | float = inf, min_leaf_samples: int = 5, min_leaf_failures: int = 2, n_features_split: int | float | str = 'sqrt', bootstrap: bool = True, kind: str = 'weibull', selection: str = 'greedy', alpha_split: float = 0.05, random_state: Any = None, feature_names: list[str] | None = None, min_split_gain: float | str = 0.0, n_jobs: int = 1)

Bases: RegressionDataFrameMixin, SerialisableMixin

Random survival forest: an ensemble of survival trees.

n_trees instances of SurvivalTree are fitted, each to an independently bootstrapped sample of the data and each considering a random subset of the features at every split. A prediction evaluates all of them and averages the fitted models their leaves return, which trades the variance of one deep tree for the bias of an average.

Constructed by fit() rather than directly. A fitted forest keeps its training data (data and Z, after dropping rows with a missing covariate) and bootstrap_indices, the rows of data each tree was grown on (with repeats); oob_log_likelihood() and feature_importances() use them to score every row with the trees that did not see it.

A forest fitted from a DataFrame (fit_from_df, or fit with a DataFrame Z) keeps the covariate names as feature_names (None when fitted from an array); its trees’ split descriptions and feature_importances() use them, and its predictions accept a DataFrame Z, read by those names.

Hf(x: int | float | ArrayLike, Z: ArrayLike | NDArray) → NDArray

Cumulative hazard averaged over the trees, as for sf().

df(x: int | float | ArrayLike, Z: ArrayLike | NDArray) → NDArray

Density averaged over the trees, as for sf().

feature_importances(n_repeats: int = 5, random_state: Any = None) → NDArray

Permutation importance of each feature, by out-of-bag log-likelihood.

For each feature, its values are shuffled among each tree’s out-of-bag rows (Breiman, 2001), which breaks its link with the outcome while keeping its distribution, and the out-of-bag log-likelihood (see oob_log_likelihood()) is computed again. A feature’s importance is the drop from the unshuffled value, averaged over n_repeats shuffles: about zero for a feature the forest does not use, positive for one it relies on, and on the scale of the log-likelihood per observation.

Each drop is the mean over the same rows before and after the shuffle, the rows whose out-of-bag log-likelihood is finite in both. A row the ensemble gives zero probability (see oob_log_likelihood()) has a log-likelihood of -inf, and a drop from or to -inf is no number at all: such a row is left out (of every feature’s importance if it is -inf unshuffled, of that shuffle’s drop if it becomes -inf only when shuffled), with one warning giving the counts. Grow more trees, or use kind="exponential", so that every row is scored.

Parameters:
  • n_repeats (int, optional) – The number of shuffles averaged for each feature. Defaults to 5.

  • random_state (None, int or numpy.random.Generator, optional) – Seeds the shuffles. None (the default) draws from NumPy’s global random state, so np.random.seed reproduces it; a seed or Generator gives its own stream.

Returns:

One importance per column of Z, indexed by the feature names (feature_labels: feature_names for a forest fitted from a DataFrame, Z0, Z1, … otherwise).

Return type:

pandas.Series

Raises:

ValueError – If n_repeats is not a positive integer, or on a forest restored with from_dict (see oob_log_likelihood()).

Examples

Only the first feature matters; the second is noise:

>>> import numpy as np
>>> from surpyval.beta.ml import RandomSurvivalForest
>>> rng = np.random.default_rng(0)
>>> Z = rng.uniform(0, 1, (200, 2))
>>> x = rng.weibull(2.0, 200) * np.where(Z[:, 0] > 0.5, 5.0, 10.0)
>>> np.random.seed(0)
>>> forest = RandomSurvivalForest.fit(
...     x, Z, n_trees=20, max_depth=1, kind="exponential"
... )
>>> forest.feature_importances(random_state=1).round(3)
Z0    0.093
Z1   -0.001
Name: importance, dtype: float64
property feature_labels: list[str]

feature_names for a forest fitted from a DataFrame, else Z0, Z1, … by column of Z.

Type:

The name of each feature

ff(x: int | float | ArrayLike, Z: ArrayLike | NDArray) → NDArray

Failure (CDF) function averaged over the trees, as for sf().

classmethod fit(x: ArrayLike | None = None, Z: ArrayLike | NDArray | None = None, c: ArrayLike | None = None, n: ArrayLike | None = None, t: ArrayLike | None = None, xl: ArrayLike | None = None, xr: ArrayLike | None = None, tl: ArrayLike | None = None, tr: ArrayLike | None = None, n_trees: int = 100, max_depth: int | float = inf, min_leaf_samples: int = 5, min_leaf_failures: int = 2, n_features_split: int | float | str = 'sqrt', bootstrap: bool = True, kind: str = 'weibull', selection: str = 'greedy', alpha_split: float = 0.05, min_split_gain: float | str = 0.0, random_state: Any = None, n_jobs: int = 1) → RandomSurvivalForest

Fit a random survival forest.

Parameters:
  • x (array_like, optional) – Event times ([left, right] rows for interval-censored observations).

  • Z (array_like) – Covariate (feature) matrix, one row per observation. Required. Rows with a missing (NaN) or infinite covariate are dropped, with a warning giving the count.

  • c (array_like, optional) – Censoring flags: 0 observed, 1 right, -1 left, 2 interval censored. Defaults to all observed.

  • n (array_like, optional) – Counts. Defaults to 1.

  • t (array_like, optional) – (N, 2) truncation bounds.

  • xl (array_like, optional) – Interval bounds, instead of 2-D x.

  • xr (array_like, optional) – Interval bounds, instead of 2-D x.

  • tl (array_like, optional) – Left and right truncation, instead of t.

  • tr (array_like, optional) – Left and right truncation, instead of t.

  • max_depth (int, optional) – Maximum depth of a tree. Defaults to unlimited.

  • min_leaf_samples (int, optional) – A split is only made if each child keeps at least this many observations. Defaults to 5.

  • min_leaf_failures (int, optional) – … and at least this many failures (rows that are not right censored, each counted n times). Defaults to 2.

  • n_features_split (int, float or str, optional) – The number of features considered at each split: an int, a fraction of the features (float), "sqrt" (the default), "log2" or "all".

  • n_trees (int, optional) – The number of trees. Defaults to 100.

  • bootstrap (bool, optional) – Fit each tree to a bootstrap resample of the data (the default); otherwise every tree sees all of it.

  • kind (str, optional) – The tree type, "weibull" (the default), "exponential" or "non-parametric"; see SurvivalTree.

  • selection (str, optional) – How each node chooses its feature: "greedy" (the default) or "ctree" (conditional inference, which also stops a tree where the data show no effect); see SurvivalTree.

  • alpha_split (float, optional) – With selection="ctree", a node splits only if the Bonferroni-adjusted p-value of its chosen feature is below alpha_split. Defaults to 0.05.

  • min_split_gain (float, "aic" or "bic", optional) – The least gain in log-likelihood a split of a "weibull" or "exponential" tree must make: a node splits only if its best cut raises the maximised log-likelihood of its working model by more than this (the two children’s against the node’s). "aic" is the kind’s degrees of freedom k (1 for "exponential", 2 for "weibull"): the split must lower Akaike’s criterion. "bic" is k log(d) / 2, with d the node’s failures (rows not right censored, counted n times; its units if it has none), as every BIC in SurPyval counts them: the split must lower the Bayesian criterion. Defaults to 0: any gain, as a forest of deep trees wants. ("aic" is the recommended setting for a single SurvivalTree.) Not used by "non-parametric" trees, whose splits are not likelihoods; stop those with selection="ctree".

  • random_state (None, int or numpy.random.Generator, optional) – Seeds the bootstrap resamples and the features drawn for each split. None (the default) draws from NumPy’s global random state, so np.random.seed reproduces the forest; a seed or Generator gives the forest a stream of its own (and each tree a child stream of it) and leaves the global one alone. The forest is the same whatever n_jobs is, except with None (see n_jobs).

  • n_jobs (int, optional) – The number of worker processes the trees are grown in, as in joblib and scikit-learn: 1 (the default) grows them one after another in this process, -1 uses every core. The trees and their predictions do not depend on it given a seed. With random_state=None and n_jobs other than 1, the trees’ feature draws come from streams seeded by one draw from NumPy’s global state (worker processes do not share it), so np.random.seed still reproduces the forest for that n_jobs, but not the n_jobs=1 forest.

Returns:

The fitted forest.

Return type:

RandomSurvivalForest

Examples

>>> import numpy as np
>>> from surpyval.beta.ml import RandomSurvivalForest
>>> rng = np.random.default_rng(0)
>>> Z = rng.uniform(0, 1, (200, 2))
>>> x = rng.weibull(2.0, 200) * np.where(Z[:, 0] > 0.5, 5.0, 10.0)
>>> c = (x > 12).astype(int)
>>> x = np.minimum(x, 12)
>>> np.random.seed(0)
>>> forest = RandomSurvivalForest.fit(
...     x, Z, c=c, n_trees=5, max_depth=1, kind="exponential"
... )
>>> forest.sf(5, [[0.2, 0.5], [0.8, 0.5]]).round(3)
array([0.561, 0.396])

A seed of its own reproduces the forest without touching NumPy’s global state:

>>> a = RandomSurvivalForest.fit(
...     x, Z, c=c, n_trees=5, max_depth=1, kind="exponential",
...     random_state=1,
... )
>>> b = RandomSurvivalForest.fit(
...     x, Z, c=c, n_trees=5, max_depth=1, kind="exponential",
...     random_state=1,
... )
>>> bool(np.array_equal(a.sf(5, Z[:3]), b.sf(5, Z[:3])))
True
classmethod fit_from_df(df: pd.DataFrame, x_col: str | None = None, Z_cols: str | list[str] | None = None, c_col: str | None = None, n_col: str | None = None, tl_col: str | None = None, tr_col: str | None = None, xl_col: str | None = None, xr_col: str | None = None, formula: str | None = None, **fit_options: Any) → Any

Fit to data held in the columns of a pandas.DataFrame.

The column names are passed in place of the arrays fit() takes, with the names of every regression fit_from_df; every other fit() option is passed to it unchanged. The columns are handed to fit() as they are, so the result, and the treatment of a missing covariate, is that of fit() on the same arrays. The fitted model keeps the covariate names as feature_names (and the formula): its split descriptions and feature importances are named by them, and it predicts from a DataFrame by them.

Parameters:
  • df (pandas.DataFrame) – The data.

  • x_col (str, optional) – Column of observed times. Required unless xl_col and xr_col are given.

  • Z_cols (str or list of str, optional) – Column(s) of the covariates, in the order fit() reads them (the columns of Z). Exactly one of Z_cols and formula must be given.

  • c_col (str, optional) – Column of censoring flags.

  • n_col (str, optional) – Column of counts.

  • tl_col (str, optional) – Columns of left / right truncation.

  • tr_col (str, optional) – Columns of left / right truncation.

  • xl_col (str, optional) – Columns of the left and right ends of each observation’s interval, in place of x_col (with c_col giving the censoring of each row, as for fit() with xl/xr).

  • xr_col (str, optional) – Columns of the left and right ends of each observation’s interval, in place of x_col (with c_col giving the censoring of each row, as for fit() with xl/xr).

  • formula (str, optional) – A formulaic formula for the covariates, e.g. "age + C(sex)", as for the regression models’ fit_from_df: categoricals get reference-level columns and the intercept is dropped. Its columns are the features.

  • **fit_options – Every other option of fit().

Returns:

The model fit() returns, with feature_names (and formula) set.

Return type:

model

Raises:

ValueError – If df is not a DataFrame, a name is not one of its columns, or the times or covariates are not given exactly once.

Examples

>>> import numpy as np
>>> import pandas as pd
>>> from surpyval.beta.ml import SurvivalTree
>>> rng = np.random.default_rng(0)
>>> df = pd.DataFrame({"z": rng.uniform(0, 1, 60)})
>>> df["x"] = rng.weibull(2, 60) * np.where(df["z"] > 0.5, 5, 20)
>>> tree = SurvivalTree.fit_from_df(
...     df, x_col="x", Z_cols="z", random_state=0
... )
>>> tree.sf(10.0, [[0.2], [0.8]]).round(3)
array([0.486, 0.029])
>>> tree.feature_names
['z']
>>> tree.sf(10.0, pd.DataFrame({"z": [0.2, 0.8]})).round(3)
array([0.486, 0.029])
classmethod from_dict(model_dict: dict) → RandomSurvivalForest

Reconstruct a fitted forest from a to_dict() dictionary.

classmethod from_json(fp: str | PathLike) → Any

Load a model from a JSON file written by to_json(), or from the JSON text it returned (a string starting with {).

hf(x: int | float | ArrayLike, Z: ArrayLike | NDArray) → NDArray

Hazard rate averaged over the trees, as for sf().

mortality(x: int | float | ArrayLike, Z: ArrayLike | NDArray) → ArrayLike

The ensemble mortality of each covariate vector: its cumulative hazard summed over the times x (the risk score used by score()).

oob_log_likelihood() → float

The mean out-of-bag log-likelihood per observation.

Each row of the training data is scored by the ensemble of the trees whose bootstrap sample left it out, so by trees that never saw it: its contribution is the log of

  • the density \(f(x)\) if it was observed at \(x\),

  • \(S(x)\) if right censored at \(x\),

  • \(1 - S(x)\) if left censored at \(x\),

  • \(S(x_l) - S(x_r)\) if interval censored in \((x_l, x_r]\),

divided by \(S(t_l) - S(t_r)\) if it is truncated to \((t_l, t_r]\), where \(S\) and \(f\) are the averages of the out-of-bag trees’ leaf survival functions and densities. A censored row’s interval is first cut to its truncation window, as in the fitters’ likelihoods. The result is the count-weighted mean over the rows, so higher is better and it estimates the expected log-likelihood of a new observation; it works for every censoring type and truncation, where the concordance of score() needs orderable event times.

A non-parametric leaf is a step function, which puts no probability at an out-of-bag event time unless the tree saw a tied one, so for this score it is read as a continuous distribution: its survival curve is joined linearly between the points where it drops, from 1 at time 0 (or at the smallest time, if it is negative), and continued past its last drop with the constant hazard it averaged up to there (Brown, Hollander and Korwar’s exponential tail). Its density is then per unit of time, on the same scale as a parametric leaf’s, so forests of different kind can be compared by this score.

A row that is in the bootstrap sample of every tree has no out-of-bag prediction; it is left out of the mean, with one warning giving the count (every row, and a NaN result, with bootstrap=False). A row the ensemble gives zero probability makes the mean -inf, with one warning giving the count of such rows: with few trees a row can land only in leaves that put no density at its time (a steep Weibull leaf grown on bunched failures, or a step leaf); grow more trees, or use kind="exponential", whose leaves give every time a density.

Returns:

The mean out-of-bag log-likelihood per observation.

Return type:

float

Raises:

ValueError – On a forest restored with from_dict, which keeps neither the training data nor the bootstrap samples.

Examples

>>> import numpy as np
>>> from surpyval.beta.ml import RandomSurvivalForest
>>> rng = np.random.default_rng(0)
>>> Z = rng.uniform(0, 1, (200, 2))
>>> x = rng.weibull(2.0, 200) * np.where(Z[:, 0] > 0.5, 5.0, 10.0)
>>> c = (x > 12).astype(int)
>>> x = np.minimum(x, 12)
>>> np.random.seed(0)
>>> forest = RandomSurvivalForest.fit(
...     x, Z, c=c, n_trees=20, max_depth=1, kind="exponential"
... )
>>> round(forest.oob_log_likelihood(), 3)
-2.518
score(x: ArrayLike, Z: ArrayLike | NDArray, c: ArrayLike, tie_tol: float = 1e-08, ties: str = 'therneau') → float

Harrell’s concordance index of the forest’s mortality scores.

The index is surpyval.metrics.concordance_index(), with its tie conventions: by default (ties="therneau", as R’s survival::concordance, lifelines and every model’s concordance) two events at the same time are not a usable pair; ties="harrell" counts them, as this method did before v0.22.

A missing (NaN) covariate or time leaves a subject’s score, and so the index, undefined: the index is NaN, not a number computed by comparing the NaN score as though it were one.

sf(x: int | float | ArrayLike, Z: ArrayLike | NDArray, ensemble_method: str = 'sf') → NDArray

Returns the ensemble survival function

Parameters:
  • x (int | float | ArrayLike) – Times, the same for every covariate vector.

  • Z (ArrayLike | NDArray) – One covariate vector (1-D), or a matrix with one covariate vector per row (2-D).

  • ensemble_method (str, optional) – Determines whether to average across terminal nodes the terminal node survival functions or cumulative hazard functions. For these respectively, ensemble_method must be “sf” or “Hf”. Defaults to “sf”.

Returns:

For a 1-D Z, the survival function at x, shaped like x (a scalar for a scalar x). For a 2-D Z, a grid of shape (n_rows,) + x.shape whose row i is the survival function for Z[i] (every row at every time). A covariate vector with a missing (NaN) value gives NaN, and leaves the other rows unaffected.

Return type:

NDArray

to_dict() → dict

Serialise the fitted forest to a plain, JSON/BSON-safe dictionary: the ensemble settings and every fitted tree. The training data is not persisted – a restored forest is a predictor, not re-fittable.

to_json(fp: str | PathLike | None = None, with_data: bool = False) → str | None

Write to_dict() to fp as strict JSON, or return it.

Parameters:
  • fp (str or os.PathLike, optional) – The file to write. Without it the JSON is returned as a string (as pandas.DataFrame.to_json does), which from_json also reads.

  • with_data (bool, optional) – Write to_dict(with_data=True), which also stores the fitted data, for the models whose to_dict takes with_data (the univariate Parametric and NonParametric); a TypeError for any other model. Defaults to False.

Survival Tree

class surpyval.beta.ml.forest.tree.SurvivalTree(data: SurpyvalData, Z: NDArray, max_depth: int | float = inf, min_leaf_samples: int = 5, min_leaf_failures: int = 2, n_features_split: int | float | str = 'sqrt', kind: str = 'weibull', selection: str = 'greedy', alpha_split: float = 0.05, random_state: Any = None, feature_names: list[str] | None = None, min_split_gain: float | str = 0.0)

Bases: RegressionDataFrameMixin, SerialisableMixin

A Survival Tree, for use in RandomSurvivalForest.

The Tree is built on initialisation. Supports the full SurPyval data model: observed, left-, right- and interval-censored observations with optional left and/or right truncation.

The tree’s kind couples the split criterion with the matching leaf model, so every split greedily improves the model the tree predicts with:

  • "weibull" (default): full-likelihood Weibull deviance split (a 2-d.f. likelihood-ratio gain, with power against scale and shape differences) with Weibull MLE leaves. Supports the full data model.

  • "exponential": exponential deviance split (Davis & Anderson, 1989; 1-d.f., splits on rate) with Exponential MLE leaves. Supports the full data model.

  • "non-parametric": for observed / right-censored data (optionally left-truncated), the risk-set log-rank split with Nelson-Aalen leaves. For data with left or interval censoring or right truncation, the Turnbull-score split – the standardised sum of each child’s log-rank scores under the node’s pooled Turnbull estimate (Finkelstein, 1986), which reduces to the log-rank scores on right-censored data; a truncated row’s score is that of its truncation-conditioned likelihood – with Turnbull leaves. Supports the full data model.

selection decides how a node chooses the feature it splits on:

  • "greedy" (default): the best cut of the kind’s criterion over every feature drawn for the split. A feature with many distinct values offers more cuts, so it is favoured even when it carries no information, and a node always splits if some cut is allowed.

  • "ctree": conditional inference (Hothorn, Hornik and Zeileis, 2006). Each feature is tested for association with the scores of the kind’s split statistic (the log-rank scores for "non-parametric"; the working model’s score contributions for "exponential" and "weibull"), by its maximally selected statistic over its cuts, whose p-value allows for the number of cuts. The feature with the smallest p-value is chosen, and the node splits only if that p-value, Bonferroni-adjusted for the number of features tested, is below alpha_split; its cut is then chosen by the kind’s criterion. This removes the preference for features with many values and stops the tree where the data show no effect. See conditional_inference.

Hf(x: int | float | ArrayLike, Z: ArrayLike | NDArray) → NDArray

Cumulative hazard at x of the leaf model each covariate vector falls in, as for sf().

apply_model_function(function_name: str, x: int | float | ArrayLike, Z: ArrayLike | NDArray) → NDArray

Evaluate function_name ("sf", "ff", "df", "hf" or "Hf") of the leaf model that each covariate vector falls in.

Parameters:
  • function_name (str) – The name of the leaf model’s function to evaluate.

  • x (int, float or array_like) – Times, the same for every covariate vector.

  • Z (array_like) – One covariate vector (1-D), or a matrix with one covariate vector per row (2-D).

Returns:

For a 1-D Z, the values at x, shaped like x (a scalar for a scalar x). For a 2-D Z, a grid of shape (n_rows,) + x.shape whose row i is the values for Z[i] – every row at every time, the one documented exception to pairing rows with times – as for RandomSurvivalForest. A covariate vector with a missing (NaN) value gives NaN, and leaves the other rows unaffected.

Return type:

ndarray

describe() → str

The tree as text: one line per split, name <= value for the left branch and name >  value for the right one (with the adjusted p-value of a selection="ctree" split), each branch’s subtree indented under it, and each leaf’s model.

df(x: int | float | ArrayLike, Z: ArrayLike | NDArray) → NDArray

Density at x of the leaf model each covariate vector falls in, as for sf().

property feature_labels: list[str]

feature_names for a tree fitted from a DataFrame (fit_from_df, or fit with a DataFrame Z), else Z0, Z1, … by column of Z. Split descriptions and the printout use them.

Type:

The name of each feature

ff(x: int | float | ArrayLike, Z: ArrayLike | NDArray) → NDArray

Failure (CDF) function at x of the leaf model each covariate vector falls in, as for sf().

classmethod fit(x: ArrayLike | None = None, Z: ArrayLike | NDArray | None = None, c: ArrayLike | None = None, n: ArrayLike | None = None, t: ArrayLike | None = None, xl: ArrayLike | None = None, xr: ArrayLike | None = None, tl: ArrayLike | None = None, tr: ArrayLike | None = None, max_depth: int | float = inf, min_leaf_samples: int = 5, min_leaf_failures: int = 2, n_features_split: int | float | str = 'sqrt', kind: str = 'weibull', selection: str = 'greedy', alpha_split: float = 0.05, min_split_gain: float | str = 0.0, random_state: Any = None) → SurvivalTree

Fit a survival tree from data in the full xcnt(+truncation) data model.

x/c/n/t follow the standard SurPyval conventions (c in {-1, 0, 1, 2}; interval-censored entries of x are [left, right] pairs). Interval bounds can alternatively be given as xl/xr, and truncation as tl/tr instead of the two-column t. kind selects the tree type (see the class docstring).

Parameters:
  • x (array_like, optional) – Event times ([left, right] rows for interval-censored observations).

  • Z (array_like) – Covariate (feature) matrix, one row per observation. Required. Rows with a missing (NaN) or infinite covariate are dropped, with a warning giving the count.

  • c (array_like, optional) – Censoring flags: 0 observed, 1 right, -1 left, 2 interval censored. Defaults to all observed.

  • n (array_like, optional) – Counts. Defaults to 1.

  • t (array_like, optional) – (N, 2) truncation bounds.

  • xl (array_like, optional) – Interval bounds, instead of 2-D x.

  • xr (array_like, optional) – Interval bounds, instead of 2-D x.

  • tl (array_like, optional) – Left and right truncation, instead of t.

  • tr (array_like, optional) – Left and right truncation, instead of t.

  • max_depth (int, optional) – Maximum depth of a tree. Defaults to unlimited.

  • min_leaf_samples (int, optional) – A split is only made if each child keeps at least this many observations. Defaults to 5.

  • min_leaf_failures (int, optional) – … and at least this many failures (rows that are not right censored, each counted n times). Defaults to 2.

  • n_features_split (int, float or str, optional) – The number of features considered at each split: an int, a fraction of the features (float), "sqrt" (the default), "log2" or "all".

  • kind (str, optional) – "weibull" (the default), "exponential" or "non-parametric"; see the class docstring.

  • selection (str, optional) – How a node chooses its feature: "greedy" (the default), the best cut over every feature, or "ctree", conditional inference; see the class docstring.

  • alpha_split (float, optional) – With selection="ctree", a node splits only if the Bonferroni-adjusted p-value of its chosen feature is below alpha_split, the size of the test of no association. Defaults to 0.05. Ignored by "greedy".

  • min_split_gain (float, "aic" or "bic", optional) – The least gain in log-likelihood a split of a "weibull" or "exponential" tree must make: a node splits only if its best cut raises the maximised log-likelihood of its working model by more than this (the two children’s against the node’s). "aic" is the kind’s degrees of freedom k (1 for "exponential", 2 for "weibull"): the split must lower Akaike’s criterion. "bic" is k log(d) / 2, with d the node’s failures (rows not right censored, counted n times; its units if it has none), as every BIC in SurPyval counts them: the split must lower the Bayesian criterion. Defaults to 0: any gain, as a forest of deep trees wants. "aic" is the recommended setting for a tree used on its own, which otherwise splits on noise until min_leaf_samples or min_leaf_failures stops it. Not used by "non-parametric" trees, whose splits are not likelihoods; stop those with selection="ctree".

  • random_state (None, int or numpy.random.Generator, optional) – Seeds the features drawn for each split (when n_features_split is less than the number of features). None (the default) draws from NumPy’s global random state, so np.random.seed reproduces the tree; a seed or Generator gives a stream of its own and leaves the global one alone.

Returns:

The fitted tree. Its sf(x, Z) (and ff, df, hf, Hf) evaluate the model of the leaf that a covariate vector Z falls in; a matrix Z gives one row per covariate vector and one column per time. A covariate vector with a missing (NaN) value gives NaN.

Return type:

SurvivalTree

Examples

Life halves when the first feature exceeds 0.5; a single split finds it:

>>> import numpy as np
>>> from surpyval.beta.ml import SurvivalTree
>>> rng = np.random.default_rng(0)
>>> Z = rng.uniform(0, 1, (200, 2))
>>> x = rng.weibull(2.0, 200) * np.where(Z[:, 0] > 0.5, 5.0, 10.0)
>>> c = (x > 12).astype(int)
>>> x = np.minimum(x, 12)
>>> tree = SurvivalTree.fit(x, Z, c=c, max_depth=1, n_features_split=2)
>>> tree.sf(5, [0.2, 0.5]).round(4), tree.sf(5, [0.8, 0.5]).round(4)
(np.float64(0.8831), np.float64(0.3168))

A matrix routes each row to its own leaf:

>>> tree.sf([2, 5], [[0.2, 0.5], [0.8, 0.5]]).round(4)
array([[0.9897, 0.8831],
       [0.8062, 0.3168]])

With conditional-inference selection, a tree grown on the same data without the effect does not split at all, where greedy search always does:

>>> x0 = rng.weibull(2.0, 200) * 10.0
>>> ctree = SurvivalTree.fit(
...     x0, Z, kind="non-parametric", n_features_split="all",
...     selection="ctree",
... )
>>> type(ctree._root).__name__
'TerminalNode'
>>> ctree = SurvivalTree.fit(
...     x, Z, c=c, kind="non-parametric", n_features_split="all",
...     selection="ctree", max_depth=1,
... )
>>> root = ctree._root
>>> int(root.split_feature_index), bool(root.p_value < 1e-10)
(0, True)
classmethod fit_from_df(df: pd.DataFrame, x_col: str | None = None, Z_cols: str | list[str] | None = None, c_col: str | None = None, n_col: str | None = None, tl_col: str | None = None, tr_col: str | None = None, xl_col: str | None = None, xr_col: str | None = None, formula: str | None = None, **fit_options: Any) → Any

Fit to data held in the columns of a pandas.DataFrame.

The column names are passed in place of the arrays fit() takes, with the names of every regression fit_from_df; every other fit() option is passed to it unchanged. The columns are handed to fit() as they are, so the result, and the treatment of a missing covariate, is that of fit() on the same arrays. The fitted model keeps the covariate names as feature_names (and the formula): its split descriptions and feature importances are named by them, and it predicts from a DataFrame by them.

Parameters:
  • df (pandas.DataFrame) – The data.

  • x_col (str, optional) – Column of observed times. Required unless xl_col and xr_col are given.

  • Z_cols (str or list of str, optional) – Column(s) of the covariates, in the order fit() reads them (the columns of Z). Exactly one of Z_cols and formula must be given.

  • c_col (str, optional) – Column of censoring flags.

  • n_col (str, optional) – Column of counts.

  • tl_col (str, optional) – Columns of left / right truncation.

  • tr_col (str, optional) – Columns of left / right truncation.

  • xl_col (str, optional) – Columns of the left and right ends of each observation’s interval, in place of x_col (with c_col giving the censoring of each row, as for fit() with xl/xr).

  • xr_col (str, optional) – Columns of the left and right ends of each observation’s interval, in place of x_col (with c_col giving the censoring of each row, as for fit() with xl/xr).

  • formula (str, optional) – A formulaic formula for the covariates, e.g. "age + C(sex)", as for the regression models’ fit_from_df: categoricals get reference-level columns and the intercept is dropped. Its columns are the features.

  • **fit_options – Every other option of fit().

Returns:

The model fit() returns, with feature_names (and formula) set.

Return type:

model

Raises:

ValueError – If df is not a DataFrame, a name is not one of its columns, or the times or covariates are not given exactly once.

Examples

>>> import numpy as np
>>> import pandas as pd
>>> from surpyval.beta.ml import SurvivalTree
>>> rng = np.random.default_rng(0)
>>> df = pd.DataFrame({"z": rng.uniform(0, 1, 60)})
>>> df["x"] = rng.weibull(2, 60) * np.where(df["z"] > 0.5, 5, 20)
>>> tree = SurvivalTree.fit_from_df(
...     df, x_col="x", Z_cols="z", random_state=0
... )
>>> tree.sf(10.0, [[0.2], [0.8]]).round(3)
array([0.486, 0.029])
>>> tree.feature_names
['z']
>>> tree.sf(10.0, pd.DataFrame({"z": [0.2, 0.8]})).round(3)
array([0.486, 0.029])
classmethod from_dict(model_dict: dict) → SurvivalTree

Reconstruct a fitted tree from a to_dict() dictionary.

classmethod from_json(fp: str | PathLike) → Any

Load a model from a JSON file written by to_json(), or from the JSON text it returned (a string starting with {).

hf(x: int | float | ArrayLike, Z: ArrayLike | NDArray) → NDArray

Hazard rate at x of the leaf model each covariate vector falls in, as for sf().

sf(x: int | float | ArrayLike, Z: ArrayLike | NDArray) → NDArray

Survival function at x of the leaf model each covariate vector falls in; Z and the result are as for apply_model_function() (a 2-D Z gives one row per covariate vector).

to_dict() → dict

Serialise the fitted tree to a plain, JSON/BSON-safe dictionary.

Only what prediction needs is stored – the tree kind, the resolved n_features_split, the selection and alpha_split it was grown with, and the recursive node structure with its fitted leaf models (and, for selection="ctree", each split’s p-value). The training data and covariate matrix are not persisted: a restored tree is a predictor, not a re-fittable object.

to_json(fp: str | PathLike | None = None, with_data: bool = False) → str | None

Write to_dict() to fp as strict JSON, or return it.

Parameters:
  • fp (str or os.PathLike, optional) – The file to write. Without it the JSON is returned as a string (as pandas.DataFrame.to_json does), which from_json also reads.

  • with_data (bool, optional) – Write to_dict(with_data=True), which also stores the fitted data, for the models whose to_dict takes with_data (the univariate Parametric and NonParametric); a TypeError for any other model. Defaults to False.

Tree Nodes

The nodes a fitted tree is built from. Users do not normally construct these directly; they are documented because a serialised tree is a nested structure of them, and because TerminalNode.model is how a leaf’s fitted distribution is reached.

class surpyval.beta.ml.forest.node.Node

Bases: ABC

The common methods between IntermediateNode and TerminalNode.

abstractmethod apply_model_function(function_name: str, x: int | float | ArrayLike, Z: NDArray) → NDArray

Evaluate function_name ("sf", "Hf", …) of the leaf model(s) reached by Z at the times x.

A 1-D Z is one covariate vector and returns that leaf’s values at x. A 2-D Z holds one covariate vector per row; each row is routed on its own and the result is an (n_rows, x.size) grid, row i being the values for Z[i].

class surpyval.beta.ml.forest.node.IntermediateNode(data: SurpyvalData, Z: NDArray, curr_depth: int, max_depth: int | float, min_leaf_samples: int, min_leaf_failures: int, n_features_split: int, split_feature_index: int, split_feature_value: float, feature_indices_in: NDArray, kind: str = 'weibull', rng: Any = None, selection: str = 'greedy', alpha_split: float = 0.05, p_value: float | None = None, min_split_gain: float | str = 0.0)

Bases: Node

A split in a survival tree: observations whose feature split_feature_index is at most split_feature_value go to left_child, the rest to right_child. Building one grows the subtree below it. In a tree grown with selection="ctree", p_value is the Bonferroni-adjusted p-value that chose the split’s feature (None otherwise).

apply_model_function(function_name: str, x: int | float | ArrayLike, Z: NDArray) → NDArray

Evaluate function_name ("sf", "Hf", …) of the leaf model(s) reached by Z at the times x.

A 1-D Z is one covariate vector and returns that leaf’s values at x. A 2-D Z holds one covariate vector per row; each row is routed on its own and the result is an (n_rows, x.size) grid, row i being the values for Z[i].

describe(feature_names: list[str] | None = None, right: bool = False) → str

The split rule as text, "temp <= 42" ("temp >  42" for the right branch with right=True), naming the feature by feature_names (Z3 for column 3 without them).

to_dict() → dict

Serialise the split rule and both child subtrees. The training data is deliberately not stored – a restored tree is a predictor, rebuilt from its structure and its leaf models, not re-fitted.

class surpyval.beta.ml.forest.node.TerminalNode(data: SurpyvalData, kind: str = 'weibull')

Bases: Node

A leaf of a survival tree. It holds the observations that reach it and fits, on first use, the leaf model given by the tree’s kind (model): a Weibull or Exponential fit, or for a non-parametric tree a Nelson-Aalen estimate (a Turnbull estimate if the leaf holds left- or interval-censored or right-truncated rows); NeverOccurs for a parametric leaf with no failures. On observed and right-censored data a parametric leaf is the maximum found as the split search finds a child’s, built from its parameters (so it has no cb() of its own).

apply_model_function(function_name: str, x: int | float | ArrayLike, Z: NDArray) → NDArray

Evaluate function_name ("sf", "Hf", …) of the leaf model(s) reached by Z at the times x.

A 1-D Z is one covariate vector and returns that leaf’s values at x. A 2-D Z holds one covariate vector per row; each row is routed on its own and the result is an (n_rows, x.size) grid, row i being the values for Z[i].

describe() → str

The leaf as text: its model (with the parameters of a parametric one) and, on a fitted tree, the number of units that reached it.

property model: Any

The leaf’s fitted model, fitted when first used.

to_dict() → dict

Serialise the leaf as its fitted model rather than its data, so the restored leaf predicts without re-fitting. NeverOccurs (the empty / all-censored leaf) is a parameterless class, stored as a string sentinel; every other leaf is a Parametric or NonParametric model with its own to_dict.

property units: float | None

The number of units (n-weighted rows) that reached the leaf when the tree was grown (kept by a restored tree; None for one saved before it was stored).