Regression Modelling with SurPyval
The time until an event — a failure, death, recovery — will almost always depend on external factors. A bearing may last longer in a cool, clean environment than in a hot, dirty one. A patient’s survival may depend on age, dosage, and comorbidities. The question regression modelling answers is: how much do these factors matter, and in what direction?
For the concepts and mathematics behind the regression families used here — the proportional-hazards, accelerated-failure-time, accelerated-life, proportional-odds and additive-hazards models, the Cox partial likelihood and its tie corrections, the Cox diagnostics, robust, stratified and frailty fits, time-varying covariates, survival trees, and prediction-validation metrics — see the Regression Analysis page. This page is the practical companion: every example runs, and each section links back to the theory it relies on. The full API reference is under Regression Modelling.
Regression survival modelling is fundamentally about capturing the relationship between covariates \(Z\) and the survival distribution. Unlike ordinary regression, we must handle censored observations — items that had not yet failed when we stopped watching them — and we want our model to remain valid as a probability distribution (survival functions must start at 1 and decay to 0).
For the rest of this page we assume the following imports:
import surpyval as surv
import numpy as np
from matplotlib import pyplot as plt
How this page is organised. It is long, and it builds up in four parts. First, how the families differ and the conventions every fitter shares (the data arguments and how predictions pair times with covariates). Second, the semi-parametric models, which leave the baseline to the data: the Cox model with its whole toolkit — tied times, delayed entry, DataFrames and formulas, time-varying covariates, the proportional-hazards check, robust errors and strata — then its accelerated-time counterpart Buckley-James and the additive-hazards model. Third, the parametric families — proportional hazards, accelerated failure time, proportional odds, their confidence bounds, accelerated life testing, time-varying covariates across families and shared frailty. Fourth, choosing between fitted models, validating their predictions, and saving them. If you already know the question you are asking, the table below points straight to the section that answers it.
If you want to … |
use |
see |
|---|---|---|
estimate hazard ratios without assuming a lifetime distribution |
|
|
predict the whole lifetime distribution, extrapolate beyond the data, or use left- or interval-censored data |
|
|
say “this factor costs a fraction of the life” (a time ratio) |
|
Accelerated Failure Time (AFT), Semi-Parametric — Buckley-James (AFT) |
carry accelerated-test results to use conditions through a physical stress-life law |
|
|
model an effect that fades as time goes on |
|
|
report an excess risk (extra failures per unit time) |
|
|
use covariates that change during follow-up, or forecast along a planned covariate path |
|
Time-Varying Covariates, Time-varying covariates across families |
check that a hazard ratio really is constant |
|
|
get honest standard errors for grouped or repeated data |
|
|
model, and predict for, the variation between groups |
|
|
remove a nuisance factor that breaks proportional hazards |
|
|
find structure you cannot specify (thresholds, interactions) |
|
|
compare, validate or store fitted models |
|
Model Selection, Validating a survival predictor, Saving and loading a fitted model |
Choosing a regression model family
There are four fundamentally different ways a covariate can affect a survival distribution through a simple link, plus the physics-driven accelerated life model. Each gives rise to a distinct model family:
Proportional Hazards (PH) — the covariate multiplies the rate of dying:
If \(\phi(Z) = 2\), an individual with that covariate value fails at twice the rate at every instant in time. The shape of the hazard curve is unchanged; only its level shifts. This is the most common choice in medical research.
Accelerated Failure Time (AFT) — the covariate stretches or compresses the time axis:
If \(\phi(Z) = 2\), an individual “ages” at twice the normal rate — reaching at age 10 the same cumulative risk that a baseline individual has at age 20. The entire survival curve shifts left or right on the (log) time axis. This is often a more natural framing in engineering and materials science.
Proportional Odds (PO) — the covariate scales the odds of surviving:
The PO formulation is natural when you think about the problem in terms of odds ratios rather than hazard ratios. A key practical difference from PH: the covariate effect attenuates over time. Early in life, hazard ratios and odds ratios behave similarly; at long follow-up times the PO effect fades as everyone converges toward failure regardless of their covariates. When you believe the PH assumption (“constant hazard ratio for all time”) is too strong, PO is often a better default.
Additive Hazards (AH) — the covariate adds to the hazard instead of multiplying it:
Here \(\beta_j\) is a risk difference — the change in the absolute hazard per unit of covariate, constant over time — rather than a hazard ratio. This is the natural scale when you care about how many extra failures per unit time a factor causes: excess risk in epidemiology, or reliability settings where hazards from independent mechanisms genuinely add. Because the effect is additive rather than exponential it is not constrained to be positive, which is both its interpretive appeal and its main caveat (discussed below).
Accelerated Life (AL) — the covariate substitutes the distribution’s life parameter:
This is the standard approach in accelerated life testing (ALT) — reliability testing under elevated stress (high temperature, voltage, humidity) to extract failure data quickly, then extrapolating back to use conditions. The stress relationship \(\phi(Z)\) is chosen from domain knowledge: Arrhenius for thermally-activated failure, Eyring for quantum-mechanical processes, Power Law for voltage or mechanical loading.
For PH, AFT, and PO the covariate function is always the log-linear form,
because the exponential guarantees \(\phi > 0\) for any covariate value and any β — no parameter constraints needed. For PH and AFT a positive β makes failure faster and a negative one slower. Proportional odds is the exception: its \(\phi\) multiplies the odds of survival, so a positive β makes failure slower. Expect a PO fit to report coefficients of the opposite sign to a PH fit on the same data.
SurPyval supports all of these families, each available as pre-built instances for every standard distribution and as factory functions for custom combinations, together with semi-parametric versions that leave the baseline unspecified, a frailty version for grouped data, and tree-based predictors:
Family |
Effect of covariates |
Ready-to-use examples |
|---|---|---|
Proportional Hazards (PH) |
Multiplies the hazard rate \(h(x|Z) = h_0(x)\,\phi(Z)\) |
|
Accelerated Failure Time (AFT) |
Scales the time axis \(H(x|Z) = H_0(\phi(Z)\,x)\) |
|
Proportional Odds (PO) |
Scales the survival odds \(O(x|Z) = O_0(x)\,\phi(Z)\) |
|
Additive Hazards (AH) |
Adds to the hazard rate \(h(x|Z) = h_0(x) + \beta'Z\) |
|
Accelerated Life (AL) |
Substitutes the life parameter with a physics-motivated function |
|
Shared frailty PH |
PH with a random multiplier shared within a group |
|
Survival trees and forests (beta) |
No link: recursive splits on the covariates |
|
Data, covariates and predictions
Every regression fitter takes the observed times x and a covariate matrix
Z with one row per observation and one column per covariate (a
one-dimensional Z is read as a single covariate, one value per row) —
plus surpyval’s usual optional arrays: the censoring flag c (0 observed,
1 right, -1 left, 2 interval censored) and counts n. What each
fitter accepts:
the parametric families (PH, AFT, PO, AH and AL) accept every censoring type, and truncation through
t(a two-column[tl, tr]array);CoxPHandProportionalOddstake observed and right-censored data, with left truncation through a 1-Dtl, and refuse left- or interval-censored rows (their likelihoods have no term for them);the Lin-Ying (
AdditiveHazards), Buckley-James and frailty fitters take observed and right-censored data only, and say so if given anything else.
A row with a missing (NaN) or infinite covariate cannot enter any of
these likelihoods, so every fitter drops it — from the times, flags, counts
and truncation too — and warns with the number of rows dropped. The same
goes for fit_from_df, with named columns or a formula. (The
time-varying-covariate fits are the exception: dropping one interval would
change a subject’s history, so they refuse a missing covariate instead.)
Predicting from a DataFrame row with a missing covariate – numeric or
categorical – gives nan for that row, in its place.
A covariate that separates the events from the survivors – a level of a
factor with no events, say – has no finite estimate: the likelihood keeps
increasing as its coefficient grows. Every fitter says so with one warning
naming the coefficient (CoxPH and FineGray as a “monotone partial
likelihood”; the parametric, additive-hazards and frailty fits as “no finite
maximum”) and returns the model where the search stopped, whose value for that
coefficient, its standard error and its bounds mean nothing. Remove or coarsen
the covariate (merge the level with another), or fit a penalised model.
Each family also has a fit_from_df that names DataFrame columns instead
(see Fitting from a DataFrame: formulas and categorical covariates).
Predictions — sf, ff, df, hf and Hf — take times and
covariates. Given one covariate row they return the curve over all the
times; given n rows and n times they pair them element-wise, one
time per row, which is what you want for scoring a data set; any other number
of rows is refused with a ValueError. For a curve per covariate row –
every time for every row, lifelines’ predict_survival_function – pass
grid=True: the result has shape (len(Z),) + x.shape, row i for
row i of Z, as the survival tree and forest return it.
A small simulated data set shows the three forms:
from surpyval import WeibullPH
rng = np.random.default_rng(42)
Z_demo = rng.binomial(1, 0.5, size=(300, 1)).astype(float)
# Weibull(10, 2) baseline; exposure multiplies the hazard by e^0.7 ~ 2
x_demo = 10 * rng.weibull(2.0, 300) * np.exp(-0.7 * Z_demo[:, 0] / 2.0)
demo = WeibullPH.fit(x=x_demo, Z=Z_demo)
print('alpha, beta, beta_0 :', demo.params.round(3))
print('one row, three times :', demo.sf([5.0, 10.0, 15.0], Z=[1.0]).round(3))
print('two rows, paired :', demo.sf([5.0, 5.0], Z=[[0.0], [1.0]]).round(3))
print('two rows, a grid :')
print(demo.sf([5.0, 10.0, 15.0], Z=[[0.0], [1.0]], grid=True).round(3))
alpha, beta, beta_0 : [10.339 2.032 0.725]
one row, three times : [0.624 0.145 0.012]
two rows, paired : [0.796 0.624]
two rows, a grid :
[[0.796 0.393 0.119]
[0.624 0.145 0.012]]
The regression models do not have a quantile function (qf). A quantile at a
given covariate value is the root of \(S(x \mid Z) = 1 - p\), which a
bracketing root-finder finds reliably because sf is monotone:
from scipy.optimize import brentq
for z in [0.0, 1.0]:
median = brentq(lambda t: demo.sf([t], Z=[z])[0] - 0.5, 1e-6, 100.0)
print(f'median life at Z = {z:g}: {median:.2f}')
median life at Z = 0: 8.63
median life at Z = 1: 6.04
With a hazard ratio of about 2 and a Weibull shape of 2, the exposed median is shorter by a factor of about \(2^{1/2}\) — exactly what the PH/AFT equivalence for the Weibull (see Accelerated Failure Time (AFT)) predicts.
Semi-Parametric — Cox Proportional Hazards
The Cox PH model is the most widely used survival regression model in any field. Its central insight is that \(\beta\) can be estimated without specifying the shape of the baseline hazard \(h_0(x)\). The baseline cancels out of a partial likelihood, leaving only the relative ordering of event times. This means you can detect and quantify covariate effects even when you have no idea what the baseline distribution looks like (the derivation is in the Regression Analysis page).
The price you pay is that the model is semi-parametric: once you have β, the baseline is estimated non-parametrically (a step function with jumps only at observed times). This means predictions can only be made within the observed time range, and extrapolation is not possible. If you need to predict far beyond your observed data, a parametric PH model is more appropriate.
In this example we use data from Krivtsov et al., testing tires to failure with seven measured characteristics. We want to know which characteristics significantly affect tire life.
from surpyval.datasets import load_tires_data
from surpyval import CoxPH
tires = load_tires_data()
x = tires['Survival']
c = tires['Censoring']
Z = tires[['Tire age', 'Wedge gauge', 'Interbelt gauge', 'EB2B', 'Peel force',
'Carbon black (%)', 'Wedge gauge×peel force']]
model = CoxPH.fit(x=x, Z=Z, c=c)
model
Semi-Parametric Regression SurPyval Model
=========================================
Type : Proportional Hazards
Kind : Cox
Parameterization : Semi-Parametric
Data : 34 units: 11 events at 8 unique times, 23 right censored
Tie method : efron
Coefficients : exp(coef) is the hazard ratio; Wald 95% intervals
coef exp(coef) se(coef) lower 95% upper 95% z p
beta_0 2.208 9.094 1.419 -0.5744 4.99 1.555 0.1199
beta_1 -10.29 3.391e-05 4.69 -19.49 -1.099 -2.194 0.02822
beta_2 -11.5 1.011e-05 4.727 -20.77 -2.238 -2.433 0.01496
beta_3 -14.74 3.967e-07 8.245 -30.9 1.42 -1.788 0.07382
beta_4 -36.46 1.468e-16 13.99 -63.88 -9.037 -2.606 0.009162
beta_5 -55.07 1.207e-24 34.2 -122.1 11.95 -1.611 0.1073
beta_6 22.42 5.465e+09 9.081 4.624 40.22 2.469 0.01354
We can immediately check which coefficients are statistically significant:
print(model.p_values)
[0.11986867 0.02822124 0.0149579 0.07381522 0.00916244 0.10728548
0.01354104]
Several covariates are not significant at the 5% level (the first, fourth and sixth). We can re-fit with only the significant ones, which also improves numerical stability — there are only 34 tires, 11 of them failures, so every extra coefficient is expensive:
Z = tires[['Wedge gauge', 'Interbelt gauge', 'Peel force',
'Wedge gauge×peel force']]
model = CoxPH.fit(x=x, Z=Z, c=c)
print(model.p_values)
model
[0.01986279 0.00978738 0.00818545 0.00839949]
Semi-Parametric Regression SurPyval Model
=========================================
Type : Proportional Hazards
Kind : Cox
Parameterization : Semi-Parametric
Data : 34 units: 11 events at 8 unique times, 23 right censored
Tie method : efron
Coefficients : exp(coef) is the hazard ratio; Wald 95% intervals
coef exp(coef) se(coef) lower 95% upper 95% z p
beta_0 -9.647 6.461e-05 4.142 -17.77 -1.528 -2.329 0.01986
beta_1 -7.446 0.0005838 2.882 -13.1 -1.797 -2.583 0.009787
beta_2 -28.71 3.386e-13 10.86 -50 -7.431 -2.644 0.008185
beta_3 19.1 1.978e+08 7.248 4.897 33.31 2.636 0.008399
The first three coefficients are negative, meaning higher gauge and peel force values reduce the hazard rate (improve life); the positive interaction term captures a counteracting combined effect.
A coefficient is a log hazard ratio per unit of its covariate. The fitted model
keeps the score and information closures of its partial likelihood
(model.jac(beta) returns the pair), so the model-based standard errors are
the square roots of the diagonal of the inverse information:
info = model.jac(model.beta)[1] # observed information
se = np.sqrt(np.diag(np.linalg.inv(info)))
for name, b, s, p in zip(Z.columns, model.beta, se, model.p_values):
print(f'{name:24s} beta = {b:7.2f} se = {s:5.2f} p = {p:.3f}')
Wedge gauge beta = -9.65 se = 4.14 p = 0.020
Interbelt gauge beta = -7.45 se = 2.88 p = 0.010
Peel force beta = -28.71 se = 10.86 p = 0.008
Wedge gauge×peel force beta = 19.10 se = 7.25 p = 0.008
The p_values are exactly the Wald tests \(2(1 - \Phi(|\beta/\text{se}|))\)
built from these standard errors. Because the model contains an interaction,
no single coefficient can be changed on its own — raising peel force also
raises the interaction column — so a hazard ratio is best computed between two
concrete tires. model.phi(Z) returns the multiplier \(e^{\beta'Z}\),
the hazard ratio against a tire with \(Z = 0\), whose hazard the baseline
h0 is, and the ratio of two multipliers is their hazard ratio at every
time. (The fit itself centres the covariates on their means, as R’s
coxph does, so a covariate far from zero cannot overflow
\(e^{\beta'Z}\); CoxPH.fit(..., center=True) keeps the baseline at
the means, model.center, and phi is then the hazard ratio against a
tire there, \(e^{\beta'(Z - \bar Z)}\).)
model.summary() gathers these into one table, as R’s summary(coxph)
and lifelines’ summary do: the coefficient, the hazard ratio
exp(coef), the standard error, 95% Wald intervals for both, z and the
p-value, one row per covariate (named by the columns for a model fitted with
fit_from_df). The model’s printout shows the same table;
summary(robust=True) uses the cluster-robust standard errors of
Cluster-robust standard errors instead. The parametric models’
summary() has the same columns, with the baseline distribution’s
parameters in rows of their own, above the coefficients.
model.summary().round(3)
| coef | exp(coef) | se(coef) | coef lower 95% | coef upper 95% | exp(coef) lower 95% | exp(coef) upper 95% | z | p | |
|---|---|---|---|---|---|---|---|---|---|
| covariate | |||||||||
| beta_0 | -9.647 | 0.000000e+00 | 4.142 | -17.766 | -1.528 | 0.000 | 2.170000e-01 | -2.329 | 0.020 |
| beta_1 | -7.446 | 1.000000e-03 | 2.882 | -13.095 | -1.797 | 0.000 | 1.660000e-01 | -2.583 | 0.010 |
| beta_2 | -28.714 | 0.000000e+00 | 10.859 | -49.997 | -7.431 | 0.000 | 1.000000e-03 | -2.644 | 0.008 |
| beta_3 | 19.103 | 1.977899e+08 | 7.248 | 4.897 | 33.309 | 133.867 | 2.922368e+14 | 2.636 | 0.008 |
Z_mean = Z.mean().values
for f in [0.9, 1.1]:
hr = model.phi(Z_mean * f) / model.phi(Z_mean)
print(f'every covariate x {f}: hazard ratio against the mean tire = '
f'{hr:.2f}')
every covariate x 0.9: hazard ratio against the mean tire = 14.43
every covariate x 1.1: hazard ratio against the mean tire = 0.07
Survival curves can be evaluated at any covariate value. Here we compare the mean tire against 10% above and below average:
plot_x = np.linspace(x.min(), x.max())
for f in [0.9, 1.0, 1.1]:
plt.step(plot_x, model.sf(plot_x, Z=Z_mean * f), label=f'{f:.0%}')
plt.legend(title='Covariate scale')
plt.xlabel('Survival time')
plt.ylabel('S(x)')
plt.show()
The step-function shape is the signature of the non-parametric baseline —
the model makes no smoothness assumptions about \(h_0(x)\). Keep the
evaluation times inside the observed range: the baseline is only estimated
there, and a Cox model’s hf returns the size of the baseline step at the
latest observed time, not a smooth hazard rate.
Tied event times
When failure times are recorded coarsely — to the day, the shift, the
inspection — several units share a time and the partial likelihood needs a tie
convention, chosen with tie_method=: 'breslow', 'efron', 'exact' or
'kalbfleisch-prentice' ('kp'). Every CoxPH fit defaults to Efron,
as R and lifelines do, and an Efron fit’s baseline hazard takes the same tie
correction (the covariate-weighted Fleming-Harrington estimator). Below, fifty units have
continuous lifetimes that were recorded only to the whole day, so up to six
share a day; each method is compared with the fit to the unrounded times,
which is the answer rounding took away:
rng = np.random.default_rng(0)
z_tie = rng.binomial(1, 0.5, 50).astype(float)
t_true = 10 * rng.weibull(2, 50) * np.exp(-0.7 * z_tie / 2) # true beta 0.7
t_day = np.ceil(t_true) # recorded to the whole day
counts = np.unique(t_day, return_counts=True)[1]
print('distinct days:', counts.size, ' largest tie:', counts.max())
print('unrounded times : beta = %.3f'
% CoxPH.fit(x=t_true, Z=z_tie).beta[0])
tie_beta = {}
for method in ['breslow', 'efron', 'exact', 'kalbfleisch-prentice']:
m = CoxPH.fit(x=t_day, Z=z_tie, tie_method=method)
tie_beta[method] = m.beta[0]
print(f'{method:22s} : beta = {m.beta[0]:.3f}')
distinct days: 19 largest tie: 6
unrounded times : beta = 1.073
breslow : beta = 0.989
efron : beta = 1.058
exact : beta = 1.065
kalbfleisch-prentice : beta = 1.127
Breslow’s approximation pulls the coefficient towards zero; Efron’s recovers
almost all of what rounding lost, and 'exact' — which averages the
likelihood over every order in which the tied failures could have happened —
lands closest to the unrounded fit, as it should for rounded continuous time.
'kalbfleisch-prentice' (alias 'kp') answers a different question: it
treats time as genuinely discrete, so its coefficient is a log odds ratio of
failing within a day rather than a log hazard ratio, and is larger here for
that reason, not because it is more accurate. Use it when time really is
discrete (a unit can only fail at an inspection), and compare it only with
other discrete-time fits. The two exact methods cost more than Efron — each
tie group is a recursion ('kalbfleisch-prentice') or a numerical integral
('exact') rather than a closed form — but both grow only polynomially with
the size of a tie group, so they remain practical on heavily tied data. With
no ties all four agree.
Delayed entry (left truncation)
A unit that only came under observation after it had already survived for a
while — equipment bought second-hand, patients enrolled some time after
diagnosis — must not be treated as if it had been watched from new: that would
credit it with survival it was never at risk of failing to show. Pass the entry
ages as tl to CoxPH.fit (for the parametric families, as the first
column of t, or tl_col in fit_from_df). Below, units whose failure
came before their entry age were never seen at all, as happens in practice:
rng = np.random.default_rng(1)
z_all = rng.normal(size=900)
T_all = 10 * rng.weibull(1.5, 900) * np.exp(-0.8 * z_all / 1.5) # beta = 0.8
entry_all = rng.uniform(0, 8, 900)
seen = T_all > entry_all # the rest failed before entry
x_le, z_le, tl_le = T_all[seen], z_all[seen].reshape(-1, 1), entry_all[seen]
t_le = np.column_stack([tl_le, np.full(seen.sum(), np.inf)])
naive = WeibullPH.fit(x=x_le, Z=z_le)
trunc = WeibullPH.fit(x=x_le, Z=z_le, t=t_le)
print('truth (alpha, beta, beta_0) : [10. 1.5 0.8]')
print('WeibullPH ignoring entry :', naive.params.round(3))
print('WeibullPH with truncation :', trunc.params.round(3))
print('Cox ignoring / with tl : %.3f / %.3f' % (
CoxPH.fit(x=x_le, Z=z_le).beta[0],
CoxPH.fit(x=x_le, Z=z_le, tl=tl_le).beta[0]))
truth (alpha, beta, beta_0) : [10. 1.5 0.8]
WeibullPH ignoring entry : [12.182 1.939 0.92 ]
WeibullPH with truncation : [10.204 1.563 0.827]
Cox ignoring / with tl : 0.878 / 0.837
Ignoring the entry ages badly distorts the baseline — the fitted life is too
long and the wear-out too steep, because the sample has been filtered towards
survivors — while the coefficient is affected less here because entry age is
unrelated to the covariate. If entry were related to the covariate the
coefficient would be biased too. Only left truncation is available for Cox; right or
interval truncation cannot be expressed in the forward partial likelihood and is
rejected, so use a parametric family (t=[tl, tr]) for those.
CoxPH.fit_from_df takes the entry ages as a column, tl_col:
import pandas as pd
entry_df = pd.DataFrame({'age': x_le, 'z': z_le[:, 0], 'entry': tl_le})
cox_df = CoxPH.fit_from_df(entry_df, x_col='age', Z_cols='z',
tl_col='entry', tie_method='breslow')
print('Cox with tl_col : %.3f' % cox_df.beta[0])
Cox with tl_col : 0.837
Fitting from a DataFrame: formulas and categorical covariates
Every family has fit_from_df, which names the columns of a pandas
DataFrame instead of passing arrays: x_col, c_col, n_col and the
covariates either as Z_cols (a list of numeric columns) or as a
formula (a formulaic
formula such as "age + site" or "age * site"). The fitted model
remembers its feature_names — and the formula’s encoding — so it can
predict directly from a DataFrame of raw covariates. Beyond those, the
parametric families take tl_col / tr_col (truncation) and init /
fixed; CoxPH.fit_from_df takes tl_col (delayed entry), tie_method
and strata_col; and the
frailty fitter requires a group_col. There is a single time column, so
interval-censored data (two time columns) go through fit.
A categorical covariate (a string or categorical column) is expanded with reference-level (treatment) coding: its first level is the baseline and each other level gets a coefficient measuring its effect relative to that reference. One level must be left out because the baseline hazard (or baseline distribution) already plays the role of the intercept — a column for every level would sum to one and be perfectly collinear with it, leaving the coefficients non-identified. Here three sites have different risks:
import pandas as pd
rng = np.random.default_rng(8)
n_df = 300
site = rng.choice(['A', 'B', 'C'], n_df)
age = rng.uniform(20, 60, n_df)
site_effect = np.select([site == 'A', site == 'B', site == 'C'],
[0.0, 0.5, -0.5]) # log-HR against site A
T_df = 50 * rng.weibull(2.0, n_df) * np.exp(-(0.03 * (age - 40)
+ site_effect) / 2.0)
cens_df = rng.uniform(20, 90, n_df)
patients = pd.DataFrame({'time': np.minimum(T_df, cens_df),
'censored': (T_df > cens_df).astype(int),
'age': age, 'site': site})
cox_df = CoxPH.fit_from_df(patients, x_col='time', c_col='censored',
formula='age + site')
for name, b in zip(cox_df.feature_names, cox_df.beta):
print(f'{name:10s} {b:6.3f}')
age 0.036
site[T.B] 0.538
site[T.C] -0.544
site[T.B] and site[T.C] are the log hazard ratios of sites B and C
against site A (true values 0.5 and -0.5), and age the log hazard ratio per
year (true 0.03). The same call works for the parametric families, and the
fitted model predicts from a DataFrame of raw covariates — the new rows are
encoded exactly as at fit time:
weib_df = WeibullPH.fit_from_df(patients, x_col='time', c_col='censored',
formula='age + site')
new = pd.DataFrame({'age': [40, 40, 40], 'site': ['A', 'B', 'C']})
print(weib_df.feature_names)
print('S(40) at age 40, sites A, B, C:',
weib_df.sf(np.full(3, 40.0), new).round(3))
['age', 'site[T.B]', 'site[T.C]']
S(40) at age 40, sites A, B, C: [0.541 0.352 0.699]
A formula beginning with 0 + asks for the full one-hot coding instead; with
a baseline distribution in the model that brings back the collinearity above,
so the last level is aliased (see below) and it is rarely what you want.
A covariate column the data cannot determine – a constant column (the
baseline is the intercept: Cox’s or Lin-Ying’s baseline hazard, the
Buckley-James intercept, or the scale of a family whose scale absorbs a
constant, as for Weibull PH or any AFT family), one
constant within each stratum of a stratified Cox fit, or a column that is a
linear combination of the others – is aliased, as R’s coxph and
lm do it: the fit runs on the other columns, whose estimates are what they
are without it, reports its coefficient, standard error and p-value as
nan, lists it in model.aliased, and predicts as though its coefficient
were 0. One warning names the columns. The columns are taken in order, so of
two collinear columns it is the later one that is aliased:
import warnings
doubled = np.column_stack([patients['age'], 2 * patients['age']])
with warnings.catch_warnings(record=True) as caught:
warnings.simplefilter('always')
aliased = CoxPH.fit(patients['time'], doubled, patients['censored'])
print(aliased.beta, aliased.aliased)
print(str(caught[0].message)[:60])
[0.03938244 nan] [1]
Covariate column(s) 1 of Z cannot be estimated: the partial
Wrapped categoricals (C(site), with
levels= or contrasts such as contr.sum) and data-dependent transforms
(scale(x), center(x), poly(x, 2), bs(x, df=3)) work too, and
are kept when the model is saved (see Saving and loading a fitted model).
A categorical level the model was not fitted with has no coefficient, so a
prediction for it is undefined: it raises a ValueError naming the column
and the level. So does a level declared with C(site, levels=[...]) (or
an unused category of a pd.Categorical column) that the fitted data has
no rows of: declaring the full list keeps the columns the same across data
splits, but nothing estimates that level’s coefficient, so the fit warns,
naming the level, and a prediction for it raises. (The level’s column is
all zeros, so it is aliased, without a second warning.) A
missing value is not a level: that row predicts nan, as for a missing
numeric covariate. This holds for every family that takes a formula,
before and after saving:
try:
weib_df.sf([40.0], pd.DataFrame({'age': [40], 'site': ['D']}))
except ValueError as err:
print(str(err).split('. ')[0])
print(weib_df.sf(np.full(2, 40.0),
pd.DataFrame({'age': [40, 40], 'site': ['B', None]})))
import warnings
with warnings.catch_warnings(record=True) as caught:
warnings.simplefilter('always')
declared = WeibullPH.fit_from_df(
patients, x_col='time', c_col='censored',
formula="age + C(site, levels=['A', 'B', 'C', 'D'])")
print([str(w.message).split('. ')[0] for w in caught
if 'no rows' in str(w.message)])
try:
declared.sf([40.0], pd.DataFrame({'age': [40], 'site': ['D']}))
except ValueError as err:
print(str(err).split('. ')[0])
Unknown categorical level(s): column 'site' has the level(s) ['D'], which are not among the levels ['A', 'B', 'C'] of the formula term 'site'
[0.35205175 nan]
['Categorical level(s) with no rows in the fitted data: column \'site\' has no rows at the level(s) [\'D\'] of the formula term "C(site, levels=[\'A\', \'B\', \'C\', \'D\'])"']
Unknown categorical level(s): column 'site' has the level(s) ['D'], which are not among the levels ['A', 'B', 'C'] of the formula term "C(site, levels=['A', 'B', 'C', 'D'])" (['D'] declared, but with no rows in the fitted data)
Time-Varying Covariates
Sometimes a covariate changes during a subject’s follow-up — a dose is
increased, a treatment begins, a machine is moved to a harsher environment.
Cox handles this through the counting-process (start-stop) format: each
subject contributes one row per interval \((\text{start}, \text{stop}]\)
on which its covariates are constant, and c is 0 only on the interval
that ends at the subject’s failure. Because each interval is exactly a
delayed-entry (left-truncated) observation, the partial likelihood fits this
format directly — splitting a subject into intervals with the same covariates
leaves the fit unchanged.
Use CoxPH.fit_tvc (arrays) or CoxPH.fit_tvc_from_df (a
start-stop DataFrame). The interval bounds follow surpyval’s xl / xr
naming and the status column c follows surpyval’s censoring convention —
c = 0 for the terminal event, c = 1 for a right-censored interval end
(a covariate change or administrative end). A subject may have at most one
c = 0 row, it must be its last interval, and its intervals must not overlap.
The covariates are Z_cols (numeric columns) or, as in fit_from_df, a
formula= that codes categorical columns such as the "yes" / "no"
columns of load_rossi_time_varying().
In the example below a covariate
stress switches from 0 to 1 at a random time for each unit and genuinely
raises the hazard once it turns on; units that fail before the switch
contribute a single interval, those that survive it contribute two:
from surpyval import CoxPH
from surpyval.univariate.regression import StepSchedule
rng = np.random.default_rng(0)
n, lam, beta = 2000, 0.5, 1.0
switch = rng.uniform(0.3, 1.5, size=n)
t_low = rng.exponential(1 / lam, size=n)
t_high = switch + rng.exponential(1 / (lam * np.exp(beta)), size=n)
T = np.where(t_low > switch, t_high, t_low)
rows = []
for i in range(n):
if T[i] <= switch[i]: # failed before the switch
rows.append((i, 0.0, T[i], 0, 0.0)) # c=0: event on one interval
else: # survived it: two intervals
rows.append((i, 0.0, switch[i], 1, 0.0)) # c=1: covariate change
rows.append((i, switch[i], T[i], 0, 1.0)) # c=0: terminal event
df = pd.DataFrame(rows, columns=['id', 'xl', 'xr', 'c', 'stress'])
model = CoxPH.fit_tvc_from_df(
df, i_col='id', xl_col='xl', xr_col='xr', c_col='c', Z_cols='stress',
)
model
Semi-Parametric Regression SurPyval Model
=========================================
Type : Proportional Hazards
Kind : Cox
Parameterization : Semi-Parametric
Data : 2000 units in 3298 start-stop intervals: 2000 events
Tie method : efron
Coefficients : exp(coef) is the hazard ratio; Wald 95% intervals
coef exp(coef) se(coef) lower 95% upper 95% z p
stress 1.052 2.863 0.07187 0.911 1.193 14.64 0
The fitted coefficient recovers the simulated log-hazard-ratio of the stress (\(\beta \approx 1\)) — a plain Cox fit that ignored the timing of the switch could not.
Writing intervals by hand is error-prone. A covariate timeline — one row per
covariate change per subject, each value holding until the subject’s next row,
the first row’s time being the entry and the last row carrying the exit time
and status — can be given instead with CoxPH.fit_tvc_timeline or
CoxPH.fit_tvc_timeline_from_df. The covariate on a subject’s last row
is ignored, as is c on every row but the last. It is expanded to the same
intervals, so the fit is identical:
timeline = []
for i in range(n):
timeline.append((i, 0.0, 0.0, 1)) # enters with stress 0
if T[i] > switch[i]:
timeline.append((i, switch[i], 1.0, 1)) # stress turns on
timeline.append((i, T[i], np.nan, 0)) # fails (Z ignored)
tl_df = pd.DataFrame(timeline, columns=['id', 'time', 'stress', 'c'])
print(tl_df.head(5))
model_tl = CoxPH.fit_tvc_timeline_from_df(
tl_df, i_col='id', x_col='time', Z_cols='stress', c_col='c',
)
print('same fit:', np.allclose(model_tl.beta, model.beta))
id time stress c
0 0 0.000000 0.0 1
1 0 0.893014 NaN 0
2 1 0.000000 0.0 1
3 1 0.623744 1.0 1
4 1 1.468383 NaN 0
same fit: True
Because survival now depends on the whole covariate path, evaluate it with
sf_tvc, describing the path as a
StepSchedule (or
(xl, Z) arrays). This is the same sf_tvc interface every regression
family exposes (see Time-varying covariates across families); with a single
constant segment it reduces exactly to sf. Here a unit stressed from
t = 1 onward has visibly lower survival than one never stressed:
t = np.linspace(0.01, 3.0, 100)
never = StepSchedule.constant([0.0])
late = StepSchedule.from_changepoints([0.0, 1.0], [[0.0], [1.0]])
plt.plot(t, model.sf_tvc(t, never), label='never stressed')
plt.plot(t, model.sf_tvc(t, late), label='stressed after t=1')
plt.legend()
plt.xlabel('Time')
plt.ylabel('S(t)')
plt.show()
The older interval-oriented
predict_tvc()
— which returns the survival and cumulative hazard at the baseline jump times
along a subject’s (xl, xr] intervals — remains available and agrees with
sf_tvc exactly.
Note
Predicting along a future covariate path. Some packages (lifelines, for
one) deliberately refuse to produce a survival curve from a
time-varying-covariate Cox model, reasoning that a subject’s future
covariate values are unknown: you cannot know \(Z(u)\) for \(u\) up
to a future time \(t\), so an observed subject’s survival past its
last record is undefined. SurPyval takes a different view because sf_tvc
answers a different question. You supply the covariate path as a plan or a
hypothesis — mission phases, a periodic duty cycle, a scheduled load
increase — and given that stated \(Z(\cdot)\) the survival
\(S(t \mid Z(\cdot))\) is exactly defined. This is scenario / what-if
evaluation under an assumed trajectory, not a claim to know the future; the
answer is only ever as good as the covariate plan you feed it. When the
future schedule is genuinely unknown, that is a reason not to specify one —
not a reason for the model to refuse an otherwise well-posed question.
Delayed entry is also supported here: a subject whose first interval starts after 0 simply enters the risk sets late, and gaps between a subject’s intervals are allowed (it is not at risk in the gap). The cluster-robust standard errors of a start-stop fit (below) cluster the rows by subject automatically.
Checking the proportional-hazards assumption
A Cox fit is only trustworthy if its central assumption holds: that each covariate multiplies the baseline hazard by a constant factor over time. When a coefficient is really drifting with time — a treatment that helps early but not late, say — the single number Cox reports is a time-average that can hide the effect entirely.
The standard check is the Grambsch-Therneau test, built on the scaled
Schoenfeld residuals. A fitted model exposes it through
check_ph().
It returns a table, as R’s cox.zph prints it: a 1-d.f. test for each
covariate and a joint GLOBAL test on the last row; a small p-value is
evidence against proportional hazards. We fit the tires model
with CoxPH.fit_from_df so the report carries the covariate names:
cols = ['Wedge gauge', 'Interbelt gauge', 'Peel force',
'Wedge gauge×peel force']
model = CoxPH.fit_from_df(tires, x_col='Survival', Z_cols=cols,
c_col='Censoring')
ph = model.check_ph()
ph.round(3)
| statistic | df | p | |
|---|---|---|---|
| covariate | |||
| Wedge gauge | 0.269 | 1 | 0.604 |
| Interbelt gauge | 0.033 | 1 | 0.855 |
| Peel force | 0.557 | 1 | 0.456 |
| Wedge gauge×peel force | 0.532 | 1 | 0.466 |
| GLOBAL | 1.148 | 4 | 0.887 |
Here every p-value is large, so there is no evidence against proportional
hazards — the Cox coefficients can be read as constant hazard ratios. (With 11
failures the test has little power, so “no evidence” is not strong evidence of
proportionality either.) The statistics match R’s cox.zph and lifelines,
including under Efron ties.
To see what a violation looks like, simulate a covariate whose effect reverses: exposed units have three times the baseline hazard before \(t = 0.5\) and a third of it afterwards. The Cox coefficient averages the two into something unremarkable, but the test is emphatic, and the scaled Schoenfeld residuals — which estimate \(\beta(t)\) at each failure — show the effect falling over time:
rng = np.random.default_rng(11)
z_rev = rng.binomial(1, 0.5, 400).astype(float)
e = rng.exponential(size=400) # unit-exponential "cumulative hazard"
# exposed: H(t) = 3t before t = 0.5, then 1.5 + (t - 0.5) / 3
t_exposed = np.where(e < 1.5, e / 3, 0.5 + 3 * (e - 1.5))
t_rev = np.where(z_rev == 1, t_exposed, e)
rev = pd.DataFrame({'x': np.minimum(t_rev, 3.0),
'c': (t_rev > 3.0).astype(int), 'z': z_rev})
m_rev = CoxPH.fit_from_df(rev, x_col='x', Z_cols='z', c_col='c')
print('averaged beta :', m_rev.beta.round(3))
for transform in ['km', 'rank', 'identity', 'log']:
p = m_rev.check_ph(transform=transform).loc['GLOBAL', 'p']
print(f'check_ph(transform={transform!r:10s}) p = {p:.1e}')
scaled = m_rev.compute_residuals('scaled_schoenfeld')[:, 0]
event_times = rev['x'][rev['c'] == 0].to_numpy() # same order as residuals
plt.plot(event_times, scaled, '.', alpha=0.4)
plt.axhline(m_rev.beta[0], color='k', label='fitted constant beta')
plt.xlabel('failure time'); plt.ylabel('scaled Schoenfeld residual')
plt.legend()
plt.show()
averaged beta : [0.394]
check_ph(transform='km' ) p = 5.5e-16
check_ph(transform='rank' ) p = 5.5e-16
check_ph(transform='identity') p = 1.1e-14
check_ph(transform='log' ) p = 2.9e-14
The residuals sit high early and low late, straddling the constant fit. The
transform argument chooses the function of time the residuals are tested
against: "km" (the default, \(1 -\) Kaplan-Meier) spreads the failures
evenly and is the usual choice; "rank", "identity" and "log" are the
alternatives from cox.zph. A violation like this one calls for
stratification (below), a time-varying covariate, or a different family.
The residuals underlying the test (and several others) are available directly
through
compute_residuals(),
with kind one of "schoenfeld", "scaled_schoenfeld",
"martingale", "deviance", "score" or "dfbeta". Schoenfeld
residuals come one row per failure, in the order the failures appear in the
fitted data (as in the plot above); the others one row per observation, in
input order (for a start-stop fit, in the fit’s internal order, sorted by
subject and entry time). They follow the tie method of a Breslow or Efron
fit; after an 'exact' or 'kalbfleisch-prentice' fit the Breslow forms
are used. Martingale
residuals plotted against a covariate reveal non-linear functional form;
deviance residuals highlight poorly-predicted individuals; dfbeta residuals
show how far each observation moves each coefficient:
martingale = model.compute_residuals('martingale')
print('martingale residuals sum to zero:',
np.isclose(martingale.sum(), 0.0))
dfbeta = model.compute_residuals('dfbeta')
most = np.abs(dfbeta).argmax(axis=0)
print('most influential tire per coefficient:', most)
martingale residuals sum to zero: True
most influential tire per coefficient: [10 13 13 13]
Schoenfeld, score and martingale residuals all sum to zero at the maximum of the partial likelihood — a useful sanity check that the fit has converged.
Cluster-robust standard errors
The model-based standard errors assume every observation is independent. When the data are clustered — several failures from the same machine, repeated events on the same subject, items drawn in grouped batches — that assumption is wrong and the naive errors are too small. The Lin-Wei sandwich (or “robust”) variance corrects for it, using the dfbeta residuals grouped by cluster.
Pass a cluster label per observation to
robust_summary()
(or robust_covariance for the matrix);
with no cluster argument each observation is its own cluster (the ordinary
robust variance):
summary = model.robust_summary()
for name, se, p in zip(summary['covariate'], summary['se'],
summary['p_value']):
print(f" {name:24s} robust SE = {se:8.3f} p = {p:.3f}")
Wedge gauge robust SE = 4.635 p = 0.037
Interbelt gauge robust SE = 2.610 p = 0.004
Peel force robust SE = 11.801 p = 0.015
Wedge gauge×peel force robust SE = 7.766 p = 0.014
To see why clustering matters, imagine every tire had been measured twice and
both rows entered the fit as if independent. Ignoring that would understate the
standard errors by exactly \(\sqrt{2}\); passing the shared cluster
label recovers the correct value. (We use Breslow ties so that duplicating the
data leaves the coefficients exactly unchanged; Efron’s correction treats the
duplicates as extra ties.)
twice = pd.concat([tires, tires], ignore_index=True)
tire_id = np.tile(np.arange(len(tires)), 2)
once = CoxPH.fit_from_df(tires, x_col='Survival', Z_cols=cols,
c_col='Censoring', tie_method='breslow')
dup = CoxPH.fit_from_df(twice, x_col='Survival', Z_cols=cols,
c_col='Censoring', tie_method='breslow')
naive_se = lambda m: np.sqrt(np.diag(np.linalg.inv(m.jac(m.beta)[1])))
print('naive SE, once / twice :', (naive_se(once) / naive_se(dup)).round(3))
print('robust SE, once :', once.robust_summary()['se'].round(3))
print('robust SE, twice, clustered by tire:',
dup.robust_summary(cluster=tire_id)['se'].round(3))
naive SE, once / twice : [1.414 1.414 1.414 1.414]
robust SE, once : [ 4.459 2.467 11.143 7.316]
robust SE, twice, clustered by tire: [ 4.459 2.467 11.143 7.316]
Stratified Cox models
When proportional hazards fails for a nuisance covariate — a study site, a batch, a device generation you do not want to model explicitly — the standard remedy is stratification: fit a separate baseline hazard for each stratum while sharing the coefficients \(\beta\). Risk sets never cross a stratum boundary, so the comparison is always within-stratum.
Pass strata (a label per observation) to CoxPH.fit, or
strata_col to CoxPH.fit_from_df. The example below is deliberately
adversarial: the baseline hazard differs by an order of magnitude across three
sites and the covariate is correlated with the site. An ordinary Cox fit is
badly confounded; the stratified fit recovers the true coefficient:
st_rng = np.random.default_rng(0)
n_st = 600
site = st_rng.integers(0, 3, n_st)
Z_st = st_rng.normal(site.astype(float), 1.0).reshape(-1, 1) # confounded
baseline = np.array([1.0, 6.0, 30.0])[site]
x_st = st_rng.exponential(baseline / np.exp(0.8 * Z_st[:, 0])) # true 0.8
c_st = (st_rng.random(n_st) < 0.15).astype(int)
pooled = CoxPH.fit(x=x_st, Z=Z_st, c=c_st)
stratified = CoxPH.fit(x=x_st, Z=Z_st, c=c_st, strata=site)
print(f"true beta = 0.80")
print(f"pooled = {pooled.beta[0]:.3f} (confounded by site)")
print(f"stratified = {stratified.beta[0]:.3f}")
true beta = 0.80
pooled = 0.069 (confounded by site)
stratified = 0.732
An observation whose stratum label is missing (None, NaN or pandas
NA) has no baseline to belong to, so it is dropped with a warning giving
the count, just as a row with a missing covariate is.
Prediction on a stratified model needs a stratum argument to pick the right
baseline — one of stratified.strata_labels, here the site codes 0, 1 and
2; sf, Hf, hf, ff and df all accept it, and refuse to
guess if it is left out. A missing label (NaN or NA) predicts nan,
as a missing covariate does:
stratified.sf(x=1.0, Z=[[0.0]], stratum=0)
np.float64(0.4554426602019849)
The residual diagnostics and robust errors above assume a single baseline, so
they are not available on a stratified model; the coefficients and their
model-based p-values are. A stratified model can also not be serialised.
Along a covariate path, sf_tvc, Hf_tvc and predict_tvc take the
same stratum argument.
Semi-Parametric — Buckley-James (AFT)
Cox leaves the baseline hazard unspecified; Buckley-James is its accelerated-time counterpart, leaving the error distribution unspecified:
with \(\varepsilon\) drawn from an arbitrary distribution estimated from the
data. It is fitted by the Buckley-James iteration — repeatedly imputing each
censored log-time by its conditional expectation under the Kaplan-Meier of the
current residuals, then re-fitting by least squares, until the coefficients stop
moving. Coefficients are reported in surpyval’s accelerated-failure sign, the
same as WeibullAFT: a positive coefficient shortens life (it is the
negative of the slope in the equation above).
from surpyval import BuckleyJames
# log T = 3 - 0.8 Z + noise, right-censored: higher Z shortens life, so
# the accelerated-failure coefficient is +0.8.
rng = np.random.default_rng(0)
Z_bj = rng.normal(size=(400, 1))
T_bj = np.exp(3.0 - 0.8 * Z_bj[:, 0] + rng.normal(0, 0.5, size=400))
cens = np.exp(3.6)
c_bj = (T_bj > cens).astype(int)
x_bj = np.minimum(T_bj, cens)
model = BuckleyJames.fit(x=x_bj, Z=Z_bj, c=c_bj)
model
Buckley-James AFT SurPyval Model
================================
Kind : Semi-Parametric AFT
Converged : True (12 iters)
Data : 400 units: 291 events at 291 unique times, 109 right censored
Coefficients (positive => accelerates failure):
beta_0 : 0.757063
The Converged line reports whether the iteration reached a fixed point
(model.converged and model.n_iter); the estimator can settle into a
two-point cycle, which surpyval detects and averages, and a fit that has not
converged within max_iter iterations (default 100, with step tolerance
tol=1e-5) warns. The coefficients are model.beta (also model.coef).
Only observed and right-censored data with positive times are accepted.
Buckley-James has no simple closed-form standard error, so uncertainty comes from a percentile bootstrap — resampling, refitting, and taking coefficient percentiles:
model.bootstrap_ci(n_boot=200, random_state=1)
array([[0.70921631, 0.8128302 ]])
Predictions use the fitted residual distribution directly,
\(S(t \mid Z) = S_\varepsilon(\log t + \beta' Z)\), so the survival curves
shift with the covariate (sf, ff and Hf are available; being a
step function, the model has no density or hazard rate). Each call takes a
single covariate row and returns the curve at every time given:
t = np.linspace(1, 60, 200)
for z in [-1.0, 0.0, 1.0]:
plt.step(t, model.sf(t, Z=[z]), where='post', label=f'Z = {z:g}')
plt.legend()
plt.xlabel('Time')
plt.ylabel('S(t)')
plt.show()
BuckleyJames.fit_from_df accepts Z_cols or a formula, exactly as the
other families do.
Semi-Parametric — Additive Hazards
The Lin & Ying additive hazards model is the additive-scale companion to Cox. Like Cox it leaves the baseline hazard \(h_0(x)\) completely unspecified, but the covariate effect is a risk difference rather than a hazard ratio:
Its practical convenience is that, unlike Cox’s iterative partial likelihood, the coefficient estimator is closed form — a ratio of sums accumulated over the risk sets — so there is no optimisation and nothing to converge. Standard errors come from the Lin-Ying sandwich estimator. We can reuse the tire data and the significant covariates from the Cox fit above:
from surpyval import AdditiveHazards
model = AdditiveHazards.fit(x=x, Z=Z, c=c)
model
Semi-Parametric Regression SurPyval Model
=========================================
Type : Additive Hazards
Kind : Lin-Ying
Parameterization : Semi-Parametric
Parameters :
beta_0 : -2.2060322340265235
beta_1 : -1.0635504695453748
beta_2 : -2.6523846942293723
beta_3 : 2.4439071709451117
print(model.p_values)
print(model.standard_errors())
[0.03400316 0.03996378 0.01927683 0.01830796]
[1.04056444 0.51776374 1.13343292 1.0358479 ]
The coefficients read as risk differences: a one-unit change in a covariate shifts the absolute hazard by \(\beta\) at every time. As in the Cox fit, higher gauge and peel-force values reduce the hazard (improving life) while the interaction term counteracts — the same story, told on the additive scale.
The model’s Hf, sf and ff use the step baseline
\(\hat H_0(t) + t\,\beta'Z\). A hazard rate needs a smooth baseline, so
hf (and df) kernel-smooth the baseline increments; the bandwidth
argument of hf controls the smoothing (by default a normal-reference rule on
the event times), and estimates near the ends of the observed range are
attenuated. The additive model has no multiplier, so phi() is not defined
for it.
Note
An additive hazard can go negative when \(\beta' Z\) is sufficiently
negative — nothing constrains \(h_0(x) + \beta' Z > 0\) — and the
Lin-Ying estimate also dips between the event times, where its baseline
drifts down by \(\beta'\bar Z(t)\). AdditiveHazards therefore
predicts with the running maximum of its estimate from time 0: Hf is
held where the estimate falls, hf is 0 there, and sf stays in
\([0, 1]\) and never rises. Where Hf is held for long (a flat
survival curve), the model says the covariate value’s hazard is negative
there: the additive model is a poor description at that covariate value,
or you are outside the range where it is well behaved. When covariate
effects are strongly protective, a proportional-hazards model — whose
exponential form keeps the hazard positive — is often the safer choice.
The parametric AH models below keep their own values where the
hazard is negative (survival above 1, a negative density), and every
prediction there warns once that it is so.
Just as Cox has parametric proportional-hazards counterparts (the next
section), there is also a parametric additive-hazards model — a parametric
baseline hazard with the same additive covariate term, h(x|Z) = h_0(x;θ) +
β'Z, fit by maximum likelihood. It is available as the AH(distribution)
factory and as pre-built WeibullAH, ExponentialAH, … instances, and
gives a smooth, extrapolatable version of what AdditiveHazards estimates
non-parametrically. Below, a Weibull wear-out hazard has an exposure that adds
0.05 failures per unit time on top of it:
from surpyval import WeibullAH
rng = np.random.default_rng(4)
z_ah = rng.binomial(1, 0.5, 500).astype(float)
# H(x) = (x / 10)^2 + 0.05 z x -> solve H(x) = E for a unit exponential E
E = rng.exponential(size=500)
x_ah = (-0.05 * z_ah + np.sqrt((0.05 * z_ah) ** 2 + 0.04 * E)) / 0.02
c_ah = (x_ah > 25).astype(int)
x_ah = np.minimum(x_ah, 25)
wah = WeibullAH.fit(x=x_ah, Z=z_ah.reshape(-1, 1), c=c_ah)
print('WeibullAH (alpha, beta, beta_0):', wah.params.round(3))
print('standard errors :', wah.standard_errors().round(3))
ly = AdditiveHazards.fit(x=x_ah, Z=z_ah.reshape(-1, 1), c=c_ah)
print('Lin-Ying beta_0 = %.3f (se %.3f)' % (ly.beta[0], ly.se[0]))
WeibullAH (alpha, beta, beta_0): [9.659 1.997 0.035]
standard errors : [0.321 0.09 0.009]
Lin-Ying beta_0 = 0.046 (se 0.011)
Both estimators recover the risk difference to within about two standard errors (0.035 and 0.046 against a true 0.05), and the Weibull baseline (scale 10, shape 2) is recovered too. The parametric fit is more efficient when its baseline is right; Lin-Ying makes no assumption about the baseline.
The positivity caveat above bites differently here. The likelihood needs
log(h) at every failure, so the optimiser only accepts parameter values
that keep \(h_0(x) + \beta'Z\) positive at every observed failure. When
the data would prefer a negative hazard — a strongly protective covariate —
the fit returns the best model that stays positive, pressed against that
boundary (the fitted hazard of the protected units is then close to zero at
their earliest failures, and the baseline is bent to compensate), and warns
that it has done so; it raises only if the optimiser cannot end at a
positive-hazard point. Treat that warning as a verdict on the model, not the
optimiser, and prefer a proportional-hazards model when effects are strongly
protective.
Parametric Proportional Hazards (PH)
When you are confident about the shape of the baseline distribution — or when you need to extrapolate beyond the observed time range — a fully parametric PH model is preferable. It estimates the same β as the Cox model, but also estimates the baseline distribution parameters, giving a smooth, continuous survival function.
SurPyval provides pre-built parametric PH instances for every standard
distribution: ExponentialPH, NormalPH, WeibullPH, GumbelPH,
LogisticPH, LogNormalPH, and GammaPH. The Weibull is the most
common choice in reliability engineering — its shape parameter lets it capture
increasing, constant, or decreasing hazard rates.
from surpyval import WeibullPH
model = WeibullPH.fit(x=x, Z=Z, c=c)
model
Parametric Regression SurPyval Model
====================================
Kind : Proportional Hazard
Distribution : Weibull
Regression Model : Log Linear [e^(beta'Z)]
Fitted by : MLE
Data : 34 units: 11 events at 8 unique times, 23 right censored
Baseline : Weibull parameters; Wald 95% intervals
estimate se lower 95% upper 95%
alpha 0.2426 0.08144 0.1256 0.4684
beta 16.06 3.951 9.914 26.01
Coefficients : exp(coef) is the hazard ratio; Wald 95% intervals
coef exp(coef) se(coef) lower 95% upper 95% z p
beta_0 -9.165 0.0001046 3.724 -16.46 -1.867 -2.461 0.01384
beta_1 -7.999 0.0003359 2.812 -13.51 -2.487 -2.845 0.004447
beta_2 -27.5 1.136e-12 9.537 -46.19 -8.812 -2.884 0.003927
beta_3 18.39 9.654e+07 6.422 5.798 30.97 2.863 0.004199
Notice the coefficients are close to the Cox model’s, each within 10% of it —
this is expected when the Weibull is a reasonable fit to the baseline. The parameters are listed in
the order model.parameter_names gives: the distribution’s own parameters
first, then one beta_j per covariate column.
If none of the pre-built distributions suit your data, the PH factory creates
a parametric PH model for any surpyval distribution:
from surpyval import LogNormal
from surpyval import PH
model = PH(LogNormal).fit(x=x, Z=Z, c=c)
model
Parametric Regression SurPyval Model
====================================
Kind : Proportional Hazard
Distribution : LogNormal
Regression Model : Log Linear [e^(beta'Z)]
Fitted by : MLE
Data : 34 units: 11 events at 8 unique times, 23 right censored
Baseline : LogNormal parameters; Wald 95% intervals
estimate se lower 95% upper 95%
mu -0.02098 0.3008 -0.6105 0.5685
sigma 0.1549 0.08069 0.05579 0.43
Coefficients : exp(coef) is the hazard ratio; Wald 95% intervals
coef exp(coef) se(coef) lower 95% upper 95% z p
beta_0 -0.3739 0.6881 2.856 -5.972 5.225 -0.1309 0.8959
beta_1 -0.3814 0.6829 1.657 -3.629 2.866 -0.2302 0.818
beta_2 0.9695 2.637 3.699 -6.281 8.22 0.2621 0.7933
beta_3 -2.152 0.1163 2.986 -8.004 3.7 -0.7207 0.4711
The log-normal is a natural choice when the log of the survival time is expected to be normally distributed — common in medical and biological data. (A log-normal PH model multiplies a log-normal hazard; it is not the same as the log-normal AFT model below, which is linear regression on log time.)
Fixed parameters, censoring and simulation
The PH, AFT, PO, AH and accelerated-life fitters accept fixed={name: value}
to hold any parameter —
a distribution parameter or a coefficient — at a known value, for instance a
Weibull shape known from experience with the failure mode. The fixed parameter
is excluded from the covariance (its standard error is zero):
fixed_shape = WeibullPH.fit(x=x, Z=Z, c=c, fixed={'beta': 15})
print(fixed_shape.parameter_names)
print(fixed_shape.params.round(3))
print(fixed_shape.standard_errors().round(3))
['alpha', 'beta', 'beta_0', 'beta_1', 'beta_2', 'beta_3']
[ 0.238 15. -8.628 -7.618 -25.952 17.27 ]
[0.082 0. 3.098 2.384 7.381 4.772]
The parametric fitters use surpyval’s full likelihood, so every observation
type can be mixed in one fit: c = -1 (left censored), c = 2 (interval
censored, with x given as [left, right] pairs) and truncation through
t. Here the demo data set from the top of the page is re-recorded as an
inspection study — each unit checked every 2 time units, so a failure is only
known to lie in the interval between inspections — and the interval-censored fit
recovers the same parameters as the exact times:
left_edge = np.floor(x_demo / 2) * 2
x_int = np.column_stack([left_edge, left_edge + 2])
x_int[left_edge == 0, 0] = 1e-6 # first interval starts at ~0
c_int = np.full(len(x_demo), 2)
interval_fit = WeibullPH.fit(x=x_int, Z=Z_demo, c=c_int)
print('exact times :', demo.params.round(3))
print('inspection intervals:', interval_fit.params.round(3))
exact times : [10.339 2.032 0.725]
inspection intervals: [10.319 1.969 0.716]
phi(Z) returns the fitted hazard multiplier \(e^{\beta'Z}\) (for AFT
it is the acceleration factor, for PO the odds multiplier, for accelerated
life the modelled life; an additive model has none). random(size, Z)
draws lifetimes from the fitted model — useful for simulation studies and for
checking a fit against its own simulated data. It exists for the PH, parametric
AH and accelerated-life families (not AFT or PO). A PH or parametric AH model
returns size draws for each covariate row, in the order given, together
with the matching covariate rows; an accelerated life model does the same for
each distinct stress, in sorted order. As everywhere in SurPyval,
random_state seeds the draw (None, the default, draws from numpy’s
global generator, so np.random.seed reproduces it):
print('hazard multipliers at Z = 0, 1:', demo.phi([[0.0], [1.0]]).round(3))
sim_x, sim_Z = demo.random(5, [[0.0], [1.0]], random_state=0)
print(sim_x.round(2))
print(sim_Z.ravel())
hazard multipliers at Z = 0, 1: [1. 2.065]
[ 6.99 11.81 18.31 20.71 4.76 2.23 5.15 4.1 5.67 1.92]
[0. 0. 0. 0. 0. 1. 1. 1. 1. 1.]
A custom covariate function
The log-linear \(e^{\beta'Z}\) is a choice, not a requirement. Some fields
use other forms — radiation epidemiology, for example, often models an
excess relative risk that grows linearly with dose,
\(\phi(z) = 1 + \beta z\), so that \(\beta\) is the extra risk per unit
dose. ProportionalHazardsFitter builds a PH fitter around any
\(\phi(Z, *\text{params})\) written with autograd.numpy. Its arguments
are a name, the baseline distribution, phi, a display name for it, the
parameter bounds and the parameter-name map (each either fixed or a function of
Z), and optionally a starting value. The bounds are how you keep
\(\phi\) positive — here \(\beta > 0\):
import autograd.numpy as anp
from surpyval import ProportionalHazardsFitter, Weibull
def linear_rr(Z, *params): # phi(Z) = 1 + beta'Z
return 1.0 + anp.dot(Z, anp.array(params))
WeibullLinearRR = ProportionalHazardsFitter(
'WeibullLinearRR', Weibull, linear_rr, "Linear [1 + beta'Z]",
phi_bounds=lambda Z: ((0, None),) * Z.shape[1],
phi_param_map=lambda Z: {f'beta_{i}': i for i in range(Z.shape[1])},
phi_init=lambda Z: np.full(Z.shape[1], 0.5),
)
rng = np.random.default_rng(0)
dose = rng.uniform(0, 4, 400)
# Weibull(10, 2) baseline; each unit of dose adds 50% to the hazard
x_rr = 10 * (-np.log(rng.uniform(size=400)) / (1 + 0.5 * dose)) ** (1 / 2)
rr = WeibullLinearRR.fit(x=x_rr, Z=dose)
print('alpha, beta, beta_0:', rr.params.round(3))
print('standard errors :', rr.standard_errors().round(3))
alpha, beta, beta_0: [10.15 2.081 0.604]
standard errors : [0.717 0.081 0.168]
The excess relative risk per unit dose, 0.60 (standard error 0.17), is within
one standard error of the true 0.5. Everything else — predictions, bounds,
fit_from_df — works as for the pre-built models, but a custom covariate
function cannot be rebuilt from a name, so such a model cannot be serialised.
Accelerated Failure Time (AFT)
The AFT model has a different and often more interpretable structure than PH. Rather than saying “this covariate increases your hazard rate by X%”, it says “this covariate makes you age X% faster”. Formally, if the baseline survival time is \(T_0\), then the survival time given covariates is:
A positive \(\beta_j\) means covariate \(z_j\) compresses time (accelerates failure). A negative \(\beta_j\) stretches time (prolongs life). The median survival time simply scales by \(e^{-\beta' Z}\) — a direct and intuitive interpretation.
The relationship to the cumulative hazard is:
An important practical note: the Weibull (and its special case the Exponential) is the one distribution for which AFT and PH are the same model. A Weibull-AFT and a Weibull-PH fit to the same data have identical likelihoods, and their coefficients differ only by the Weibull shape, \(\beta_{PH} = \text{shape} \times \beta_{AFT}\). For every other distribution — log-normal included — AFT and PH are genuinely distinct models.
from surpyval import WeibullAFT
model = WeibullAFT.fit(x=x, Z=Z, c=c)
model
Parametric Regression SurPyval Model
====================================
Kind : Accelerated Failure Time
Distribution : Weibull
Regression Model : Log Linear [exp(beta'Z)]
Fitted by : MLE
Data : 34 units: 11 events at 8 unique times, 23 right censored
Baseline : Weibull parameters; Wald 95% intervals
estimate se lower 95% upper 95%
alpha 0.2426 0.08144 0.1256 0.4684
beta 16.06 3.951 9.914 26.01
Coefficients : exp(coef) is the acceleration factor; Wald 95% intervals
coef exp(coef) se(coef) lower 95% upper 95% z p
beta_0 -0.5708 0.5651 0.1953 -0.9536 -0.1879 -2.922 0.003479
beta_1 -0.4981 0.6077 0.154 -0.8 -0.1962 -3.234 0.001223
beta_2 -1.713 0.1804 0.4739 -2.642 -0.784 -3.614 0.0003012
beta_3 1.145 3.142 0.3041 0.5489 1.741 3.765 0.0001665
Checking the Weibull equivalence numerically against the WeibullPH fit:
ph_fit = WeibullPH.fit(x=x, Z=Z, c=c)
shape = model.params[1]
print('shape x AFT coefficients:', (shape * model.params[2:]).round(3))
print('PH coefficients :', ph_fit.params[2:].round(3))
print('neg log-likelihoods :', round(model.neg_ll(), 4),
round(ph_fit.neg_ll(), 4))
shape x AFT coefficients: [ -9.165 -7.999 -27.503 18.385]
PH coefficients : [ -9.165 -7.999 -27.503 18.385]
neg log-likelihoods : -6.0204 -6.0204
The AFT factory works with any distribution. Log-Normal AFT is a
particularly common choice — it corresponds to ordinary linear regression on
\(\log T\) with censored observations, and is the parametric counterpart
of the Buckley-James model above:
from surpyval import LogNormalAFT
model_ln = LogNormalAFT.fit(x=x, Z=Z, c=c)
model_ln
Parametric Regression SurPyval Model
====================================
Kind : Accelerated Failure Time
Distribution : LogNormal
Regression Model : Log Linear [exp(beta'Z)]
Fitted by : MLE
Data : 34 units: 11 events at 8 unique times, 23 right censored
Baseline : LogNormal parameters; Wald 95% intervals
estimate se lower 95% upper 95%
mu -1.795 0.4082 -2.595 -0.9947
sigma 0.1104 0.02259 0.07393 0.1649
Coefficients : exp(coef) is the acceleration factor; Wald 95% intervals
coef exp(coef) se(coef) lower 95% upper 95% z p
beta_0 -0.6283 0.5335 0.2483 -1.115 -0.1416 -2.53 0.0114
beta_1 -0.6613 0.5162 0.222 -1.097 -0.2261 -2.978 0.002898
beta_2 -2.054 0.1283 0.5439 -3.12 -0.9877 -3.776 0.0001594
beta_3 1.307 3.694 0.352 0.6169 1.996 3.713 0.0002051
For distributions not in the pre-built list:
from surpyval import Gamma
from surpyval import AFT
model = AFT(Gamma).fit(x=x, Z=Z, c=c)
model
Parametric Regression SurPyval Model
====================================
Kind : Accelerated Failure Time
Distribution : Gamma
Regression Model : Log Linear [exp(beta'Z)]
Fitted by : MLE
Data : 34 units: 11 events at 8 unique times, 23 right censored
Baseline : Gamma parameters; Wald 95% intervals
estimate se lower 95% upper 95%
alpha 87.71 35.3 39.86 193
beta 503.6 265.7 179 1416
Coefficients : exp(coef) is the acceleration factor; Wald 95% intervals
coef exp(coef) se(coef) lower 95% upper 95% z p
beta_0 -0.623 0.5363 0.246 -1.105 -0.1409 -2.533 0.01132
beta_1 -0.6413 0.5266 0.2167 -1.066 -0.2165 -2.959 0.003085
beta_2 -2.009 0.1342 0.5444 -3.076 -0.9418 -3.69 0.0002244
beta_3 1.285 3.613 0.3524 0.5939 1.975 3.645 0.0002671
We can visualise the “time shift” interpretation by plotting survival curves. The entire curve moves left (shorter life) or right (longer life) on the time axis as the covariates change — a hallmark of the AFT model:
model = WeibullAFT.fit(x=x, Z=Z, c=c)
Z_mean = Z.mean().values
plot_x = np.linspace(x.min(), x.max())
for f in [0.9, 1.0, 1.1]:
plt.plot(plot_x, model.sf(plot_x, Z=Z_mean * f), label=f'{f:.0%}')
plt.legend(title='Covariate scale')
plt.xlabel('Survival time')
plt.ylabel('S(x)')
plt.show()
Compare this with the Cox PH survival curves above. Under PH the curves are powers of one another, \(S(x \mid Z) = S_0(x)^{\phi(Z)}\); under AFT they are the same curve shifted along the log-time axis, so on a log-time plot they are parallel. Neither family lets two curves cross. (For these Weibull fits the two descriptions coincide, as shown above.)
Proportional Odds (PO)
The proportional odds model is less common than PH or AFT but has an important niche: it is the right model when you believe the relative odds of failure are constant across covariate values, rather than the relative hazard rates.
The survival odds at time \(x\) are \(O(x) = S(x) / F(x)\). A PO model assumes these are scaled by \(\phi(Z)\):
Rearranging, the survival function is:
Because \(e^{\beta'Z}\) multiplies the odds of surviving, a positive coefficient here means a longer life — the opposite of the PH and AFT sign convention. To compare with a PH fit, negate the PO coefficients.
The key difference from PH becomes clear at long follow-up: as \(x \to \infty\), \(S_0(x) \to 0\), and the ratio \(F_0 + \phi S_0 \to 1\), so the covariate effect fades away. Everyone eventually fails, and the PO model respects that by letting the hazard ratio converge to 1 over time. Under PH, the hazard ratio is constant forever — a stronger and often unrealistic assumption for long studies.
PO is the natural companion to the Logistic and Log-Logistic distributions: with a logistic baseline the odds multiplier is a shift of the location, and with a log-logistic baseline it is a rescaling of time, so the model is also an AFT model (see the Regression Analysis page).
from surpyval import LogisticPO
model = LogisticPO.fit(x=x, Z=Z, c=c)
model
Parametric Regression SurPyval Model
====================================
Kind : Proportional Odds
Distribution : Logistic
Regression Model : Log Linear [exp(beta'Z)]
Fitted by : MLE
Data : 34 units: 11 events at 8 unique times, 23 right censored
Baseline : Logistic parameters; Wald 95% intervals
estimate se lower 95% upper 95%
mu -0.5108 0.3816 -1.259 0.2371
sigma 0.05339 0.01308 0.03304 0.0863
Coefficients : exp(coef) is the survival odds ratio; Wald 95% intervals
coef exp(coef) se(coef) lower 95% upper 95% z p
beta_0 11.28 7.939e+04 4.618 2.231 20.33 2.443 0.01456
beta_1 9.966 2.129e+04 3.705 2.704 17.23 2.69 0.007152
beta_2 32.98 2.102e+14 11.11 11.21 54.75 2.969 0.00299
beta_3 -21.83 3.322e-10 7.438 -36.4 -7.246 -2.934 0.003345
The tire coefficients have the opposite signs to the PH and AFT fits — the same
conclusions, in the survival-odds convention. The PO factory accepts any
distribution:
from surpyval import Weibull
from surpyval import PO
model = PO(Weibull).fit(x=x, Z=Z, c=c)
model
Parametric Regression SurPyval Model
====================================
Kind : Proportional Odds
Distribution : Weibull
Regression Model : Log Linear [exp(beta'Z)]
Fitted by : MLE
Data : 34 units: 11 events at 8 unique times, 23 right censored
Baseline : Weibull parameters; Wald 95% intervals
estimate se lower 95% upper 95%
alpha 0.00952 0.01455 0.0004763 0.1903
beta 0.7047 0.1846 0.4217 1.177
Coefficients : exp(coef) is the survival odds ratio; Wald 95% intervals
coef exp(coef) se(coef) lower 95% upper 95% z p
beta_0 10.75 4.662e+04 4.533 1.865 19.63 2.371 0.01773
beta_1 9.357 1.158e+04 3.624 2.255 16.46 2.582 0.009816
beta_2 30.88 2.578e+13 10.72 9.873 51.89 2.881 0.003963
beta_3 -20.47 1.283e-09 7.203 -34.59 -6.356 -2.842 0.004477
The Weibull baseline tells the same story as the logistic one: every coefficient has the opposite sign to the PH fit, as the survival-odds convention requires. With 11 failures and four covariates, though, none of these fits is well determined, and Model Selection below shows the data cannot separate the PO description from the PH/AFT one.
The fading effect is easiest to see on data simulated from a PO model. Below, a log-logistic baseline has its survival odds multiplied by \(e^{1}\) for exposed units. The hazard ratio of exposed to unexposed units starts near \(e^{-1} \approx 0.37\) and climbs towards 1:
from surpyval import LogLogistic
rng = np.random.default_rng(2)
z_po = rng.binomial(1, 0.5, 400).astype(float)
U = rng.uniform(size=400)
# invert F(x|z) = U for survival odds (x/10)^-3 * e^z
x_po = 10 * (np.exp(z_po) * U / (1 - U)) ** (1 / 3)
po = PO(LogLogistic).fit(x=x_po, Z=z_po.reshape(-1, 1))
print('alpha, beta, beta_0:', po.params.round(3))
times = np.array([1.0, 5.0, 10.0, 20.0, 40.0])
print('hazard ratio at', times, ':',
(po.hf(times, [1.0]) / po.hf(times, [0.0])).round(3))
alpha, beta, beta_0: [10.004 3.087 1.035]
hazard ratio at [ 1. 5. 10. 20. 40.] : [0.355 0.381 0.524 0.839 0.976]
A practical rule of thumb: if the Kaplan-Meier curves for different covariate
groups converge at long times (rather than remaining parallel on the log-hazard
scale), PO is likely a better fit than PH. Proportional odds has no
random.
A fitted PO model can be evaluated along a covariate path that changes over
time. Its hazard, \(h_0(x) / (F_0(x) + e^{\beta' Z} S_0(x))\), depends
only on the time and the covariate value at that time, so sf_tvc adds up
the cumulative-hazard increment of each constant-covariate segment exactly as
it does for PH, and for the same reason fit_tvc fits start-stop data (see
Time-varying covariates across families). Here the
unit is unexposed until \(t = 10\) and exposed afterwards:
from surpyval.univariate.regression import StepSchedule
switch = StepSchedule.from_changepoints([0.0, 10.0], [[0.0], [1.0]])
at = np.array([5.0, 10.0, 15.0, 30.0])
print('never exposed :', po.sf(at, [0.0]).round(3))
print('exposed from 10:', po.sf_tvc(at, switch).round(3))
print('always exposed :', po.sf(at, [1.0]).round(3))
never exposed : [0.895 0.5 0.223 0.033]
exposed from 10: [0.895 0.5 0.303 0.059]
always exposed : [0.96 0.738 0.446 0.087]
Up to \(t = 10\) the unit follows the unexposed curve. After the switch its hazard becomes the exposed one, but its survival does not jump up to the exposed curve: it starts from the survival it has already reached (0.5) and falls more slowly from there, ending between the other two curves.
Semi-Parametric — Proportional Odds
ProportionalOdds is to the PO models what CoxPH is to the PH ones:
the covariates multiply the survival odds as above, but the baseline is left
to the data. Its failure odds \(G_0(x) = F_0(x) / S_0(x)\) are a
non-decreasing step function that jumps at the event times:
Unlike Cox’s, this baseline does not drop out of the likelihood, so it is
estimated together with \(\beta\) by nonparametric maximum likelihood
(Murphy, Rossini and van der Vaart, 1997): for each \(\beta\) the jumps
of \(G_0\) are solved for exactly, and \(\beta\) maximises the
resulting profile likelihood, whose curvature gives the standard errors.
The model takes observed and right-censored data, with left truncation
through a 1-D tl. The sign convention is the parametric PO models’: a
positive coefficient means a longer life. (R’s timereg::prop.odds and
mets::logitSurv model the odds of failure, so their coefficients are the
negatives of these.)
Fitted to the log-logistic data above, it recovers the coefficient the
parametric PO(LogLogistic) fit found, with nearly the same standard
error, and a baseline that tracks the log-logistic one without assuming it:
from surpyval import ProportionalOdds
spo = ProportionalOdds.fit(x=x_po, Z=z_po.reshape(-1, 1))
spo
Semi-Parametric Regression SurPyval Model
=========================================
Type : Proportional Odds
Kind : NPMLE (Murphy, Rossini & van der Vaart)
Parameterization : Semi-Parametric
Data : 400 units: 400 events at 400 unique times
Coefficients : exp(coef) is the survival odds ratio; Wald 95% intervals
coef exp(coef) se(coef) lower 95% upper 95% z p
beta_0 1.007 2.738 0.1792 0.6562 1.359 5.622 1.887e-08
print('semi-parametric beta_0: %.3f (se %.3f)' % (spo.beta[0], spo.se[0]))
print('PO(LogLogistic) beta_0: %.3f (se %.3f)'
% (po.params[2], po.standard_errors()[2]))
t = np.array([5.0, 10.0, 20.0])
print('baseline survival, semi-parametric:', spo.sf(t, [0.0]).round(3))
print('baseline survival, log-logistic :', po.sf(t, [0.0]).round(3))
semi-parametric beta_0: 1.007 (se 0.179)
PO(LogLogistic) beta_0: 1.035 (se 0.180)
baseline survival, semi-parametric: [0.911 0.495 0.115]
baseline survival, log-logistic : [0.895 0.5 0.105]
Both put the true coefficient of 1 well inside their intervals (1.007 with
a standard error of 0.179, against 1.035 and 0.180), and the two baseline
survival curves are within 0.02 of each other and 0.025 of the truth
\(1 / (1 + (x/10)^3)\): 0.889, 0.5 and 0.111. The parametric fit is the
better choice when its baseline is right; ProportionalOdds is the one to
use when the shape of the baseline is what you cannot commit to.
The fitted model has the Cox model’s interface: sf, ff and Hf at
covariate rows (grid=True for a curve per row), hf and df as the
jumps at the event times, summary(), param_cb, concordance and
to_dict / from_dict; fit_from_df takes Z_cols or a
formula. Before the first event time the survival is 1, and after the
last observed time it holds its last value. A covariate that separates the
events from the survivors (a level with no events) leaves the likelihood
without a maximum, and the fit warns “No finite maximum”.
Confidence Bounds
A point estimate is only half the story. The parametric regression models (PH,
AFT, PO, AH and AL) carry the full parameter covariance — the inverse of the
Hessian of the negative log-likelihood at the fit — so every coefficient and every
predicted curve comes with an interval. The Hessian is the exact one the fit
computes (with autograd) to check that it reached a maximum; a model without one
(an accelerated-life fit, an AFT time-varying fit, a fit with no finite maximum,
or one whose Hessian is not positive definite there) uses a numerical Hessian
instead. After a fit, the parameter covariance (covariance()) and standard
errors are available directly, and are computed once:
from surpyval import WeibullPH
rng = np.random.default_rng(0)
Z_cb = rng.normal(size=(300, 1))
x_cb = 10.0 * (
-np.log(rng.uniform(size=300)) / np.exp(Z_cb[:, 0] * 0.8)
) ** (1 / 2.0)
m_cb = WeibullPH.fit(x=x_cb, Z=Z_cb, c=np.zeros(300, dtype=int))
print(m_cb.parameter_names)
print(m_cb.standard_errors())
['alpha', 'beta', 'beta_0']
[0.27985641 0.09444288 0.0693556 ]
param_cb gives a Wald confidence bound on a single parameter, computed on a
scale chosen from the parameter’s support (log for a positive scale, natural for
an unbounded coefficient) so the interval always stays valid:
m_cb.param_cb('beta_0') # 95% CI for the covariate coefficient
array([0.70158804, 0.973457 ])
cb propagates the parameter covariance through a predicted function by the
delta method, returning a confidence band. The band on sf, ff and
Hf is formed on the baseline family’s probability-plot scale, as for the
univariate models (Parametric Estimation): \(\ln H\) for a
Weibull, Exponential, Rayleigh or Gumbel baseline, the normal quantile of
\(F\) for a Normal or LogNormal one, and the logit of \(F\) for the
rest. A model with its coefficients fixed at 0 then gives the univariate band,
and the band rises with time wherever the shape’s own interval excludes 0
(before v0.22 every regression band was on the logit of the survival, and on
small samples it could turn back in a tail). Here is the survival at a
covariate value with its 95% band:
x_grid = np.linspace(1, 40, 200)
band = m_cb.cb(x_grid, Z=[0.5], on='sf') # (n, 2): [lower, upper]
sf = m_cb.sf(x_grid, Z=[0.5])
plt.plot(x_grid, sf, 'b', label='S(x | Z=0.5)')
plt.fill_between(x_grid, band[:, 0], band[:, 1], alpha=0.2,
label='95% confidence band')
plt.legend()
plt.xlabel('Time')
plt.ylabel('S(x)')
plt.show()
cb takes a single covariate vector Z and accepts on='ff',
'Hf', 'hf' or 'df', one- or two-sided bounds via bound=, and
any alpha_ci. A one-sided bound is what a reliability demonstration
usually needs — for example the lower 90% bound on the reliability at time 5:
m_cb.cb([5.0], Z=[0.5], on='sf', bound='lower', alpha_ci=0.1)
array([0.65790528])
The convenience method model.plot() draws the fitted survival at the mean
covariate, with this band, against a non-parametric estimate of the pooled
data (the exponentiated Nelson-Aalen estimate) — a quick visual check, though
the pooled curve ignores the covariates. The bounds here are Wald /
delta-method bounds; the likelihood-ratio bounds available for univariate
parametric fits are not implemented for the regression models.
The other families quantify uncertainty their own way: Cox through the
information matrix (p_values, and jac as shown earlier) and the robust
sandwich; Lin-Ying through its sandwich standard_errors(); Buckley-James by
bootstrap_ci; and the frailty model (below) through standard_errors()
and param_cb.
Accelerated Life (AL)
Accelerated life testing (ALT) is a branch of reliability engineering where products are tested under elevated stress conditions — higher temperature, voltage, humidity, or load — to generate failure data faster than would be possible at normal operating conditions. The failures observed at high stress are then extrapolated back to normal conditions using a physical model for how the stress affects the life of the product.
This is fundamentally different from the regression models above. In PH, AFT, and PO, the covariates are measured characteristics of each unit (e.g. tire gauge, patient age). In AL, the covariate is a controlled experimental condition (stress level), and there are typically only two or three distinct levels. The relationship between stress and life is not statistical but physical, and the choice of life model reflects domain knowledge about the failure mechanism.
The AL model substitutes the life parameter \(\theta\) of a distribution with a stress function \(\phi(Z)\):
For example, in a Weibull AL model the scale parameter \(\alpha\) becomes \(\phi(Z)\), while the shape parameter \(\beta\) is estimated globally across all stress levels (the assumption being that the failure mechanism is the same at all stresses, just faster or slower). Which parameter is the “life” for each distribution — and how the Exponential, Log-Normal and Gamma convert between a life and their own parameter — is tabulated on the Regression Analysis page. Weibull, Exponential, Normal, Log-Normal, Gamma, Gumbel and Logistic are supported.
Available life models
The life models are in surpyval.life_models (from surpyval import
life_models); until v0.22 they were importable from surpyval itself,
where the exponential one was ExponentialLifeModel. The choice of life
model depends on the physical failure mechanism. The
letters in each formula are the parameter names the fitted model reports, and
\(Z_1, Z_2\) are the two columns of Z for the two-stress models
(\(Z_0, Z_1, \ldots\) its columns for GeneralLogLinear, numbered from
0 as its coefficients are):
Life model |
Formula \(\phi(Z)\) |
Typical use |
|---|---|---|
|
\(b \cdot e^{a/Z}\) (Arrhenius) |
Thermally-activated (chemical, diffusion, electromigration) |
|
\(Z^{-1} e^{-(b - a/Z)}\) |
Temperature, from reaction-rate (transition-state) theory: Arrhenius with a \(1/Z\) pre-factor |
|
\(1 / (a \cdot Z^n)\) |
Voltage, electrical field, mechanical fatigue |
|
\(a \cdot Z^n\) |
The same law written directly as a life; with \(n < 0\) life falls as stress rises |
|
\(a + b \cdot Z\) |
Simple first-order approximation; valid over narrow stress ranges |
|
\(c \cdot e^{a/Z_1} e^{b/Z_2}\) |
Two thermal stresses |
|
\(c \cdot Z_1^m Z_2^n\) |
Two non-thermal stresses |
|
\(c \cdot e^{a/Z_1} Z_2^n\) |
One thermal + one non-thermal |
|
\(Z e^{c - a/Z}\), the reciprocal of Eyring |
Inverse Eyring relationship |
|
\(1 / (b \cdot e^{a/Z})\), the reciprocal of Arrhenius |
Inverse Arrhenius relationship |
|
\(c \cdot e^{\beta_0 Z_0 + \beta_1 Z_1 + \cdots}\), one
|
Any number of stresses, each entering as given (pass |
A note on units: the stress variable \(Z\) for Arrhenius and Eyring should
be in Kelvin (absolute temperature), not Celsius. The accelerated life fitter
takes the same c, n, t, init and fixed arguments as the
other parametric families (the life-model parameters can be held with
fixed too), and has a fit_from_df.
Using the factory
The example simulates a classic temperature test: twenty units at each of 85 °C, 105 °C and 125 °C, with an activation energy of 0.7 eV, and a test that is stopped at 6,000 hours so that most of the coolest units are still running (right censored):
from surpyval import AcceleratedLife, Weibull, life_models
# Discrete stress levels — three temperatures in Kelvin
stress = np.repeat([358., 378., 398.], 20) # 85°C, 105°C, 125°C
Ea, k = 0.7, 8.617e-5 # activation energy eV, Boltzmann constant eV/K
rng = np.random.default_rng(42)
true_life = 1.4e-6 * np.exp(Ea / (k * stress)) # Arrhenius, in hours
T_al = true_life * rng.weibull(2.5, stress.size)
test_end = 6000.0
x_al = np.minimum(T_al, test_end)
c_al = (T_al > test_end).astype(int) # still running at the end
print('censored at each stress:',
[int(c_al[stress == s].sum()) for s in (358., 378., 398.)])
# Weibull + Arrhenius (the Exponential life model): the most common ALT
# model
model_arr = AcceleratedLife(Weibull, life_models.Exponential).fit(
x_al, Z=stress, c=c_al)
model_arr
censored at each stress: [15, 0, 0]
Parametric Regression SurPyval Model
====================================
Kind : Accelerated Life
Distribution : Weibull
Regression Model : Exponential
Fitted by : MLE
Data : 60 units: 45 events at 45 unique times, 15 right censored
Baseline : Weibull parameters; Wald 95% intervals
estimate se lower 95% upper 95%
beta 2.89 0.3316 2.308 3.619
alpha: L(Z) of the Exponential life model
Life model : Wald 95% intervals
estimate se lower 95% upper 95%
a 7815 574.1 6690 8940
b 2.79e-06 4.181e-06 -5.405e-06 1.098e-05
Notice that the Weibull shape parameter \(\beta\) is estimated globally —
it is the same for all stress levels — while the scale parameter \(\alpha\)
varies with stress via the Arrhenius relationship. This is the key assumption of
ALT: the failure mechanism does not change with stress, only the rate. The
report shows alpha as L(Z), not as a value: the life parameter
(model_arr.life_parameter) is replaced by the life model at each stress, so
it is not estimated. Its slot in params (named by model_arr.parameter_names)
holds a placeholder 1 that carries no information: it is listed in
model_arr.fixed, is not counted as a parameter in the AIC, and param_cb
refuses it. The Arrhenius parameter a is \(E_a / k_B\), so the fit
estimates the activation energy directly:
print('activation energy (eV) : %.3f' % (model_arr.params[2] * k))
print('95% CI on a, in eV :', (model_arr.param_cb('a') * k).round(3))
activation energy (eV) : 0.673
95% CI on a, in eV : [0.576 0.77 ]
# Power law — a common choice for voltage or load acceleration
model_power = AcceleratedLife(Weibull, life_models.Power).fit(
x_al, Z=stress, c=c_al
)
model_power
Parametric Regression SurPyval Model
====================================
Kind : Accelerated Life
Distribution : Weibull
Regression Model : Power
Fitted by : MLE
Data : 60 units: 45 events at 45 unique times, 15 right censored
Baseline : Weibull parameters; Wald 95% intervals
estimate se lower 95% upper 95%
beta 2.892 0.3326 2.309 3.623
alpha: L(Z) of the Power life model
Life model : Wald 95% intervals
estimate se lower 95% upper 95%
a 2.605e+56 2.306e+57 -4.26e+57 4.781e+57
n -20.56 1.488 -23.47 -17.64
Over a narrow range of temperatures a steep power law mimics Arrhenius (hence the extreme exponent), and the two fit the test data about equally well. They part company as soon as they extrapolate — here by about 20% at a use temperature only 30 °C below the coolest test — which is why the life model should come from the physics rather than from the fit statistics alone:
use = 328.0 # 55°C use condition
for name, m in [('Arrhenius', model_arr), ('Power', model_power)]:
print(f'{name:10s} AIC = {m.aic():7.2f}'
f' characteristic life at 55°C = {m.phi([use])[0]:7.0f} h')
print('true characteristic life at 55°C = %7.0f h'
% (1.4e-6 * np.exp(Ea / (k * use))))
Arrhenius AIC = 725.44 characteristic life at 55°C = 62089 h
Power AIC = 725.92 characteristic life at 55°C = 49846 h
true characteristic life at 55°C = 79830 h
Both models under-predict the true use life. The fitted activation energy,
0.67 eV against a true 0.70 eV, is well within its confidence interval, but an
error in the slope of the stress-life line is multiplied by the distance of the
extrapolation — a reminder to report the uncertainty of an extrapolated life
(for example model_arr.cb at the use stress), not just its point estimate.
The power law, whose form is wrong for this mechanism, is further off still.
To use the fitted model for extrapolation, pass the operating stress to any of the survival functions:
# Predict life at the 55°C use condition, outside the tested range
x_pred = np.linspace(0, 200000, 500)
Z_use = np.array([[use]]) # operating condition
band = model_arr.cb(x_pred, Z=Z_use, on='sf') # 95% delta-method band
plt.plot(x_pred, model_arr.sf(x_pred, Z=Z_use), label='Predicted at 55°C (328K)')
plt.fill_between(x_pred, band[:, 0], band[:, 1], alpha=0.2, label='95% band')
plt.plot(x_pred, Weibull.sf(x_pred, 1.4e-6 * np.exp(Ea / (k * use)), 2.5),
'k--', label='true reliability')
plt.xlabel('Time (hours)')
plt.ylabel('Reliability')
plt.legend()
plt.title('Extrapolated life at operating conditions')
plt.show()
The band is wide — sixty units tested for at most 6,000 hours say only so much about lives of tens of thousands of hours — and it is the band, which here contains the true curve, rather than the point estimate that should drive a decision.
Two stresses at once
The dual life models take a two-column Z, one column per stress. A
common design tests every combination of two temperatures and two voltages;
PowerExponential then combines an Arrhenius term in temperature with a
power law in voltage, \(c\, e^{a/Z_1} Z_2^{n}\):
from surpyval.life_models import PowerExponential
rng = np.random.default_rng(0)
temp = np.repeat([358., 378., 358., 378.], 25) # kelvin
volts = np.repeat([10., 10., 20., 20.], 25)
true_life2 = 2e-4 * np.exp(Ea / (k * temp)) * volts ** -1.5
x_2s = true_life2 * rng.weibull(2.5, 100)
model_2s = AcceleratedLife(Weibull, PowerExponential).fit(
x_2s, Z=np.column_stack([temp, volts]))
for name, value in zip(model_2s.parameter_names, model_2s.params):
if name != model_2s.life_parameter: # alpha is given by the life model
print(f'{name:5s} = {value:.4g}')
print('activation energy (eV): %.3f' % (model_2s.params[3] * k))
beta = 2.389
c = 0.0004321
a = 7807
n = -1.445
activation energy (eV): 0.673
The fit separates the two effects — an activation energy of 0.67 eV against
the true 0.7, and a voltage exponent n of -1.44 against the true -1.5 —
because the design varies each stress while the other is held fixed. Had voltage been raised only
together with temperature, the two columns would be collinear and no fit
could tell their effects apart. The fit then says so: where the terms the
log-life is linear in (\(1/Z\) for an exponential term, \(\log Z\)
for a power term) are collinear, or one is constant, the later stress’s
parameter is aliased – nan, with one warning naming its column – and
the others are those of the fit without it, as for a regression coefficient
the data cannot determine. Two equal stress columns of DualPower, say,
give the Power fit, with n aliased.
Creating a custom life model
If none of the built-in life models matches your failure physics, you can define
your own by subclassing LifeModel. Its constructor takes a name, a
phi_param_map naming the parameters in order, and their phi_bounds. The
two methods you must implement are:
phi(Z, *params)— the stress relationship. Useautograd.numpyso that gradients are available for the optimiser.phi_init(life, Z)— a closed-form or least-squares initialiser for the model parameters.lifeis a vector of estimated life parameters at each unique stress level;Zis the corresponding stress values. Good initialisation is important for convergence.
from surpyval import AcceleratedLife
from surpyval.life_models import LifeModel
from surpyval import Weibull
import autograd.numpy as anp
class InverseSquareRoot(LifeModel):
"""Life proportional to 1/sqrt(Z) — a simple custom example."""
def __init__(self):
super().__init__(
name="InverseSquareRoot",
phi_param_map={"a": 0},
phi_bounds=((0, None),),
)
def phi(self, Z, *params):
a = params[0]
return a / anp.sqrt(Z)
def phi_init(self, life, Z):
# life ~ a / sqrt(Z) => a ~ life * sqrt(Z)
a_est = float(anp.mean(life * anp.sqrt(Z.flatten())))
return [a_est]
model_custom = AcceleratedLife(Weibull, InverseSquareRoot()).fit(
x_al, Z=stress, c=c_al)
model_custom
Parametric Regression SurPyval Model
====================================
Kind : Accelerated Life
Distribution : Weibull
Regression Model : InverseSquareRoot
Fitted by : MLE
Data : 60 units: 45 events at 45 unique times, 15 right censored
Baseline : Weibull parameters; Wald 95% intervals
estimate se lower 95% upper 95%
beta 1.068 0.134 0.8353 1.366
alpha: L(Z) of the InverseSquareRoot life model
Life model : Wald 95% intervals
estimate se lower 95% upper 95%
a 7.357e+04 1.027e+04 5.344e+04 9.37e+04
This life model is deliberately wrong for Arrhenius data — life cannot fall steeply enough with temperature — and the fit shows it: to reconcile the three stress levels it inflates the scatter within each (the Weibull shape drops well below the true 2.5), and its AIC is far worse:
print('AIC, Arrhenius : %.1f' % model_arr.aic())
print('AIC, InverseSquareRoot: %.1f' % model_custom.aic())
AIC, Arrhenius : 725.4
AIC, InverseSquareRoot: 834.6
Time-varying covariates across families
The counting-process machinery shown for Cox is not unique to it. Wherever the
hazard at time \(t\) depends only on \(t\) and the covariate at
\(t\), the cumulative hazard is additive over disjoint time intervals and
a time-varying-covariate subject factorises exactly into one left-truncated
observation per constant-covariate interval — so the parametric proportional
hazards (PH), additive hazards (AH) and proportional odds (PO)
families fit start-stop data with the same fit_tvc / fit_tvc_timeline (and _from_df) methods and
the same i / xl / xr / c convention as Cox. (Keyword arguments
such as fixed= and init= are passed through to the ordinary fit.)
As for Cox, the _from_df methods take the covariates as Z_cols or as a
formula=, which codes categorical columns; the model keeps the formula’s
encoding, so it predicts from a DataFrame with the same design.
Fitting the truncated likelihood takes a few seconds for these 2,000 subjects,
noticeably longer than Cox:
from surpyval import WeibullPH
ph = WeibullPH.fit_tvc_from_df(
df, i_col='id', xl_col='xl', xr_col='xr', c_col='c', Z_cols='stress',
)
ph.params
array([1.97631277, 1.01510781, 1.01542508])
The data were simulated with an exponential baseline of rate 0.5 (a Weibull with scale 2 and shape 1) and \(\beta = 1\), which the fit recovers.
Accelerated failure time also fits start-stop data through the same
fit_tvc interface (and fit_tvc_timeline and the _from_df forms of
both). AFT rescales the time axis rather than the hazard, so a
subject’s likelihood depends on its accumulated accelerated age
\(\psi = \sum e^{\beta'z}\,(b - a)\) across intervals and cannot be
reshaped into independent left-truncated rows the way PH/AH/PO can; WeibullAFT
fits it with a dedicated accumulated-age likelihood instead, but the call is
identical (it accepts fixed= but not init=):
from surpyval import WeibullAFT
aft = WeibullAFT.fit_tvc_from_df(
df, i_col='id', xl_col='xl', xr_col='xr', c_col='c', Z_cols='stress',
)
aft.params
array([1.99104598, 1.00639531, 1.02502731])
The true baseline is a Weibull of shape 1, for which accelerated failure time and proportional hazards are the same model with \(\beta_{PH} = \text{shape} \times \beta_{AFT}\), so the AFT coefficient is about 1 as well.
Because the accelerated age is integrated from time zero, the AFT fit needs each subject’s whole covariate history: every subject’s first interval must start at 0 and its intervals must be contiguous. Delayed entry or a gap would mean ageing at an unknown rate over the unobserved stretch, so instead of guessing, the fit refuses and points to Cox, which handles both exactly:
late_entry = df.copy()
late_entry.loc[late_entry.index[0], 'xl'] = 0.1 # subject 0 enters at 0.1
try:
WeibullAFT.fit_tvc_from_df(late_entry, i_col='id', xl_col='xl',
xr_col='xr', c_col='c', Z_cols='stress')
except ValueError as err:
print(err)
subject 0 enters observation at 0.1 > 0: the accelerated failure time TVC likelihood integrates the covariate path from time 0, and the pre-entry covariates are unobserved. Start each subject's first interval at 0, or use Cox TVC (CoxPH.fit_tvc), which handles delayed entry exactly.
Proportional odds fits start-stop data exactly as PH and AH do: its hazard
\(h_0(t) / (F_0(t) + \phi S_0(t))\), with \(\phi = e^{\beta' z}\),
also depends only on the time and the current covariate, so a subject’s path
likelihood is the product of its per-interval, delayed-entry PO likelihoods.
Here units follow a WeibullPO model (scale 10, shape 2, and
\(\beta = 1\) on the survival odds) whose covariate switches from 0 to 1 at
a random time. Each failure time is drawn by inverting the path survival: the
unexposed curve up to the switch, then the exposed hazard from the survival
already reached.
from surpyval import WeibullPO
rng = np.random.default_rng(1)
n_po = 2000
switch = rng.uniform(2, 12, n_po)
u = rng.uniform(size=n_po)
def po_sf(t, z):
s0 = np.exp(-(t / 10) ** 2)
return np.exp(z) * s0 / (1 - s0 + np.exp(z) * s0)
def po_time(s, z):
# invert po_sf: the baseline survival at which S(t | z) = s
s0 = s / (s + np.exp(z) * (1 - s))
return 10 * np.sqrt(-np.log(s0))
T = po_time(u, 0.0) # the failure time if never exposed ...
late = T > switch # ... unless the unit outlives its switch
sw = switch[late]
T[late] = po_time(u[late] * po_sf(sw, 1.0) / po_sf(sw, 0.0), 1.0)
po_df = pd.DataFrame({
'id': np.r_[np.arange(n_po), np.flatnonzero(late)],
'xl': np.r_[np.zeros(n_po), switch[late]],
'xr': np.r_[np.where(late, switch, T), T[late]],
'c': np.r_[late.astype(int), np.zeros(late.sum(), dtype=int)],
'z': np.r_[np.zeros(n_po), np.ones(late.sum())],
})
po_tvc = WeibullPO.fit_tvc_from_df(
po_df, i_col='id', xl_col='xl', xr_col='xr', c_col='c', Z_cols='z',
)
po_tvc.params
array([10.13780126, 1.99846469, 0.9906272 ])
The fit recovers the scale, shape and coefficient. Its log-likelihood is
exactly the sum, over units, of log sf_tvc at the exit time along the unit’s
own path plus log hf at each failure, so the fitted model and the
path evaluation below are the same model.
Evaluating a covariate path. Every family that has a closed form along a
step path — Cox, parametric PH, AH and PO, and accelerated failure
time (AFT) — exposes the same sf_tvc(x, Z, xl=None, given=None) (plus
the matching Hf_tvc). Pass either (xl, Z) arrays or a
StepSchedule, and
given= for conditional survival \(S(x \mid \text{survived to } g)\).
The path is measured from time 0: a schedule that starts later has its first
value held back to 0, and the part of a schedule before 0 is ignored (the value
in force at 0 applies from there), so a constant path gives sf(x, Z)
wherever it starts. Any time is a valid query, 0 and below included.
Proportional hazards, additive hazards and proportional odds accumulate a
cumulative hazard over the segments: in each of them the hazard at time
\(t\) depends only on \(t\) and the covariate at \(t\), so
\(H(x) = \sum [H(b \mid z) - H(a \mid z)]\) over the segments
\((a, b]\). For proportional odds, with \(\phi = e^{\beta' z}\), the
constant-covariate cumulative hazard is
\(H(t \mid z) = H_0(t) - \ln\phi + \ln(F_0(t) + \phi S_0(t))\). AFT
instead accumulates an accelerated age
\(\psi(x) = \sum e^{\beta'z}\,(b - a)\) and evaluates the baseline once at
\(\psi\). A single constant segment gives sf(x, Z) in every family,
including baselines defined below zero (Normal, Gumbel, Logistic),
for which the value in force at time 0 is taken to hold before it as well.
fit_tvc treats a subject observed from time 0 the same way (its first
interval is not left-truncated), so for PH, AH and PO a constant covariate
split into intervals reproduces the ordinary fit.
An accelerated life model whose life parameter scales time (Weibull,
Exponential, Gamma, LogNormal) accumulates an age like AFT, at the rate
\(1 / L(V)\) (see Bounds, mean life and accelerated life along a path); one whose life parameter is a
location (Normal, Gumbel, Logistic) raises NotImplementedError.
Describing the covariate path
A StepSchedule is a
piecewise-constant covariate path. The covariate
must be a step function — that is precisely what keeps each family’s
cumulative form exact — so a schedule can be built structurally (a constant,
change-points, explicit intervals, or a repeating duty cycle) or from a
step-valued expression in t. Every schedule’s covariate rows have one
column per covariate, so a multivariate path is just a wider Z:
# structural
StepSchedule.constant([1.0]) # never changes
StepSchedule.from_changepoints([0.0, 500.0], [[0.0], [1.0]]) # a switch
StepSchedule.from_intervals([0, 1, 2], [1, 2, 3], [[0.], [1.], [0.]])
StepSchedule.cyclic([0, 8], [[1.0], [0.0]], period=24) # 8-on/16-off
# expression in t, materialised to a horizon
duty = StepSchedule.from_expression("1.0 if t % 24 < 8 else 0.0",
horizon=96)
trend = StepSchedule.from_expression("0.3 * 2 ** floor(t / 1000)",
horizon=5000)
trend
StepSchedule(step, 6 segment(s), p=1)
An expression is proved piecewise-constant statically, from its syntax tree,
before it is ever evaluated: t may reach the value only through a quantizer
(floor, ceil, round, trunc, //) or a comparison. A genuinely continuous covariate
(0.3 + 1e-4 * t, sin(t)) is rejected with StepValuedError rather
than silently returning a wrong answer — a covariate that varies continuously
would break the exactness of the segment sum; describe one with a
CovariatePath instead (see Continuously varying covariates). (surpyval owns only this
step-valued guarantee; sandboxing an untrusted expression string is the
calling application’s responsibility.) The expression is sampled on a grid of
spacing resolution (default 1) up to horizon, so the resolution must be
no coarser than the narrowest step; beyond the horizon the last value is held.
For several covariates pass a list of expressions, one per covariate
(StepSchedule.from_expression(["...", "..."], horizon=...)), and t0
starts the path somewhere other than 0 (a model still evaluates it from 0, as
above). The expressions may use t, numbers, arithmetic, comparisons,
and / or / not, a if cond else b, the constants pi, e,
tau and inf, and the functions floor, ceil, round,
trunc, abs, min and max, with their keyword arguments
(round(t / 10, ndigits=1)); each means what it does in Python (and and
or return an operand, so (t > 50) and 2.0 or 1.0 is 2 after
t = 50). Anything else is refused, and so is a keyword a function cannot
take.
try:
StepSchedule.from_expression("0.3 + 1e-4 * t", horizon=100)
except surv.StepValuedError as err:
print(err)
expression '0.3 + 1e-4 * t' is not step-valued: t reaches the value continuously. Quantize it (e.g. floor(t / dt)) or use a comparison so the covariate stays piecewise-constant.
Evaluating the fitted parametric model along a path is then identical to the Cox case:
t = np.linspace(0.1, 3.0, 100)
never = ph.sf_tvc(t, StepSchedule.constant([0.0]))
late = ph.sf_tvc(t, StepSchedule.from_changepoints([0.0, 1.0],
[[0.0], [1.0]]))
plt.plot(t, never, label='never stressed')
plt.plot(t, late, label='stressed after t=1')
plt.legend()
plt.xlabel('Time')
plt.ylabel('S(t)')
plt.show()
The array form gives the segment start times as xl and one covariate row per
segment, and given= conditions on survival to an age along the same path:
\(S(x) / S(g)\) after given, and 1 at and before it (survival to
those times is certain). (A univariate model’s cs(x, given) takes the
further time x instead: it is \(S(given + x)/S(given)\).)
at = np.array([0.5, 1.5, 2.5, 3.5])
pulse = dict(Z=[[0.0], [1.0], [0.0]], xl=[0.0, 1.0, 2.0]) # on for 1 < t < 2
print('S(t) :', ph.sf_tvc(at, **pulse).round(3))
print('S(t | T > 1) :', ph.sf_tvc(at, **pulse, given=1.0).round(3))
print('same as a ratio :', (ph.sf_tvc(at[1:], **pulse)
/ ph.sf_tvc([1.0], **pulse)).round(3))
S(t) : [0.781 0.3 0.114 0.068]
S(t | T > 1) : [1. 0.495 0.188 0.112]
same as a ratio : [0.495 0.188 0.112]
Continuously varying covariates
Many stress profiles are not steps. A ramp-stress test raises the load
steadily, a thermal cycle rises and falls every day, and a measured load is
sampled densely. Describe such a path with a
CovariatePath and pass it to
the same sf_tvc / Hf_tvc:
CovariatePath.from_points(times, values, period=None)draws straight lines between(time, value)points. A time given twice is a jump.CovariatePath.from_callable(func, p=1, breakpoints=None, period=None)wraps a vectorised function of time. List any kinks or jumps inbreakpointsso the integration lines up with them.
With period either one repeats. The type of Z picks the method: a
StepSchedule is still summed exactly, and along a CovariatePath the
model’s hazard is integrated,
\(H(t) = \int_0^t h\bigl(u \mid Z(u)\bigr)\, du\), by adaptive
Gauss-Kronrod quadrature to a relative error of about \(10^{-10}\) on
\(H\). If that target is missed (a path that oscillates without limit,
say), one RuntimeWarning says at how many of the query times. A path
that would need more than a million quadrature panels, such as a fast cycle
over a long horizon, raises a ValueError. Cox needs no quadrature: its
baseline hazard is a step function, so only the covariate just before each
baseline jump counts, and the result is exact.
Here the fitted ph model is evaluated along a ramp-stress profile: the
stress rises from 0 to 1 over the first time unit and then steps down to
0.5 and holds there.
from surpyval import CovariatePath
ramp = CovariatePath.from_points([0.0, 1.0, 1.0], [0.0, 1.0, 0.5])
t = np.array([0.5, 1.0, 2.0, 3.0])
print('ramp, then hold at 0.5:', ph.sf_tvc(t, ramp).round(4))
print('never stressed :', ph.sf(t, [0.0]).round(4))
print('stressed at 1 always :', ph.sf(t, [1.0]).round(4))
ramp, then hold at 0.5: [0.7237 0.4184 0.1789 0.076 ]
never stressed : [0.7805 0.606 0.3634 0.2171]
stressed at 1 always : [0.5046 0.2509 0.0612 0.0147]
Before the step the survival is that of a stress that has been climbing:
between the two constant curves, and closer to the unstressed one early on.
Approximating the ramp by a StepSchedule works too, but only in the
limit. With each step at the ramp’s midpoint value the error falls as the
square of the step width, while the path gives the limit directly:
def midpoint_steps(n_steps):
e = np.linspace(0.0, 1.0, n_steps + 1)
mid = ramp(0.5 * (e[:-1] + e[1:])).ravel()
return StepSchedule.from_changepoints(np.r_[e[:-1], 1.0],
np.r_[mid, 0.5])
for n_steps in (10, 100, 1000):
err = np.max(np.abs(ph.sf_tvc(t, midpoint_steps(n_steps))
- ph.sf_tvc(t, ramp)))
print(f'{n_steps:5d} steps: largest error in S(t) {err:.1e}')
10 steps: largest error in S(t) 1.7e-04
100 steps: largest error in S(t) 1.8e-06
1000 steps: largest error in S(t) 1.9e-08
given= conditions on survival to an age along the same path, here the
end of the ramp. The hazard is integrated from given on, so nothing is
subtracted:
print('S(t | survived the ramp):', ph.sf_tvc(t[1:], ramp, given=1.0).round(4))
S(t | survived the ramp): [1. 0.4276 0.1816]
A callable describes a smooth cycle. Over one time unit the stress rises
from 0 to 1 and falls back, and period=1 repeats it. The survival is
lower than at a constant stress of 0.5, the cycle’s average. This is
because the hazard multiplier \(e^{\beta z}\) is convex, so the hours at
high stress cost more than the hours at low stress save:
cycle = CovariatePath.from_callable(
lambda u: 0.5 - 0.5 * np.cos(2 * np.pi * u), period=1.0)
print('daily cycle :', ph.sf_tvc(t, cycle).round(4))
print('constant 0.5 :', ph.sf(t, [0.5]).round(4))
daily cycle : [0.6437 0.4115 0.1664 0.0668]
constant 0.5 : [0.6625 0.4351 0.1861 0.079 ]
What the path means depends on the family, just as for steps.
Proportional hazards, additive hazards and proportional odds take the
hazard at time \(t\) to be that of the covariate at \(t\). For
proportional odds this is
\(h_0(t) / (F_0(t) + \phi(Z(t))\, S_0(t))\), the limit of the step sum
and the model fit_tvc fits. Accelerated failure time follows Nelson’s
cumulative-exposure model: the path accumulates an accelerated age
\(\psi(t) = \int_0^t e^{\beta' Z(u)}\, du\), and
\(S(t) = S_0(\psi(t))\). A model fitted on fixed covariates and
evaluated along a path assumes its family’s time-varying form is right, and
for the same ramp the two forms give different answers unless the baseline
is exponential. An accelerated life model follows cumulative exposure as
AFT does, when its life parameter scales time (below).
Note
A CovariatePath evaluates a known, external path only. For
now it evaluates an already-fitted model along a path you supply: a
planned load, a test profile, ambient conditions. For the result to be
a survival probability, the path must not depend on the unit’s own
failure process (an external covariate). Fitting still uses steps. A
measured covariate is known only at its sample times, so fit it in
start-stop form with fit_tvc, as above; the step approximation’s
error shrinks with the square of the step width. A covariate driven by
the unit itself, such as a degradation signal read from it, needs a
joint longitudinal-survival model, which SurPyval does not provide.
Bounds, mean life and accelerated life along a path
cb_tvc(x, Z, xl=None, given=None, on='sf', alpha_ci=0.05,
bound='two-sided') puts confidence bounds on sf, ff or Hf
along a step schedule or a CovariatePath. They are the bounds of cb
carried along the path: a Wald bound on the baseline family’s
probability-plot scale (as for cb), with its standard error propagated
from the fitted covariance by the delta method, so a constant path gives
cb. Along a
CovariatePath the quadrature mesh is adapted once, at the fitted
parameters, and then held fixed while the parameters are perturbed. The
function the delta method differentiates is then smooth in the parameters,
and the cost is \(2k + 1\) passes along the path for \(k\)
parameters.
t2 = np.array([0.5, 1.0, 2.0])
print('S(t) along the ramp:', ph.sf_tvc(t2, ramp).round(4))
print('95% bounds:')
print(ph.cb_tvc(t2, ramp).round(4))
print('given the ramp survived:', ph.cb_tvc([2.0], ramp, given=1.0).round(4))
S(t) along the ramp: [0.7237 0.4184 0.1789]
95% bounds:
[[0.7064 0.7401]
[0.3996 0.437 ]
[0.165 0.1933]]
given the ramp survived: [[0.4036 0.4513]]
In a simulation of 1,000 fits each of a WeibullPH and a WeibullAFT
model (100 units, censored at a fixed time), the 95% bounds covered the true
survival in 94.6% to 96.2% of the fits. That held along a ramp, along a step
schedule, and conditional on survival to an age partway along the ramp.
mean_tvc(Z, xl=None, given=None) is the mean life along a path, and,
with given, the mean residual life of a unit that has survived to that
age along it. The integral to infinity is adaptive Gauss-Kronrod over panels
that grow geometrically, and its nodes are simply more query times of
sf_tvc. Each round of refinement is then one pass along the path, not an
integral for every node, and a mean takes a few milliseconds. A step schedule
is integrated as the matching piecewise-constant path.
print('mean life, along the ramp :', round(ph.mean_tvc(ramp), 4))
print('mean life, never stressed :',
round(ph.mean_tvc(StepSchedule.constant([0.0])), 4))
print('mean remaining life, given 1.0 :', round(ph.mean_tvc(ramp, given=1.0), 4))
mean life, along the ramp : 1.2091
mean life, never stressed : 1.9641
mean remaining life, given 1.0 : 1.1721
A path can stop units failing. If the hazard dies away, for example because
the stress falls to a level at which nothing fails, the survival levels off
above 0. A fraction of the units then never fails and the mean is infinite.
mean_tvc then returns inf with a warning that gives the survival
where the integration stopped, as a fitted univariate model’s mean() does
for a limited-failure population.
Accelerated life along a path. An accelerated life model sets a distribution’s life parameter to \(L(V)\), a function of the stress. Where that parameter scales time, \(S(t \mid V) = S_1(t / L(V))\), with \(S_1\) the distribution at unit life. This holds for the Weibull \(\alpha\), the Exponential and Gamma rates \(1 / L\) and the LogNormal’s \(e^{\mu}\). A changing stress then ages the unit by Nelson’s cumulative exposure,
the AFT form with the rate \(1 / L(V)\). That is the classical model of a step-stress or ramp-stress accelerated test. Here the Arrhenius model fitted above is evaluated along a test in which the temperature is ramped from 85 °C to 125 °C over 4,000 hours:
ramp_T = CovariatePath.from_points([0.0, 4000.0], [358.0, 398.0])
t_al = np.array([1000.0, 2000.0, 3000.0, 4000.0])
print('S(t), ramped :', model_arr.sf_tvc(t_al, ramp_T).round(3))
print('S(t), held at 85 °C :', model_arr.sf(t_al, [358.0]).round(3))
print('S(t), held at 125 °C:', model_arr.sf(t_al, [398.0]).round(3))
print('mean life, ramped : %.0f h' % model_arr.mean_tvc(ramp_T))
S(t), ramped : [0.995 0.905 0.423 0.005]
S(t), held at 85 °C : [0.998 0.984 0.951 0.891]
S(t), held at 125 °C: [0.302 0. 0. 0. ]
mean life, ramped : 2822 h
A step schedule gives a step-stress test the same way. For the location
families (Normal, Gumbel, Logistic) the life parameter \(\mu\) shifts the
distribution rather than rescaling time. A change of stress then carries no
accumulated age from one level to the next, so those models raise
NotImplementedError along a path.
A family that rescales time (AFT, and accelerated life) also accumulates the
same age over every period of a periodic path, so along one it integrates a
single period, \(\psi(t) = k\,\Psi_P + \psi(t - kP)\) with
\(k = \lfloor t / P \rfloor\). Ten million cycles then cost no more than
one. A hazard family has no such shortcut, because its baseline ages from
one period to the next, so a fast cycle over a long horizon still needs a
panel per period and raises the ValueError above.
Worked example: forecasting equipment on a duty cycle
Putting the two halves together — fitting on a time-varying covariate and then evaluating a future path — answers a question that comes up constantly in reliability but that a semi-parametric time-varying Cox fit (as in lifelines) cannot: my equipment runs on a duty cycle; how long will it last, and how much longer if I change that cycle? The fit tells you how much the load matters; the forecast turns that into a survival curve for any schedule you might run next, including ones the fleet has never yet seen.
Consider a fleet of pumps that alternate week-long high-load and low-load shifts. High load wears a pump faster. Pumps enter service in different weeks, so their cycles are offset — and that phase spread is exactly what lets the load effect be identified. We observe each pump in start-stop format, one row per week, with the load on that week as the covariate:
Note
When a time-varying covariate earns its keep. The load coefficient is
estimated from contrast within the risk set: at each failure, does the
pump that failed carry a heavier load than the others still running? If every
pump ran the identical cycle locked to one calendar — all high-load in the
same weeks — there is no such contrast, load is perfectly confounded with
time, and the coefficient is not estimable (a Cox fit goes singular; a
parametric fit returns a silently biased number). Operational spread across
the fleet — staggered commissioning, different shift patterns, idle periods —
is what makes the effect measurable. When your equipment really is all
operated the same way there is nothing for a time-varying covariate to
latch onto, and the right tool is an ordinary Weibull fit to the
failure times: it characterises reliability under that one
operating regime (it just cannot forecast a different cycle, because the
data never varied the load).
import pandas as pd
rng = np.random.default_rng(1)
n_pumps, horizon = 40, 60
alpha, shape, load_effect = 30.0, 2.0, 0.9 # ground truth for the demo
rows = []
for pump in range(n_pumps):
phase = pump % 2 # starts on a high or a low shift
energy = rng.exponential(1.0) # latent failure threshold
H = 0.0
for week in range(horizon):
load = 1.0 if (week + phase) % 2 == 0 else 0.0
dH = (((week + 1) / alpha) ** shape - (week / alpha) ** shape) \
* np.exp(load_effect * load)
failed = (H + dH) >= energy
rows.append((pump, week, week + 1, 0 if failed else 1, load))
if failed:
break
H += dH
pumps = pd.DataFrame(rows, columns=['pump', 'xl', 'xr', 'c', 'load'])
pumps.head()
| pump | xl | xr | c | load | |
|---|---|---|---|---|---|
| 0 | 0 | 0 | 1 | 1 | 1.0 |
| 1 | 0 | 1 | 2 | 1 | 0.0 |
| 2 | 0 | 2 | 3 | 1 | 1.0 |
| 3 | 0 | 3 | 4 | 1 | 0.0 |
| 4 | 0 | 4 | 5 | 1 | 1.0 |
Fitting is the ordinary fit_tvc_from_df call. A WeibullPH recovers both
the baseline wear-out (the Weibull alpha and beta) and the load effect
as a coefficient; exp(beta_load) is the hazard ratio of a high-load week
against a low-load one:
from surpyval import WeibullPH
model = WeibullPH.fit_tvc_from_df(
pumps, i_col='pump', xl_col='xl', xr_col='xr', c_col='c', Z_cols='load',
)
print('parameters :', np.round(model.params, 3))
print('load hazard ratio exp(beta) : %.2f' % np.exp(model.params[-1]))
parameters : [35.22 1.879 0.914]
load hazard ratio exp(beta) : 2.49
Now the part that motivates the whole exercise. Because the fit is fully
parametric, sf_tvc will evaluate the survival curve along any step
schedule you hand it — not just paths the pumps actually ran. That makes it a
planning tool: pose each candidate duty cycle as a
StepSchedule and
read off the survival it implies. Here we compare the current 50/50 cycle
against an eased cycle (one high week in four) and a punishing always-high
regime:
from surpyval.univariate.regression import StepSchedule
horizon_f = 60
cycles = {
'current: 1 week on / 1 week off': '1.0 if t % 2 < 1 else 0.0',
'eased: 1 week on / 3 weeks off': '1.0 if t % 4 < 1 else 0.0',
'always high load': '1.0',
}
t = np.linspace(0.5, horizon_f, 200)
for label, expr in cycles.items():
sched = StepSchedule.from_expression(expr, horizon=horizon_f)
sf = model.sf_tvc(t, sched)
median = t[np.argmax(sf < 0.5)] if np.any(sf < 0.5) else np.nan
print('%-34s median life ~ %4.1f weeks' % (label, median))
plt.plot(t, sf, label=label)
plt.legend()
plt.xlabel('Weeks in service')
plt.ylabel('S(t)')
plt.title('Forecast pump survival under alternative duty cycles')
plt.show()
current: 1 week on / 1 week off median life ~ 21.7 weeks
eased: 1 week on / 3 weeks off median life ~ 24.7 weeks
always high load median life ~ 17.8 weeks
The eased cycle buys a few weeks of median life over the current one, and the
always-high regime costs several — a quantitative answer to “should we throttle
back?” that falls straight out of the fitted model. None of these three curves
required running a pump on that schedule; they are the model’s forecast under a
covariate plan you supply. A semi-parametric time-varying Cox model has no
baseline hazard beyond the last observed event time, so it cannot produce a
survival curve into the future at all — which is why this scenario forecasting
needs the parametric fit_tvc / sf_tvc pair.
Shared-frailty models
When your data come in groups — lots from the same supplier, units at the
same site, repeated failures of one repairable machine — the members of a group
tend to fail more alike than units from different groups, because they share
something you did not measure. A shared-frailty model captures that with a
random hazard multiplier u shared within each group, on top of a
proportional-hazards baseline. Fit it with the Frailty(distribution)
factory (or a pre-built instance: WeibullFrailty, ExponentialFrailty,
LogNormalFrailty, GammaFrailty), passing a groups
label per observation (see Shared Frailty Models):
from surpyval import WeibullFrailty
rng = np.random.default_rng(1)
n_groups, per = 60, 6
rows_x, rows_c, rows_z, rows_g = [], [], [], []
for g in range(n_groups):
u = rng.gamma(1 / 0.6, 0.6) # group frailty, mean 1
for _ in range(per):
z = rng.normal()
t = 20 * (-np.log(rng.uniform()) / (np.exp(0.8 * z) * u)) ** (1 / 1.8)
rows_x.append(min(t, 40)); rows_c.append(0 if t <= 40 else 1)
rows_z.append(z); rows_g.append(g)
model = WeibullFrailty.fit(
x=np.array(rows_x), c=np.array(rows_c),
Z=np.array(rows_z).reshape(-1, 1), groups=np.array(rows_g),
)
print(model)
print("theta 95% CI:", np.round(model.param_cb("theta"), 3))
print("theta standard error: %.3f" % model.standard_errors()["theta"])
Shared-Frailty Regression SurPyval Model
========================================
Distribution : Weibull
Frailty : gamma (mean 1, variance theta); Kendall's tau = 0.142
Groups : 60 (observations 360, events 287)
Baseline : Weibull parameters; Wald 95% intervals
estimate se lower 95% upper 95%
alpha 22.87 1.336 20.4 25.64
beta 1.734 0.08827 1.57 1.916
Coefficients : exp(coef) is the hazard ratio given the frailty; Wald 95% intervals
coef exp(coef) se(coef) lower 95% upper 95% z p
beta_0 0.7793 2.18 0.07731 0.6277 0.9308 10.08 6.82e-24
Frailty variance : Wald 95% interval
estimate se lower 95% upper 95%
theta 0.331 0.09886 0.1843 0.5943
theta 95% CI: [0.184 0.594]
theta standard error: 0.099
The frailty variance theta (also model.frailty_variance) quantifies
the between-group spread. Its interval is built on the log scale, so it can
never include zero and says how precisely \(\theta\) is known rather than
whether it is positive; the evidence of real heterogeneity is an estimate
well clear of zero relative to its standard error — here about three standard
errors. (The data were simulated with
\(\theta = 0.6\). A variance of a random effect is hard to pin down — with
60 groups of six, this sample happens to land low, and the interval only just
misses the truth — while the coefficient, true value 0.8, is recovered well.)
The per-group posterior frailties — an empirical-Bayes estimate for each
observed group, shrunk toward 1 — are on model.frailties, keyed by group
label (as a string), and model.standard_errors() gives the Wald standard
errors of every parameter as a dictionary keyed by name. Every estimate is
also in one vector, model.params, in the order of model.parameter_names:
the baseline’s parameters, the coefficients, then theta.
Prediction comes in two flavours. The default is marginal (population
averaged), the right curve for a new unit from an unknown group; passing
group= conditions on an observed group’s posterior frailty, for another unit
from a group you have already seen, and frailty= conditions on any frailty
value you supply:
t = np.linspace(1, 40, 100)
z = np.array([0.0])
plt.plot(t, model.sf(t, z), 'k-', label='marginal (new group)')
g0 = model.group_labels[0]
plt.plot(t, model.sf(t, z, group=g0), 'b--',
label=f'conditional on group {g0}')
plt.plot(t, model.sf(t, z, frailty=2.0), 'r:', label='frailty = 2')
plt.legend(); plt.xlabel('Time'); plt.ylabel('S(t)')
plt.show()
The coefficient is a within-group effect: inside any one group a unit of \(z\) multiplies the hazard by \(e^{\beta}\). Across the population the frail groups fail first, so the marginal hazard ratio starts at \(e^{\beta}\) and shrinks over time:
times = np.array([1.0, 5.0, 10.0, 20.0, 40.0])
print('exp(beta) :', np.exp(model.beta).round(3))
print('marginal hazard ratio :',
(model.hf(times, [1.0]) / model.hf(times, [0.0])).round(3))
exp(beta) : [2.18]
marginal hazard ratio : [2.176 2.122 2.007 1.751 1.406]
The frailty is Gamma-distributed by default. family="lognormal" takes a
log-normal frailty instead, \(u = e^{w}\) with \(w\) normal of mean 0
and variance theta, as R’s frailtypack and coxme define it (so
theta is then the variance of \(\log u\)). It has no closed form, and
each group’s likelihood is integrated by adaptive Gauss-Hermite quadrature.
The two families put different weight in the tail of the frailty, so they can
disagree about how much of the spread is between groups; fitting both and
comparing their AIC is the usual check. frailty_variance (the variance of
the frailty scaled to mean 1) and kendall_tau (the dependence it induces
between two units of one group) are on one scale for both. On the kidney
catheter data (two infection times for each of 38 patients):
from surpyval import Frailty, Weibull
from surpyval.datasets import load_kidney
kidney = load_kidney()
kidney['female'] = (kidney['sex'] == 2).astype(float)
kidney['censored'] = 1 - kidney['status']
fits = {
family: Frailty(Weibull, family=family).fit_from_df(
kidney, x_col='time', c_col='censored', group_col='id',
Z_cols=['age', 'female'])
for family in ('gamma', 'lognormal')
}
for family, fit in fits.items():
print('%-9s theta %.3f Var(u)/E(u)^2 %.3f tau %.3f AIC %.2f '
'female %.2f' % (family, fit.theta, fit.frailty_variance,
fit.kendall_tau, fit.aic(), fit.beta[1]))
gamma theta 0.510 Var(u)/E(u)^2 0.510 tau 0.203 AIC 674.38 female -1.91
lognormal theta 0.593 Var(u)/E(u)^2 0.809 tau 0.196 AIC 676.06 female -1.63
The gamma frailty fits slightly better (its AIC is about 1.7 lower: weak evidence), with the same within-patient dependence (Kendall’s tau of about 0.2) but a larger effect of sex. Women’s lower infection rate holds under either family, so that conclusion does not depend on the choice; its size does.
fit_from_df names the columns instead (group_col for the groups, and
Z_cols or a formula for the covariates), and the fitted model then
predicts from a DataFrame:
lots = pd.DataFrame({'x': rows_x, 'c': rows_c, 'z': rows_z,
'lot': [f'L{g:02d}' for g in rows_g]})
by_lot = WeibullFrailty.fit_from_df(lots, x_col='x', group_col='lot',
Z_cols='z', c_col='c')
print(by_lot.feature_names, by_lot.group_labels[:3])
print(by_lot.sf([10.0], pd.DataFrame({'z': [0.0]})),
by_lot.sf([10.0], [0.0], group='L00'))
['z'] [np.str_('L00'), np.str_('L01'), np.str_('L02')]
[0.79510981] [0.7538654]
Omit Z entirely for a pure random-effects survival model (grouped data, no
covariates). Frailty(dist) takes any baseline distribution
(WeibullFrailty, ExponentialFrailty, LogNormalFrailty and
GammaFrailty are pre-built, each with a Gamma frailty – the name is the
baseline’s), on observed and right-censored data, and at least two groups are
required. When the data show
little between-group variation the estimate of theta goes to its boundary
at zero, and the frailty fit then coincides with the ordinary WeibullPH fit
(the same baseline, coefficients and likelihood). The frailty model has the
same neg_ll(), aic() and bic() as the parametric families, counting
theta as one more parameter, so the two fits can be compared directly —
here on grouped data with no frailty at all:
rng = np.random.default_rng(2)
z_ff = rng.normal(size=300)
x_ff = 20 * (-np.log(rng.uniform(size=300)) / np.exp(0.8 * z_ff)) ** (1 / 1.8)
g_ff = np.repeat(np.arange(30), 10) # 30 groups, but no frailty
no_frailty = WeibullFrailty.fit(x=x_ff, Z=z_ff.reshape(-1, 1), groups=g_ff)
ph_ff = WeibullPH.fit(x=x_ff, Z=z_ff.reshape(-1, 1))
print('theta : %.1g' % no_frailty.theta)
print('neg log-likelihood : frailty %.4f, PH %.4f'
% (no_frailty.neg_ll(), ph_ff.neg_ll()))
print('AIC : frailty %.2f, PH %.2f'
% (no_frailty.aic(), ph_ff.aic()))
theta : 1e-17
neg log-likelihood : frailty 1103.3934, PH 1103.3934
AIC : frailty 2214.79, PH 2212.79
The likelihoods agree and the frailty model pays 2 AIC units for its unused
theta: report the proportional-hazards model, since a variance on its
boundary has no meaningful Wald interval (param_cb('theta') is then
[0, inf]).
A Cox baseline. CoxFrailty is the same shared gamma frailty with the
baseline hazard left unspecified, as in CoxPH – the semi-parametric
member of the family, as CoxPH is of WeibullPH. For a given theta
it is fitted by EM over the frailties: each group’s posterior mean frailty
(closed form for the gamma), then a CoxPH fit with the log-frailties as
offsets and the frailty-weighted Breslow baseline. theta maximises the
profile of the integrated likelihood. This is the fit of R’s
coxph(Surv(time, status) ~ ... + frailty(id, dist = "gamma")), with
Efron’s ties by default (tie_method="breslow" for Breslow’s); on the
kidney data it gives R’s coefficients, standard errors, frailties and
I-likelihood:
from surpyval import CoxFrailty
cox_frailty = CoxFrailty.fit_from_df(
kidney, x_col='time', c_col='censored', group_col='id',
Z_cols=['age', 'female'])
print(cox_frailty)
Shared-Frailty Cox Regression SurPyval Model
============================================
Baseline : unspecified (Cox); efron ties
Frailty : gamma (mean 1, variance theta); Kendall's tau = 0.1694
Groups : 38 (observations 76, events 58)
Data : 76 units: 58 events at 50 unique times, 18 right censored
Coefficients : exp(coef) is the hazard ratio given the frailty; Wald 95% intervals
coef exp(coef) se(coef) lower 95% upper 95% z p
age 0.005222 1.005 0.01165 -0.01762 0.02806 0.4481 0.6541
female -1.583 0.2053 0.4484 -2.462 -0.7044 -3.531 0.000414
Frailty variance : profile-likelihood standard error; Wald 95% interval
estimate se lower 95% upper 95%
theta 0.4078 0.2381 0.1299 1.28
I-likelihood : -181.6386 (Cox partial likelihood -184.3446)
The model predicts as the parametric one does: the marginal curve by
default, a patient’s own with group=. The baseline (x, h0,
H0) is a step function, of a unit at Z = 0 with frailty 1. Twice the
gain of the I-likelihood over the Cox partial likelihood
(loglik_no_frailty, its value at theta = 0) tests for a frailty;
theta is on its boundary under the null, so the p-value is half the
chi-square one:
from scipy.stats import chi2
lr = 2 * (cox_frailty.loglik - cox_frailty.loglik_no_frailty)
print('LR statistic %.2f, p = %.3f' % (lr, chi2.sf(lr, 1) / 2))
woman = pd.DataFrame({'age': [45.0], 'female': [1.0]})
print(cox_frailty.sf([30, 100], woman).round(3),
cox_frailty.sf([30, 100], woman, group=21).round(3))
LR statistic 5.41, p = 0.010
[0.752 0.577] [0.968 0.936]
Model Selection
With several competing models it is useful to compare them on information
criteria. AIC penalises log-likelihood by the number of parameters (favouring
simpler models); BIC additionally penalises by sample size (favouring even
simpler models with larger datasets). Lower is better for both. In surpyval
the parameter count \(k\) is the number of estimated parameters — held
(fixed) parameters and the accelerated-life placeholder are not counted —
and the sample size of the BIC (and of AICc) is the number of observed
failures: exact, left- and interval-censored, the same rule as for every other
model (see Comparing models: information criteria), so a regression and the univariate
fit of the same data use the same sample size.
For the tires data, we can compare the three statistical regression families with a Weibull baseline, and try a second baseline for AFT and PO:
from surpyval import WeibullAFT, WeibullPH, LogNormalAFT, LogisticPO
from surpyval import PO
from surpyval import Weibull
models = {
'WeibullPH': WeibullPH.fit(x=x, Z=Z, c=c),
'WeibullAFT': WeibullAFT.fit(x=x, Z=Z, c=c),
'WeibullPO': PO(Weibull).fit(x=x, Z=Z, c=c),
'LogNormalAFT': LogNormalAFT.fit(x=x, Z=Z, c=c),
'LogisticPO': LogisticPO.fit(x=x, Z=Z, c=c),
}
for name, m in models.items():
print(f'{name:12s} AIC={m.aic():6.2f} BIC={m.bic():6.2f}')
WeibullPH AIC= -0.04 BIC= 2.35
WeibullAFT AIC= -0.04 BIC= 2.35
WeibullPO AIC= 0.83 BIC= 3.22
LogNormalAFT AIC= 3.23 BIC= 5.62
LogisticPO AIC= 0.16 BIC= 2.55
The PH and AFT rows are identical — for a Weibull baseline they are the same model (see Accelerated Failure Time (AFT)). Proportional odds comes within one AIC unit of them with either baseline, and the log-normal AFT is about three units behind, so the choice of baseline matters here as much as the choice of family. Differences of a unit or two are not meaningful — with 11 failures the data cannot separate these descriptions — so compare each family at its best baseline before ruling it out, and let the purpose and the diagnostics decide between close contenders.
A note of caution: AIC and BIC compare how well a model fits the observed data, not whether the model’s assumptions are correct. A PH model with a lower AIC than a PO model does not mean PH is the “true” model — it means PH uses its parameters more efficiently on this dataset. If the proportional hazards assumption is violated (e.g. survival curves cross), a lower-AIC PH model can still give misleading predictions. Goodness-of-fit diagnostics like Schoenfeld residuals (for PH) or log-log survival plots should accompany any model comparison. Information criteria are also only comparable between models fitted to the same data by full likelihood: a Cox model’s partial likelihood, the Lin-Ying estimator and Buckley-James have no comparable likelihood, so compare those on held-out predictions instead (next section).
Validating a survival predictor
Information criteria compare models on the data they were fit to. To judge how
well a model predicts, score it on held-out data. Two right-censored-standard
metrics live in surpyval.metrics.validation (importable from
surpyval.metrics). They take a matrix of predicted survival probabilities,
so they work for any model; survival_probability builds that matrix from
the parametric families, CoxPH, AdditiveHazards and the
surpyval.beta.ml tree and forest (see Comparison Tests and Validation Metrics).
For a model
whose sf takes one covariate vector at a time (BuckleyJames), stack
model.sf(times, Z[i]) row by row. A DataFrame of covariates is passed to
the model’s sf as it is, so a model fitted with fit_from_df (named
columns or a formula, string levels included) is scored from one; a row
with a missing covariate scores nan (see Missing values).
Both handle censoring by inverse-probability-of-censoring weighting (IPCW), so
a subject censored before the horizon does not silently bias the score. The
censoring distribution is the reverse Kaplan-Meier with an event taken to come
before a censoring at the same time, and an event at \(x_i\) is weighted by
\(1/\hat G(x_i-)\), as in R’s pec. Without ties between event and
censoring times the numbers equal scikit-survival’s; with such ties they differ
slightly, because scikit-survival weights the event by \(1/\hat G(x_i)\). A
score that would need \(\hat G\) where a separate training set’s estimate
has fallen to zero (a horizon past its censoring support) is nan.
The Brier score
BS(t)is the weighted mean squared error between the predicted survivalS(t | Z)and the survival indicator; the integrated Brier score (IBS) averages it over a time grid. Lower is better; a useful model scores below the marginal Kaplan-Meier reference.The time-dependent AUC (Uno’s cumulative/dynamic estimator) measures discrimination as a function of the horizon — the probability that a subject who has failed by
twas assigned a higher risk than one still event-free. 0.5 is chance, 1.0 is perfect.
The helper survival_probability() builds the predicted
survival matrix S(times | Z_i) from a fitted model. We fit a Cox model on a
training set and score it on an independent test set:
from surpyval import CoxPH, KaplanMeier
from surpyval.metrics import (
survival_probability, brier_score, integrated_brier_score, auc_td,
)
def make(seed, n=600):
r = np.random.default_rng(seed)
Z = r.normal(0, 1, (n, 2))
t = r.exponential(1.0 / np.exp(1.2 * Z[:, 0] - 0.8 * Z[:, 1]))
cens = r.exponential(np.median(t) * 3)
return np.minimum(t, cens), (cens < t).astype(int), Z
x_tr, c_tr, Z_tr = make(1)
x_te, c_te, Z_te = make(2)
cox = CoxPH.fit(x=x_tr, Z=Z_tr, c=c_tr)
times = np.quantile(x_te[c_te == 0], [0.25, 0.5, 0.75])
S = survival_probability(cox, Z_te, times)
_, bs = brier_score(x_te, c_te, S, times, x_train=x_tr, c_train=c_tr)
ibs = integrated_brier_score(x_te, c_te, S, times, x_train=x_tr, c_train=c_tr)
_, auc = auc_td(x_te, c_te, 1 - S, times)
print('Brier score :', np.round(bs, 3))
print('IBS :', round(ibs, 3))
print('AUC :', np.round(auc, 3))
Brier score : [0.072 0.115 0.142]
IBS : 0.114
AUC : [0.802 0.822 0.828]
The AUC around 0.8 shows the two covariates discriminate well. To see that the IBS is meaningful, compare it against the marginal Kaplan-Meier — a model that ignores the covariates entirely. The Cox model should score lower:
km = KaplanMeier.fit(x_tr, c_tr)
S_km = np.tile([km.sf([t])[0] for t in times], (len(x_te), 1))
ibs_km = integrated_brier_score(
x_te, c_te, S_km, times, x_train=x_tr, c_train=c_tr
)
print(f'IBS Cox = {ibs:.3f} marginal KM = {ibs_km:.3f}')
IBS Cox = 0.114 marginal KM = 0.161
Concordance
Harrell’s concordance index is the fraction of comparable pairs of subjects
that a risk score ranks in the right order (the one that failed first has the
higher score), with 0.5 for chance and 1 for perfect; the pair and tie rules
are on the Regression Analysis page. Two deaths at the same time are
not a pair by default (Therneau’s convention, as R’s survival and
lifelines); ties="harrell" counts them, as Harrell’s original definition
does. Every regression model has a
concordance method: with no arguments it scores the data the model was
fitted to, and given x, c and Z it scores those, such as a test
set. For any other score there is
surpyval.metrics.concordance_index(x, c, risk),
where the scores are mortality-like — higher means expected to fail
earlier. For a proportional hazards model the linear predictor
\(\beta'Z\) is exactly such a score, and it is the one CoxPH uses:
from surpyval.metrics import concordance_index
print('C, Cox on the training set: %.3f' % cox.concordance())
print('C, Cox on the test set : %.3f' % cox.concordance(x_te, c_te, Z_te))
print('C, a random score : %.3f' % concordance_index(
x_te, c_te, np.random.default_rng(0).normal(size=len(x_te))))
C, Cox on the training set: 0.807
C, Cox on the test set : 0.789
C, a random score : 0.505
For proportional odds, where a higher linear predictor means a longer life,
negate it first; for any model, the predicted failure probability
\(1 - S(t \mid Z)\) at a fixed time is also a valid risk score. The
concordance method makes that choice for each family (the linear
predictor for Cox, the frailty models and the Lin-Ying additive model; the
cumulative hazard at the median time scored for the parametric families,
which ranks exactly as the linear predictor with its sign; minus the linear
predictor for Buckley-James). The pairs are counted in
\(O(n \log n)\), so 50,000 subjects take a fraction of a second.
Concordance only measures ranking; pair it with the Brier score, which also
checks that the predicted probabilities are right.
Survival trees and random survival forests (beta)
When you do not know how the covariates act — thresholds, interactions, effects
that only appear in combination — a tree-based predictor can find the structure
itself. SurvivalTree recursively splits the data on one covariate at a time
and fits a survival model in each leaf; RandomSurvivalForest averages many
trees grown on bootstrap resamples. Both live in surpyval.beta.ml: they
are tested and usable, but their interface may still change (see
Machine Learning (beta)). The theory is in the Regression Analysis page.
The simulated data below have a risk that depends on an interaction: units
fail faster only when \(z_0 > 0.5\) and \(z_1 < 0.5\); the third
covariate is noise. A single shallow tree, allowed to consider every covariate
at each split (n_features_split='all'), finds the interaction on its own.
The tree kind couples the split rule with the leaf model:
'non-parametric' uses the log-rank statistic and Nelson-Aalen leaves for
observed, right-censored and left-truncated data, and its score form under the
pooled Turnbull estimate, with Turnbull leaves, for left- and interval-censored
and right-truncated data (with any truncation); 'weibull' (the default)
and 'exponential' use a likelihood split and parametric leaves. Every kind
accepts every kind of censoring and truncation; the parametric ones need an
optimiser at each candidate split when there is left or interval censoring or
truncation, at a much higher computational cost:
from surpyval.beta.ml import SurvivalTree, RandomSurvivalForest
def make_tree_data(seed, n=300):
r = np.random.default_rng(seed)
Z = r.uniform(0, 1, (n, 3))
risky = (Z[:, 0] > 0.5) & (Z[:, 1] < 0.5) # an interaction
t = 10 * r.weibull(1.5, n) * np.exp(-1.5 * risky / 1.5)
cens = r.uniform(5, 25, n)
return np.minimum(t, cens), (cens < t).astype(int), Z
xt_tr, ct_tr, Zt_tr = make_tree_data(1)
xt_te, ct_te, Zt_te = make_tree_data(2)
np.random.seed(0) # trees draw their candidate features at random
tree = SurvivalTree.fit(x=xt_tr, Z=Zt_tr, c=ct_tr, max_depth=2,
kind='non-parametric', n_features_split='all')
print(tree)
print('S(5), risky unit :', tree.sf([5.0], [0.9, 0.1, 0.5]).round(3))
print('S(5), ordinary unit:', tree.sf([5.0], [0.1, 0.9, 0.5]).round(3))
SurvivalTree(kind='non-parametric', selection='greedy')
|--- Z0 <= 0.589819
| |--- Z1 <= 0.0684096
| | |--- leaf: Nelson-Aalen, 10 units
| |--- Z1 > 0.0684096
| | |--- leaf: Nelson-Aalen, 170 units
|--- Z0 > 0.589819
| |--- Z1 <= 0.49846
| | |--- leaf: Nelson-Aalen, 57 units
| |--- Z1 > 0.49846
| | |--- leaf: Nelson-Aalen, 63 units
S(5), risky unit : [0.2]
S(5), ordinary unit: [0.765]
The root splits on \(z_0\) near 0.5, and the right-hand branch then splits
on \(z_1\) near 0.5 — the interaction, recovered without being specified.
Printing a tree shows each split, the left branch (<=) and then the right
(>) with its subtree indented under it, and each leaf’s model. Fitted from
arrays, the covariates are named by their column of Z (Z0, Z1,
…); a tree or forest fitted from a DataFrame keeps the column names as
feature_names and uses them instead (see below).
A forest averages many such trees, each grown on a bootstrap sample and
considering a random subset of n_features_split covariates at each split.
Its sf(x, Z), like a tree’s, returns a grid for a covariate matrix — one
row per covariate row, one column per time — unlike the element-wise
regression models, and its score(x, Z, c)
is the concordance of its mortality score (with the same ties option and
default as concordance_index). Trees and forests follow the
package’s missing-value rule: a row with a missing covariate is
dropped at fit time, with one warning giving the count, and predicts nan
(it is not sent down either branch of a split). The trees are grown one after
another; n_jobs=-1 grows them in parallel on every core, and given a
random_state gives the same forest. On held-out data it is compared with a
Cox model on the same metrics:
np.random.seed(0)
rsf = RandomSurvivalForest.fit(x=xt_tr, Z=Zt_tr, c=ct_tr, n_trees=10,
max_depth=3, n_features_split=2,
kind='non-parametric')
print('forest sf grid shape:', rsf.sf([3.0, 6.0], Zt_te[:4]).shape)
cox_t = CoxPH.fit(x=xt_tr, Z=Zt_tr, c=ct_tr)
grid = np.array([3.0, 6.0, 9.0])
scores = {}
for name, m in [('forest', rsf), ('Cox', cox_t)]:
S_m = survival_probability(m, Zt_te, grid)
ibs_m = integrated_brier_score(xt_te, ct_te, S_m, grid,
x_train=xt_tr, c_train=ct_tr)
C = rsf.score(xt_te, Zt_te, ct_te) if m is rsf else \
m.concordance(xt_te, ct_te, Zt_te)
scores[name] = ibs_m, C
print(f'{name:6s} IBS = {ibs_m:.3f} C = {C:.3f}')
forest sf grid shape: (4, 2)
forest IBS = 0.184 C = 0.680
Cox IBS = 0.197 C = 0.665
With ten shallow trees the forest already edges out a Cox model that cannot
represent the interaction; more and deeper trees usually widen the gap, at a
proportional cost in time. Setting kind='weibull' (the default) gives
parametric leaves and handles left and interval censoring and truncation. On
observed and right-censored data like these its split search costs the same
order as the log-rank’s, because each candidate child’s Weibull maximum
likelihood is found directly (the scale in closed form, the shape from the
one-dimensional profile likelihood); its leaves are Weibull fits, made when
the forest first predicts. With left or interval censoring or truncation every
candidate needs an optimiser, and it is much slower. Fitted trees and forests
serialise like every other model (next section).
A forest can also be validated without held-out data. Each tree is grown
without about a third of the rows, so every row can be scored by the trees
that never saw it. oob_log_likelihood() does this with the row’s full
likelihood — the density for an observed failure, \(S(x)\) for a
right-censored row, \(F(x)\) for a left-censored one,
\(S(x_l) - S(x_r)\) for an interval, each over the truncation
probability — and returns the mean per observation, so higher is better and
it works for every kind of censoring and truncation (the concordance needs
event times that can be ordered). A non-parametric leaf is a step function,
which puts no probability exactly at a time it did not see, so for this score
its survival curve is joined linearly between its drops and continued past
the last one with its average hazard. That makes its density a density per
unit of time, on the same scale as a parametric leaf’s, so forests of different
kind can be compared. feature_importances(random_state=...) shuffles
one covariate at a time among the out-of-bag rows and reports how much the
score drops, as a pandas.Series keyed by covariate name. With very few
trees a row can land only in leaves that give it zero probability, which makes
the score \(-\infty\); both methods then warn with the number of such rows,
and the importances are computed over the rows scored before and after each
shuffle. More trees, or kind="exponential", remove the problem. Here the forest
is fitted with fit_from_df, so the names are the DataFrame’s columns:
import pandas as pd
df_tr = pd.DataFrame(Zt_tr, columns=['z0', 'z1', 'z2'])
df_tr['time'], df_tr['censored'] = xt_tr, ct_tr
oob = {}
for depth in [0, 3]: # depth 0: every tree is one leaf
np.random.seed(0)
rsf_oob = RandomSurvivalForest.fit_from_df(
df_tr, x_col='time', c_col='censored',
Z_cols=['z0', 'z1', 'z2'], n_trees=30, max_depth=depth,
n_features_split=2, kind='non-parametric')
oob[depth] = rsf_oob.oob_log_likelihood()
print(f'max_depth={depth}: OOB log-likelihood {oob[depth]:.3f}')
importance = rsf_oob.feature_importances(random_state=1)
print(importance.round(3))
max_depth=0: OOB log-likelihood -2.644
max_depth=3: OOB log-likelihood -2.484
z0 0.083
z1 0.061
z2 -0.009
Name: importance, dtype: float64
The splits raise the out-of-bag log-likelihood above that of the pooled estimate, and the two covariates of the interaction carry the importance while the noise covariate \(z_2\) has almost none. A row that happens to be in every tree’s sample has no out-of-bag score; it is left out, with a warning giving the count. A restored forest keeps no training data, so these methods need the fitted one.
Because the score is a likelihood, it validates forests on data that concordance cannot handle. Below, units are only inspected every two time units, so every failure is interval censored (or left censored, before the first inspection, or right censored, still running at the last); units with \(z_0 > 0.5\) wear out about twice as fast. A non-parametric forest splits such data with the log-rank scores of the pooled Turnbull estimate:
r_ic = np.random.default_rng(5)
Z_ic = r_ic.uniform(0, 1, (200, 3))
T_ic = 10 * r_ic.weibull(1.5, 200) * np.where(Z_ic[:, 0] > 0.5, 0.5, 1.0)
inspections = np.arange(0.0, 22.0, 2.0)
k_ic = np.minimum(np.searchsorted(inspections, T_ic),
inspections.size - 1)
c_ic = np.where(k_ic == 1, -1, np.where(T_ic > 20, 1, 2))
x_ic = [inspections[j] if cj == -1 else 20.0 if cj == 1
else [inspections[j - 1], inspections[j]]
for j, cj in zip(k_ic, c_ic)]
oob_ic = {}
for depth in [0, 2]:
np.random.seed(0)
rsf_ic = RandomSurvivalForest.fit(
x=x_ic, Z=Z_ic, c=c_ic, n_trees=30, max_depth=depth,
n_features_split=2, kind='non-parametric')
oob_ic[depth] = rsf_ic.oob_log_likelihood()
print(f'max_depth={depth}: OOB log-likelihood {oob_ic[depth]:.3f}')
importance_ic = rsf_ic.feature_importances(random_state=1)
print(importance_ic.round(3))
max_depth=0: OOB log-likelihood -2.111
max_depth=2: OOB log-likelihood -2.086
Z0 0.042
Z1 -0.008
Z2 -0.005
Name: importance, dtype: float64
Again the splits beat the pooled Turnbull estimate out of bag, and the
importance falls on \(z_0\) alone (Z0: this forest was fitted from
arrays).
Both the tree and the forest take random_state: None (the default)
draws the bootstrap samples and the candidate covariates from NumPy’s global
generator, so np.random.seed reproduces them as above, while a seed gives
the fit a stream of its own that leaves the global one alone.
Conditional-inference trees
By default a node takes the best cut over every covariate it considers
(selection='greedy'). A continuous covariate offers a cut between every
pair of its values where a two-valued one offers one, so by chance alone its
best cut tends to look better: greedy search prefers covariates with many
values whether or not they matter, and it always finds a split to make.
selection='ctree' chooses the covariate first, by a p-value that allows
for the number of cuts each covariate had to choose from, and splits only if
the smallest p-value, multiplied by the number of covariates tested
(Bonferroni), is below alpha_split (0.05 by default); the cut on that
covariate is then chosen as usual. The theory is in
Regression Analysis. It works with every kind and every kind of
censoring.
Below, a two-valued covariate \(z_0\) shortens life by 30% and three
continuous covariates are noise. Over 40 simulated data sets, each tree
makes one split (max_depth=1); then the same again with no effect at
all:
def make_mixed_data(seed, effect, n=150):
r = np.random.default_rng(seed)
Z = np.column_stack([r.integers(0, 2, n), r.uniform(0, 1, (n, 3))])
t = 10 * r.weibull(1.5, n) * np.where(Z[:, 0] == 1, effect, 1.0)
cens = r.uniform(3, 25, n)
return np.minimum(t, cens), (cens < t).astype(int), Z
tallies = {}
for effect in [0.7, 1.0]:
tally = {'greedy': [0, 0, 0], 'ctree': [0, 0, 0]}
for seed in range(40):
xm, cm, Zm = make_mixed_data(seed, effect)
for selection in tally:
tm = SurvivalTree.fit(x=xm, Z=Zm, c=cm, max_depth=1,
kind='non-parametric',
n_features_split='all',
selection=selection)
j = getattr(tm._root, 'split_feature_index', None)
tally[selection][2 if j is None else int(j > 0)] += 1
tallies[effect] = tally
print(f'effect {effect}:')
for selection, (on_z0, on_noise, none) in tally.items():
print(f' {selection:6s} split on z0: {on_z0:2d} on noise: '
f'{on_noise:2d} no split: {none:2d}')
effect 0.7:
greedy split on z0: 22 on noise: 18 no split: 0
ctree split on z0: 22 on noise: 3 no split: 15
effect 1.0:
greedy split on z0: 1 on noise: 39 no split: 0
ctree split on z0: 0 on noise: 2 no split: 38
With the effect, greedy search splits on a noise covariate in 18 of the 40 data sets, ctree in 3; ctree declines to split in 15, where the evidence does not reach the 5% level. Without an effect, greedy search always splits (39 times on noise), while ctree leaves 38 of the 40 trees as a single leaf, close to the 95% its level promises.
A conditional-inference tree also needs no depth limit: it stops where the
data show no further effect. On the interaction data from above it grows
exactly the two splits of the interaction, and each split keeps the
adjusted p-value that chose it (p_value). In a forest, stopping early
keeps the trees from fitting noise:
ctree = SurvivalTree.fit_from_df(df_tr, x_col='time', c_col='censored',
Z_cols=['z0', 'z1', 'z2'],
kind='non-parametric',
n_features_split='all',
selection='ctree')
print(ctree)
oob_sel = {}
for selection in ['greedy', 'ctree']:
rsf_sel = RandomSurvivalForest.fit(
x=xt_tr, Z=Zt_tr, c=ct_tr, n_trees=30, n_features_split=2,
kind='non-parametric', selection=selection, random_state=0)
oob_sel[selection] = rsf_sel.oob_log_likelihood()
print(f'{selection:6s} forest: OOB log-likelihood '
f'{oob_sel[selection]:.3f}')
SurvivalTree(kind='non-parametric', selection='ctree')
|--- z0 <= 0.589819 (p = 0.000358)
| |--- leaf: Nelson-Aalen, 180 units
|--- z0 > 0.589819
| |--- z1 <= 0.49846 (p = 2.45e-06)
| | |--- leaf: Nelson-Aalen, 57 units
| |--- z1 > 0.49846
| | |--- leaf: Nelson-Aalen, 63 units
greedy forest: OOB log-likelihood -2.514
ctree forest: OOB log-likelihood -2.496
The unrestricted greedy trees grow until their leaves are too small to split, and the forest built from them scores a little lower out of bag than the one built from conditional-inference trees.
Truncated data are split the same way. Below, a failure is recorded only if it happened before the unit’s truncation time (retrospective sampling, so long lives are under-represented), and that time comes later for units with \(z_1 > 0.5\), which therefore show longer recorded lives although \(z_1\) has no effect on survival; \(z_0 > 0.5\) halves the life. A truncated unit’s score is that of its likelihood given its window, so the tree splits on \(z_0\) and not on \(z_1\), and its Turnbull leaves estimate the untruncated survival (at \(t = 3\), 0.848 and 0.628 for the true distributions):
r_rt = np.random.default_rng(0)
Z_rt = r_rt.uniform(0, 1, (400, 3))
T_rt = 10 * r_rt.weibull(1.5, 400) * np.where(Z_rt[:, 0] > 0.5, 0.5, 1.0)
tr_rt = r_rt.uniform(2, 20, 400) + 10 * (Z_rt[:, 1] > 0.5)
seen = T_rt <= tr_rt # the units we get to see
tree_rt = SurvivalTree.fit(x=T_rt[seen], Z=Z_rt[seen], tr=tr_rt[seen],
kind='non-parametric', n_features_split='all',
selection='ctree')
print(tree_rt)
s_rt = tree_rt.sf([3.0], [[0.2, 0.5, 0.5], [0.8, 0.5, 0.5]])[:, 0]
print('S(3), z0 = 0.2 and 0.8:', s_rt.round(3))
SurvivalTree(kind='non-parametric', selection='ctree')
|--- Z0 <= 0.475945 (p = 6.93e-11)
| |--- leaf: Turnbull, 140 units
|--- Z0 > 0.475945
| |--- leaf: Turnbull, 201 units
S(3), z0 = 0.2 and 0.8: [0.862 0.596]
The likelihood kinds ('weibull' and 'exponential') can also stop on
the size of the gain itself. Their split is chosen by the rise in the working
model’s maximised log-likelihood, and noise always gives some rise, so by
default a tree keeps splitting until min_leaf_samples or
min_leaf_failures stops it: what a forest of deep trees wants, but not a
tree used on its own. min_split_gain sets the least gain (in
log-likelihood units) a split must make: a number, 'aic' (the kind’s
degrees of freedom \(k\), 1 for the exponential and 2 for the Weibull: the
split must lower Akaike’s criterion) or 'bic' (\(k \log(d) / 2\), with
\(d\) the node’s failures). 'aic' is the recommended setting for a
single tree. Neither is a test – each split is the best of many cuts, so
noise clears the AIC penalty more often than once in a while – and
selection='ctree' remains the stop with a stated error rate. On the
no-effect data from above:
def n_leaves(node):
if hasattr(node, 'left_child'):
return n_leaves(node.left_child) + n_leaves(node.right_child)
return 1
leaves = {}
for gain in [0.0, 'aic', 'bic']:
leaves[gain] = [
n_leaves(SurvivalTree.fit(
x=xm, Z=Zm, c=cm, kind='exponential', n_features_split='all',
min_split_gain=gain)._root)
for xm, cm, Zm in (make_mixed_data(seed, 1.0) for seed in range(10))
]
print(f'min_split_gain={gain!r:5}: leaves {leaves[gain]}')
min_split_gain=0.0 : leaves [24, 23, 24, 22, 23, 26, 23, 24, 22, 23]
min_split_gain='aic': leaves [3, 1, 7, 8, 1, 9, 4, 15, 10, 10]
min_split_gain='bic': leaves [1, 1, 1, 1, 1, 1, 4, 1, 1, 1]
Saving and loading a fitted model
A fitted parametric regression model can be serialised to a plain dictionary or a JSON file and rebuilt later — so you can fit once and reuse the model without the training data on hand. This works for the fixed-form parametric families — Accelerated Failure Time, Proportional Hazards, Proportional Odds and (parametric) Additive Hazards — and for Accelerated Life models built on a built-in life model.
import tempfile, os
from surpyval import WeibullAFT
from surpyval.univariate.regression import ParametricRegressionModel
model = WeibullAFT.fit(x=x, Z=Z, c=c)
# to a dictionary (JSON-serialisable) ...
blob = model.to_dict()
# ... and back
restored = ParametricRegressionModel.from_dict(blob)
# the restored model predicts identically
import numpy as np
t = np.array([1.0, 5.0, 20.0])
Z_use = np.asarray(Z)[0]
print("match:", np.allclose(model.sf(t, Z_use), restored.sf(t, Z_use)))
match: True
Use to_json / from_json for a file directly:
path = os.path.join(tempfile.mkdtemp(), "aft.json")
model.to_json(path)
reloaded = ParametricRegressionModel.from_json(path)
print(reloaded)
Parametric Regression SurPyval Model
====================================
Kind : Accelerated Failure Time
Distribution : Weibull
Regression Model : Log Linear [exp(beta'Z)]
Fitted by : MLE
Data : 34 units: 11 events at 8 unique times, 23 right censored
Baseline : Weibull parameters; Wald 95% intervals
estimate se lower 95% upper 95%
alpha 0.2426 0.08144 0.1256 0.4684
beta 16.06 3.951 9.914 26.01
Coefficients : exp(coef) is the acceleration factor; Wald 95% intervals
coef exp(coef) se(coef) lower 95% upper 95% z p
beta_0 -0.5708 0.5651 0.1953 -0.9536 -0.1879 -2.922 0.003479
beta_1 -0.4981 0.6077 0.154 -0.8 -0.1962 -3.234 0.001223
beta_2 -1.713 0.1804 0.4739 -2.642 -0.784 -3.614 0.0003012
beta_3 1.145 3.142 0.3041 0.5489 1.741 3.765 0.0001665
If you don’t know (or don’t want to hard-code) which class wrote a file, the
package-level readers surpyval.from_json / surpyval.from_dict dispatch
on the serialised content itself and work for every serialisable SurPyval
model (see Saving and Loading Models):
import surpyval
print(surpyval.from_json(path))
Parametric Regression SurPyval Model
====================================
Kind : Accelerated Failure Time
Distribution : Weibull
Regression Model : Log Linear [exp(beta'Z)]
Fitted by : MLE
Data : 34 units: 11 events at 8 unique times, 23 right censored
Baseline : Weibull parameters; Wald 95% intervals
estimate se lower 95% upper 95%
alpha 0.2426 0.08144 0.1256 0.4684
beta 16.06 3.951 9.914 26.01
Coefficients : exp(coef) is the acceleration factor; Wald 95% intervals
coef exp(coef) se(coef) lower 95% upper 95% z p
beta_0 -0.5708 0.5651 0.1953 -0.9536 -0.1879 -2.922 0.003479
beta_1 -0.4981 0.6077 0.154 -0.8 -0.1962 -3.234 0.001223
beta_2 -1.713 0.1804 0.4739 -2.642 -0.784 -3.614 0.0003012
beta_3 1.145 3.142 0.3041 0.5489 1.741 3.765 0.0001665
Storing models in MongoDB
to_dict emits documents of native Python types only — string keys, lists,
floats, ints — so its output is BSON-safe and every serialisable model can be
stored in MongoDB directly. On the way back, surpyval.from_dict ignores
the _id field MongoDB adds and restores the right class from the document
itself:
collection.insert_one(model.to_dict())
doc = collection.find_one({"distribution": "Weibull"})
model = surpyval.from_dict(doc)
Every document also carries a "schema" version stamped by to_dict.
It changes only when a document’s shape changes incompatibly, so models
stored today stay recognisable to future SurPyval versions; a document
written by a newer schema than the installed SurPyval understands is
refused with a clear error rather than misread.
If the fitted model carried a computable parameter covariance, it is stored in
the dictionary, so the reloaded model can also produce confidence bounds
(cb, param_cb, standard_errors) without the original data, and it
survives repeated save/load cycles. Only the
prediction/inference state is serialised — the empirical overlay in plot
needs the fitted data, so re-fit if you need that. An Accelerated Life model is
rebuilt from its distribution and life-model names, so the built-in life models
round-trip but a custom LifeModel subclass (like InverseSquareRoot
above) cannot be rebuilt from a name and raises NotImplementedError when
serialised:
al_restored = surpyval.from_dict(model_arr.to_dict())
print('Arrhenius model restored:',
np.allclose(al_restored.sf([1e4], [use]), model_arr.sf([1e4], [use])))
try:
model_custom.to_dict()
except NotImplementedError as err:
print('custom life model:', type(err).__name__)
Arrhenius model restored: True
custom life model: NotImplementedError
A model fitted from a DataFrame keeps its covariate names, and a formula model
keeps everything its formula learned from the data – the levels of each
categorical, in order (so the reference level is unchanged), and the fitted
statistics of transforms such as scale(), poly() and bs() – so the
restored model still predicts from a DataFrame of raw covariates, exactly as
the original does.
The semi-parametric regression models save and load the same way, each on
its own result class: Cox proportional hazards
(SemiParametricRegressionModel), the Lin-Ying additive-hazards model
(AdditiveHazardsModel), the semi-parametric proportional odds model
(ProportionalOddsModel), and the Buckley-James AFT (BuckleyJamesModel).
Because their baseline is nonparametric, what is stored is the fitted
coefficients plus the baseline step arrays (or, for Buckley-James, the residual
survival), so the reloaded model predicts identically:
from surpyval import CoxPH
from surpyval.univariate.regression import SemiParametricRegressionModel
cox = CoxPH.fit(x=x, Z=Z, c=c)
cox_reloaded = SemiParametricRegressionModel.from_dict(cox.to_dict())
print("match:", np.allclose(cox.sf(t, Z_use), cox_reloaded.sf(t, Z_use)))
match: True
Cox’s time-varying-covariate prediction (predict_tvc), the additive model’s
covariance, and Buckley-James’s bootstrap_ci (which keeps a copy of the fit
data) all survive the round-trip. The residual diagnostics, check_ph and
the robust variance need the training data and are not available on a restored
Cox model, and a stratified Cox model cannot be serialised at all. Be aware that
a serialised Buckley-James model contains its training data.
The frailty model, survival trees and random survival forests serialise too —
a tree as its node structure with each leaf’s fitted model, a forest as its
trees — and all of them restore through surpyval.from_dict:
frailty_back = surpyval.from_dict(by_lot.to_dict())
forest_back = surpyval.from_dict(rsf.to_dict())
print(type(frailty_back).__name__, np.allclose(
frailty_back.sf([10.0], [0.0]), by_lot.sf([10.0], [0.0])))
print(type(forest_back).__name__, np.allclose(
forest_back.sf([5.0], Zt_te[0]), rsf.sf([5.0], Zt_te[0])))
FrailtyModel True
RandomSurvivalForest True
A restored tree or forest is a predictor only: it keeps no training data and cannot be re-fitted.