Regression Analysis
The time until some event happens will, almost certainly, be impacted by factors. For example, when considering how long a machine will last before failure an engineer will want to account for the operational conditions. It may operate in a humid environment, or it may operate at a higher rate. The question is then, how do we account for these variations, or ‘covariates’, on the time until failure?
Regression analysis is the process of capturing the effect that covariates have on the item. That is, we use data on other factors to ‘regress’ onto the survival distribution. The purpose of this type of regression is so that you can ask, and answer, questions like “what effect will increasing X have on the survival time?”
Survival regression differs from ordinary (least-squares) regression in two ways. First, the response is a time that is often not fully observed: an item still running when the test stopped is right censored, one that was only inspected occasionally is interval censored, and one that only entered the study after it had already survived a while is left truncated (see Types of Data). Second, the answer is a whole distribution of lifetimes for each covariate value, not a single predicted number. Every model on this page is therefore defined by how the covariates change a survival distribution, and every one is fitted by a likelihood that knows about censoring and truncation.
Surpyval covers five families of regression model, distinguished by how the covariates act on the distribution:
Proportional Hazards (they multiply the hazard),
Accelerated Failure Time (they scale the time axis),
Accelerated Life (they replace the life parameter with a physical stress-life relationship),
Proportional Odds (they multiply the odds of survival), and
Additive Hazards (they add to the hazard).
Around these sit a number of variations — semi-parametric versions that leave the baseline unspecified (Cox, Lin-Ying, Buckley-James), a random-effects (frailty) version for grouped data, stratified and time-varying-covariate versions, and tree-based predictors that make no structural assumption at all. There are special cases when several of these coincide, however, it is important to understand the difference between them in general. I detail the differences in the following sections; the worked, runnable versions of everything here are on the Regression Modelling with SurPyval page, and the complete API is under Regression Modelling.
Notation. Throughout, \(T\) is the (random) lifetime and \(x\) an observed time or an age at which a function is evaluated. The sections on the Cox model and later, which follow a process along the time axis, also write \(t\) for a point on that axis and \(t_k\) for the \(k\)-th distinct failure time; truncation bounds are always subscripted, \(t_l\) and \(t_r\) (the tl/tr of surpyval’s t). The covariates of one unit are a row vector \(Z = (z_1, \dots, z_p)\) and the coefficients a column vector \(\beta\), so \(\beta' Z = \beta_1 z_1 + \dots + \beta_p z_p\) is a single number, the linear predictor. The survival function is \(S(x) = P(T > x)\), the CDF \(F = 1 - S\), the density \(f\), the hazard \(h = f / S\) and the cumulative hazard \(H = -\log S\). A subscript \(0\) marks the baseline — the distribution of a unit whose covariates are all zero. Censoring follows the surpyval convention: c = 0 observed, c = 1 right censored, c = -1 left censored and c = 2 interval censored.
The table below is the one-line summary of each family — what the covariates do, and how to read a coefficient. The sign convention matters: in the PH, AFT and additive families a positive coefficient means a shorter life, while in the proportional odds family it means a longer one, and in accelerated life the stress-life function is written directly in units of life.
Family |
Definition |
Reading a coefficient |
|---|---|---|
Proportional hazards |
\(h(x \mid Z) = e^{\beta' Z} h_0(x)\) |
\(e^{\beta_j}\) is the hazard ratio per unit of \(z_j\); positive shortens life. |
Accelerated failure time |
\(S(x \mid Z) = S_0(e^{\beta' Z} x)\) |
\(e^{\beta_j}\) is the factor by which a unit of \(z_j\) speeds up ageing; positive shortens life. |
Accelerated life |
\(\text{life} = \phi(Z)\), a stress-life model |
The parameters of a physical law (activation energy, power-law exponent, …). |
Proportional odds |
\(\frac{S(x \mid Z)}{F(x \mid Z)} = e^{\beta' Z} \frac{S_0(x)}{F_0(x)}\) |
\(e^{\beta_j}\) is the survival-odds ratio; positive lengthens life. |
Additive hazards |
\(h(x \mid Z) = h_0(x) + \beta' Z\) |
\(\beta_j\) is an absolute change in hazard (a rate); positive shortens life. |
The families really are different, and the quickest way to see it is to look at the ratio of the hazard of a “treated” unit to the baseline hazard over time. The cell below takes one baseline (a log-normal) and applies a covariate effect of the same size through each of three families. Under proportional hazards the ratio is flat by definition; under proportional odds it starts at \(e^{-\beta' Z}\) and fades to 1; under accelerated failure time it changes shape with the baseline:
import numpy as np
import matplotlib.pyplot as plt
from surpyval import LogNormal
mu, sigma = 2.0, 0.5 # baseline log-normal
phi = 2.0 # the covariate effect, exp(beta'Z)
x = np.linspace(0.5, 30, 300)
S0, F0 = LogNormal.sf(x, mu, sigma), LogNormal.ff(x, mu, sigma)
h0 = LogNormal.hf(x, mu, sigma)
hr_ph = np.full_like(x, phi) # e^{b'Z} h0 / h0
hr_aft = phi * LogNormal.hf(phi * x, mu, sigma) / h0 # phi h0(phi x) / h0
hr_po = 1.0 / (F0 + (1 / phi) * S0) # odds of survival x 1/phi
plt.plot(x, hr_ph, label='proportional hazards')
plt.plot(x, hr_aft, label='accelerated failure time')
plt.plot(x, hr_po, label='proportional odds')
plt.axhline(1.0, color='grey', lw=0.5)
plt.xlabel('x'); plt.ylabel('h(x | Z) / h0(x)'); plt.legend()
plt.show()
(The proportional-odds curve uses a survival-odds multiplier of \(1/2\) so that, like the other two, the covariate is harmful.) Deciding which of these shapes matches your data is the central modelling choice; the last section, Validating a survival predictor, closes with a guide to making it.
Proportional Hazards Model
Every distribution can be described by its hazard rate \(h(x)\) — the instantaneous risk of failing at age \(x\), given survival to \(x\) — and the density, CDF and survival function all follow from it (see the Handy References - Aide-mémoire page). So a natural way to let covariates act on a lifetime is to let them act on the hazard. A proportional hazards model multiplies it:
where \(h_0\) is the baseline hazard and \(\phi(Z) > 0\) a function of the covariates. A concrete picture: pumps run at high load (\(z = 1\)) wear twice as fast as pumps at low load (\(z = 0\)), so \(\phi(1) = 2\), \(\phi(0) = 1\). At every age, a high-load pump is twice as likely to fail in the next hour as a low-load pump of the same age; the shape of the hazard over age — rising for wear-out, flat for random failures — is the same for both, only its level differs.
Which \(\phi\)? The simplest idea, \(\phi(Z) = 1 + \beta' Z\) (a linear relative risk), turns negative for a protective covariate with a large enough value, and a negative hazard is impossible, so its coefficients have to be constrained. The standard choice avoids that:
It is positive for every \(\beta\) and \(Z\); a positive coefficient raises the hazard and a negative one lowers it; and any number of covariates enter through the one linear predictor \(\beta' Z\). This log-linear form is used by the Cox model and by every pre-built parametric PH model in surpyval. Another function can be supplied when the subject matter calls for one (see the custom-covariate-function example on the Regression Modelling with SurPyval page), but the log-linear form is the one whose coefficients have the clean reading below.
What a coefficient means. Compare two units that differ by one unit in \(z_j\) and agree on everything else. The ratio of their hazards is
at every age \(x\) — that is what “proportional” means. \(e^{\beta_j}\) is the hazard ratio: \(\beta_j = 0.7\) roughly doubles the instantaneous risk of failure, \(\beta_j = -0.7\) roughly halves it. Because the hazards stay a constant multiple apart, the survival curves are powers of one another, \(S(x \mid Z) = S_0(x)^{e^{\beta' Z}}\), and so they never cross. The multiplier only has meaning relative to the baseline: the baseline is the unit with \(Z = 0\), so centring a covariate (subtracting its mean) changes the baseline but not \(\beta\).
For a parametric proportional hazards model (a known baseline such as the Weibull, with parameters \(\theta\)) surpyval estimates \(\theta\) and \(\beta\) together by maximum likelihood. Integrating the hazard gives \(H(x \mid Z) = \phi(Z) H_0(x)\), and every other function follows:
These are all the likelihood of the next section needs. (The Cox model further down is semi-parametric — its baseline is left unspecified and its coefficients are estimated by a partial likelihood instead.) Pre-built versions exist for the Exponential, Weibull, Normal, Gumbel, Logistic, Log-Normal and Gamma baselines (WeibullPH and friends), and the PH(distribution) factory builds one from other surpyval distributions; see Parametric Regression Models.
The details on fitting proportional hazards model is detailed more in the Regression Modelling with SurPyval page.
Maximum likelihood with censoring and truncation
Every parametric regression family in surpyval — proportional hazards, accelerated failure time, proportional odds, parametric additive hazards and accelerated life — is fitted by the same likelihood. The families differ only in the covariate-aware \(f(x \mid Z)\), \(S(x \mid Z)\) and \(F(x \mid Z)\) they plug in; the treatment of the data is shared, which is why all of them accept the same censoring and truncation.
Each row \(i\) contributes the probability of what was actually seen, weighted by its count \(n_i\):
observed at \(x_i\) (
c = 0): the density \(f(x_i \mid Z_i)\);right censored at \(x_i\) (
c = 1): \(S(x_i \mid Z_i)\) — it was still working;left censored at \(x_i\) (
c = -1): \(F(x_i \mid Z_i)\) — it had already failed;interval censored in \((x_{l,i}, x_{r,i}]\) (
c = 2): \(F(x_{r,i} \mid Z_i) - F(x_{l,i} \mid Z_i)\).
If the row could only have been observed inside a truncation window \((t_{l,i}, t_{r,i})\) — for example a unit that entered the study at age \(t_{l,i}\) having already survived that long — the contribution is divided by the probability of landing in the window, \(F(t_{r,i} \mid Z_i) - F(t_{l,i} \mid Z_i)\). Putting it together,
The truncation term is where a naive analysis goes wrong: ignoring delayed entry treats every late entrant as if it had been watched from birth, so it over-represents long lives. For a window bounded on one side only, surpyval evaluates \(\log S(t_l \mid Z)\) or \(\log F(t_r \mid Z)\) directly rather than as a difference of CDFs, which keeps the term finite even when the probability underflows.
The estimate \((\hat\theta, \hat\beta)\) maximises \(\log L\) numerically (the baseline parameters \(\theta\) are optimised on a transformed scale that respects their support). Parameters can be held at known values with fixed=, e.g. a Weibull shape known from experience. The information criteria reported with a fit, \(\text{AIC} = 2k - 2\log L\) and \(\text{BIC} = k \log d - 2 \log L\), use \(d\), the number of observed failures (exact, left- and interval-censored, weighted by their counts; the number of units if there is none), as the sample size of the BIC, the rule of every SurPyval model (see Comparing models: information criteria), and count as \(k\) the parameters that were estimated — a parameter held at a known value with fixed is not counted.
Covariates far from zero. A log-linear family evaluates \(e^{\beta' Z}\), which overflows for a covariate far from zero — a calendar year, a date as a day count — and used to send the optimiser, silently, to a wrong answer. fit(..., center=True) fits every family on the covariates centred at their (count-weighted) means \(\bar Z\) and reports the baseline there: the model keeps \(\bar Z\) as center, says so in its summary, and uses \(Z - \bar Z\) in every prediction, so nothing depends on where the covariates’ zero is. By default the baseline is reported at \(Z = 0\), as always. Where the baseline family is closed under the change the link makes between the two points, the fit still runs on centred covariates and moves the baseline to \(Z = 0\) exactly, with \(s = \beta'\bar Z\): for proportional hazards the cumulative hazard at \(Z = 0\) is \(e^{-s}\) times that at the means, which a Weibull or Rayleigh scale, an Exponential rate or a Gumbel location absorbs; for accelerated failure time time is rescaled by \(e^{s}\), which every scale family absorbs (the Weibull, Exponential, Gamma, Log-Normal, Log-Logistic, Rayleigh and Exponentiated-Weibull scale, and both parameters of the Normal, Logistic and Gumbel); for proportional odds the survival odds are multiplied by \(e^{-s}\), which a Log-Logistic scale or a Logistic location absorbs. For these the model is the same wherever the covariates’ zero is; if the baseline at \(Z = 0\) cannot be represented (\(\alpha e^{s/\beta}\) overflows for a Weibull PH with \(s\) in the thousands), the fit refuses and says to pass center=True. For the other pairs (a Log-Normal, Gamma, Normal or Logistic PH baseline, a PO baseline other than those two) the model with its baseline at the means is a different model from the one at \(Z = 0\), with a different maximum likelihood, so the default fit is at \(Z = 0\) on the covariates as given, as before, and is not checked: on covariates far from zero it can fail or stop at a poor answer, and center=True is the way round that. The additive hazards models are fitted as defined too: \(h_0(x) + \beta'(Z - \bar Z)\) is a different model from \(h_0(x) + \beta' Z\).
Uncertainty. The covariance of the estimates is approximated by the inverse of the observed information — the Hessian of \(-\log L\) at the optimum (the exact one the fit computes, or a numerical one where it has none), \(\widehat{\text{Cov}} = \mathcal{I}(\hat\theta, \hat\beta)^{-1}\) (fixed parameters get a zero row and column). From it:
a Wald interval on a single parameter is formed on a scale that respects its support — the natural scale for an unbounded coefficient, the log of the distance from the bound for a positive parameter such as a Weibull scale — so the interval can never leave the parameter space;
a band on a predicted curve at covariate \(Z\) uses the delta method: the gradient \(g\) of, say, \(S(x \mid Z)\) with respect to all parameters gives \(\text{se} = \sqrt{g' \,\widehat{\text{Cov}}\, g}\), and the band is built on the baseline family’s probability-plot scale — \(\log H\) for a Weibull, Exponential, Rayleigh or Gumbel baseline, the normal quantile of \(F\) for a Normal or LogNormal one, the logit of \(F\) for the rest, as for the univariate models — from the cumulative hazard, with the bands on \(S\), \(F = 1 - S\) and \(H = -\log S\) read off from it, so that \(S\) and \(F\) stay in \((0, 1)\) and \(H\) positive, and on the log scale for \(h\) and \(f\) (so they stay positive).
Both are large-sample approximations. They are least trustworthy with few failures or with a parameter near a boundary, in which case the information matrix may not be invertible and the covariance is reported as unavailable.
Semi-Parametric
Earlier pages covered ‘parametric’ and ‘non-parametric’ survival models, so what is ‘semi-parametric’? A semi-parametric model is a survival model with a non-parametric baseline and a parametric function that acts on that baseline. Recall that a proportional hazards model is
The covariate function \(\phi\) has to be parametric — it is what the coefficients describe — but nothing forces the baseline hazard \(h_0\) to be: it can be left completely unspecified and estimated non-parametrically. A parametric covariate effect on a non-parametric baseline is a ‘semi-parametric’ model.
By far the most common of any regression model of any kind (parametric, non-parametric, and semi-parametric of all the accelerated life, proportional hazard, and accelerated time) is the Cox Proportional Hazard model [Cox1972reg], it is a semi-parametric model.
The Cox model is used in a wide variety of fields. It has been used in criminology to study the recidivism of parolees, in engineering to understand the factors affecting tire reliability, and in medical science to understand factors affecting cancer and other diseases, among many many other applications. The wide use of the model shows the utility the model has and the broad applicability to solve problems.
The partial likelihood. Cox’s insight is that \(\beta\) can be estimated without ever writing down \(h_0\). Order the distinct failure times \(t_1 < t_2 < \dots\), and at each one ask: given that exactly one unit failed at \(t_k\), what is the probability that it was the one that actually did? Every unit \(j\) still at risk has hazard \(h_0(t_k) e^{\beta' Z_j}\), so that probability is
where \(Z_{(k)}\) is the covariate of the unit that failed and \(R_k\) is the risk set at \(t_k\). The unknown baseline cancels. Multiplying over failure times gives the partial likelihood,
with counts \(n_j\) as weights. Only the order of the failure times matters, which is why the Cox model cannot tell you anything about the shape of the baseline — and why it does not need to.
Risk sets, censoring and delayed entry. Censored units never appear in the numerator, but they do sit in the risk sets of every failure time up to their censoring time — that is how they contribute information. A unit that entered observation late (left truncation, tl in surpyval) is at risk only after it entered. surpyval uses the standard (entry, exit] convention: unit \(j\) is in \(R_k\) when \(t_{l,j} < t_k \le x_j\), so a unit entering exactly at a failure time is not at risk for it, and a unit is at risk at its own failure or censoring time. Right and interval truncation cannot be expressed in this forward-in-time comparison, so the Cox fitter accepts left truncation only. For the same reason it needs to know, for every unit, whether it was still at risk at each failure time, which a left- or interval-censored observation does not say: the Cox fitter is for observed (c = 0) and right-censored (c = 1) data, and refuses left- (c = -1) and interval-censored (c = 2) rows; use a parametric family for those.
A tiny example makes the formula concrete. Four units fail in turn at times 1, 2, 3 and 4; the first and third are “exposed” (\(z = 1\)). The partial log-likelihood at \(\beta = 0.3\), computed by hand, matches the value surpyval optimises:
from surpyval import CoxPH
x_toy = np.array([1.0, 2.0, 3.0, 4.0])
z_toy = np.array([1.0, 0.0, 1.0, 0.0])
b = 0.3
by_hand = 0.0
for k in range(4): # failure k: risk set is units k..3
risk_set = z_toy[k:]
by_hand += b * z_toy[k] - np.log(np.exp(b * risk_set).sum())
cox_toy = CoxPH.fit(x=x_toy, Z=z_toy.reshape(-1, 1))
print('by hand :', round(by_hand, 6))
print('surpyval (-neg_ll):', round(-cox_toy.neg_ll(np.array([b])), 6))
by hand : -3.010776
surpyval (-neg_ll): -3.010776
Tied failure times. The argument above assumes one failure at each time. With ties (times rounded to days, inspections, genuinely discrete time) there are several conventions, chosen with method=:
Breslow [Breslow1974reg] treats the \(d_k\) tied units as if each failed against the full risk set: the log term becomes \(d_k \log \sum_{j \in R_k} n_j e^{\beta' Z_j}\). Simple and fast, but it biases \(\hat\beta\) towards zero when ties are heavy.
Efron [Efron1977reg] assumes the tied failures happened in some unknown order and removes, on average, a fraction of their weight from the risk set for each successive one. Writing \(\mathcal{R}_k = \sum_{j \in R_k} n_j e^{\beta' Z_j}\) and \(\mathcal{D}_k\) for the same sum over the tied failures, the term is \(\sum_{l=0}^{d_k - 1} \log\bigl(\mathcal{R}_k - \tfrac{l}{d_k}\mathcal{D}_k\bigr)\). It is much closer to the exact answer at almost no cost, and it is the default in R and lifelines.
Exact (
'exact') sums the sequential contribution over every ordering of the tied failures — appropriate when the ties come from rounding a continuous time. surpyval evaluates the sum over orderings as the one-dimensional integral of DeLong, Guirguis and So [DeLong1994reg], \(\int_0^\infty \prod_{j \in D_k} \bigl(1 - e^{-t\, e^{\beta' Z_j} / W_k}\bigr)\, e^{-t}\, dt\) with \(W_k\) the total risk score of the units at risk that did not fail at \(t_k\) (the form SAS uses forTIES=EXACT), by quadrature, so its cost grows only linearly with the size of a tie group. Like the next method it needs integer countsn, because each row is expanded into that many tied units.Kalbfleisch-Prentice (
'kalbfleisch-prentice'or'kp') is the exact discrete-time (conditional logistic) likelihood [KalbfleischPrentice2002reg], for time that really is discrete: the denominator sums the product of the risk scores over every subset of the risk set of size \(d_k\), \(\sum_{|D| = d_k,\, D \subseteq R_k} \prod_{j \in D} e^{\beta' Z_j}\). surpyval evaluates it, its score and its information by the recursion of Gail, Lubin and Rubinstein [Gail1981reg] (the method of R’scoxph(ties="exact")), in time proportional to \(|R_k|\, \min(d_k, |R_k| - d_k)\) rather than by listing subsets, so it has no limit on the size of a tie group; it is still slower than Efron.
The first three estimate the same continuous-time hazard ratio and differ only in how well they approximate it; with light ties Breslow and Efron are close to the exact answer and Efron is the better of the two. Kalbfleisch-Prentice is a different model: in discrete time \(e^{\beta}\) is an odds ratio of the per-period failure probabilities, which is larger in magnitude than the hazard ratio when many units fail per period, so its \(\hat\beta\) is not directly comparable with the other three. With no ties all four are identical. Every CoxPH fit defaults to Efron.
Estimation and uncertainty. The score \(U(\beta) = \partial \ell / \partial \beta = \sum_k \bigl(Z_{(k)} - \bar Z_k\bigr)\), where \(\bar Z_k\) is the risk-weighted mean covariate over \(R_k\), is solved for zero by Newton-Raphson with step-halving on the partial log-likelihood, the algorithm of R’s coxph and lifelines (falling back to a root finder, and then a direct minimisation, should it fail). The observed information \(\mathcal{I}(\hat\beta) = -\partial^2 \ell / \partial\beta\,\partial\beta'\) is the risk-weighted covariance of \(Z\) summed over failures; its inverse is the covariance of \(\hat\beta\), and the reported p_values are the Wald tests \(2\bigl(1 - \Phi(|\hat\beta_j| / \text{se}_j)\bigr)\). When the design is degenerate (for example a covariate that never varies inside a risk set) the information is singular and the standard error is reported as unavailable.
The baseline and predictions. Once \(\hat\beta\) is known the baseline is recovered non-parametrically by the Breslow estimator — at each distinct time, the number of failures divided by the total risk score at risk,
As R’s coxph, lifelines and scikit-survival do, surpyval computes all of this with the covariates centred on their (count-weighted) means \(\bar Z\): every \(Z\) above is \(Z - \bar Z\). The partial likelihood depends on the covariates only through their differences within a risk set, so \(\hat\beta\) and every prediction are unchanged by the shift, but a covariate far from zero — a calendar year, a date as a day count — can no longer overflow \(e^{\beta' Z}\). The baseline this gives is that of a unit at the means; by default it is reported at \(Z = 0\), \(\hat h_0 e^{-\hat\beta'\bar Z}\) (R’s basehaz(fit, centered = FALSE)), with phi(Z) \(= e^{\hat\beta' Z}\), and predictions combine the two on the log scale. Where that baseline is not representable — covariates so far from zero that \(e^{-\hat\beta'\bar Z}\) over- or underflows — the fit refuses and says to pass center=True, which keeps the baseline at the means (model.center, R’s basehaz(fit)) and makes phi(Z) \(= e^{\hat\beta' (Z - \bar Z)}\), the hazard ratio against a unit there. The Fine-Gray fit does the same.
This is a step function that jumps only at observed failure times: it is 0 (so \(\hat S = 1\)) before the first failure, and beyond the last observed time the curve is simply held flat, so a Cox model predicts only within the range of the data and cannot be used to extrapolate. A parametric PH model is the tool for that. Note too that the Cox model’s hf is not a smooth hazard rate: it returns the jump of the step at the latest time \(\le t\) on the fitted baseline’s time grid, times \(e^{\hat\beta' Z}\). That grid holds every distinct observed time, censoring times included, and the jump at a censoring time is zero, so hf is \(\hat h_0(t_k) e^{\hat\beta' Z}\) from a failure time \(t_k\) up to the next observed time and 0 after it. After an Efron fit the baseline takes the same tie correction as the likelihood: the \(d_k\) tied failures leave the risk set a fraction at a time, so the step at \(t_k\) is \(\sum_{l=0}^{d_k - 1} 1 / \bigl(\mathcal{R}_k - \tfrac{l}{d_k}\mathcal{D}_k\bigr)\) rather than \(d_k / \mathcal{R}_k\). It is the covariate-weighted form of the Fleming-Harrington estimator (see Non-Parametric Estimation), just as the Breslow estimator is the covariate-weighted Nelson-Aalen, and it is what R’s survival uses after an Efron fit. The two baselines differ only at tied times; after any other tie method the baseline is Breslow’s.
Accelerated Failure Time
An accelerated failure time (AFT) model applies the covariate function somewhere else. Instead of multiplying the hazard, it multiplies age: a unit with covariates \(Z\) that has been running for a time \(x\) is as worn as a baseline unit that has been running for
its accelerated age. Think of a battery cycled at a higher temperature: each hour at 60 °C does the chemical damage of, say, three hours at 25 °C, so \(\phi = 3\) and every quantity of the baseline is simply read off three times further along the time axis. It is called accelerated failure time because the covariates speed up (or slow down) the clock. The survival and CDF are the baseline’s evaluated at the accelerated age; the density and hazard carry an extra \(\phi(Z)\) factor — the Jacobian of the change of variables \(x \to \phi(Z) x\):
surpyval fits these by the maximum likelihood of the previous section; pre-built versions are WeibullAFT, LogNormalAFT and friends (Exponential, Normal, Gumbel, Logistic, Log-Normal, Gamma, Weibull), and AFT(distribution) builds one for any distribution.
What a coefficient means. surpyval uses \(\phi(Z) = e^{\beta' Z}\), so a unit with covariates \(Z\) “ages” \(e^{\beta' Z}\) times faster than the baseline: it reaches at time \(x\) the state a baseline unit reaches at \(e^{\beta' Z} x\). Equivalently the lifetime itself is scaled, \(T = T_0\, e^{-\beta' Z}\), so every quantile — the median, the B10 life — is multiplied by the time ratio \(e^{-\beta' Z}\). As with PH, a positive coefficient shortens life. Taking logs,
which is a linear regression of log-time on the covariates, with an error distribution fixed by the baseline: a log-normal baseline gives normal errors (censored linear regression on \(\log T\)), a Weibull baseline gives extreme-value errors, a log-logistic one logistic errors. This is the most direct interpretation of any family — “this covariate costs you 20% of your life” — and the reason AFT is often preferred in engineering, where a stress speeding up a physical process is exactly the mechanism.
Where AFT and PH meet. For a Weibull baseline with shape \(k\), \(H_0(x) = (x / \alpha)^k\), so
which is a proportional hazards model with \(\beta_{PH} = k\,\beta_{AFT}\). The Weibull (and its special case the Exponential) is the only distribution that is both PH and AFT [Bagdonavicius]: a Weibull PH fit and a Weibull AFT fit to the same data have the same likelihood and differ only in how the coefficients are scaled. For every other baseline — log-normal included — the two families are genuinely different models, and the hazard-ratio plot at the top of this page shows how.
Accelerated Life
An accelerated life model is, in many cases, simply the inverse of an accelerated time model. However, there are some cases where they are different. Consider an accelerated failure time model with a normal baseline:
where \(\Phi\) is the CDF of the standard normal distribution. Here \(\mu\) is the expected life of a baseline unit, and the covariates rescale time. But what we usually want to know is how the covariates change the expected life itself. So instead substitute the expected life directly:
An accelerated life model is, therefore, simply a model where the life parameter of a distribution is substituted with a function of the covariates, that is, it ‘accelerates’ the expected life, as opposed to accelerating time as per an accelerated time model. This is the standard framework of accelerated life testing (ALT) [Meeker1998]: units are tested at elevated stress — temperature, voltage, humidity, load — so that they fail quickly, and a physical stress-life relationship carries the result back to use conditions.
For each of the distributions in Surpyval their life parameter that varies is as per the following table. The built-in stress-life functions \(\phi(Z)\) are written as a life (a time). For the Exponential, Gamma and Log-Normal, whose parameter is a rate or a log-location, surpyval converts the life to that parameter:
Distribution |
Life Param |
Weibull |
alpha |
Exponential |
1./lambda ( |
Normal |
mu |
LogNormal |
mu (set to \(\log \phi(Z)\), so \(\phi\) is the median life) |
Gamma |
1./beta (the rate |
Gumbel |
mu |
Logistic |
mu |
LogLogistic |
Not Avail |
ExpoWeibull |
Not Avail |
Uniform |
Not Avail |
Beta |
Not Avail |
Given the simple substitution into the life parameter, surpyval uses MLE to calculate the parameters: the remaining distribution parameters (for a Weibull, the shape) are shared across all stress levels — the assumption that the failure mechanism is the same at every stress and only its speed changes — and the life-model parameters replace the life parameter. The life parameter itself is kept in the parameter vector as a fixed placeholder (reported as 1.0), since its value now comes from \(\phi(Z)\); it is not estimated, so it is not counted in the \(k\) of the AIC and BIC, and an accelerated life model’s AIC is directly comparable with that of any other parametric family. To start the search, surpyval fits the distribution separately at each distinct stress and regresses those lives on the stress, so the data need at least two distinct stress levels.
The built-in stress-life relationships, in surpyval.life_models (all are LifeModel instances, and a custom one can be written by subclassing LifeModel). The letters are the parameter names surpyval reports; \(Z_1, Z_2\) are the two columns of a two-stress Z:
Life model |
\(\phi(Z)\) |
Typical use |
|---|---|---|
|
\(a Z^{n}\) |
Mechanical load, voltage (with \(n < 0\)) |
|
\(1 / (a Z^{n})\) |
The inverse power law of voltage endurance and fatigue |
|
\(b\, e^{a / Z}\) |
Arrhenius: temperature (in kelvin) with \(a = E_a / k_B\) |
|
\(1 / (b\, e^{a / Z})\) |
Reciprocal of the above |
|
\(\frac{1}{Z} e^{-(b - a/Z)}\) |
Temperature, with a \(1/Z\) pre-factor from reaction-rate theory |
|
\(Z\, e^{c - a/Z}\), the reciprocal of Eyring |
Reciprocal of the above |
|
\(a + b Z\) |
First-order approximation over a narrow stress range |
|
\(c\, e^{a / Z_1} e^{b / Z_2}\) |
Two thermal-type stresses (e.g. temperature and humidity) |
|
\(c\, Z_1^{m} Z_2^{n}\) |
Two non-thermal stresses |
|
\(c\, e^{a / Z_1} Z_2^{n}\) |
One thermal and one non-thermal stress |
|
\(c\, e^{\beta' Z}\), one \(\beta_j\) per column of |
Any number of stresses (transformed as needed, e.g. \(1/T\)) |
Accelerated life versus AFT. For a Weibull, substituting a log-linear life \(\alpha(Z) = e^{a + b' Z}\) gives \(S(t \mid Z) = \exp\bigl(-(t e^{-b'Z} / e^{a})^{k}\bigr)\), which is exactly an AFT model with \(\beta = -b\). So for scale-family distributions the two coincide under a log-linear link, with opposite signs: an accelerated life coefficient says how much life a unit of stress buys, an AFT coefficient how much faster it ages. They part company for location-family distributions (Normal, Gumbel, Logistic), where accelerated life shifts the location but AFT rescales time, and whenever the stress-life relationship is not log-linear — which is the point of having physically motivated life models. The distinction follows [Bagdonavicius]; see also Handy References - Aide-mémoire.
The practical caution is extrapolation. The fitted life model is used precisely outside the tested stresses, so its form is an assumption about physics that the data can only weakly check: Arrhenius and a power law can fit three test temperatures equally well and still disagree substantially at use conditions, and the disagreement grows the further the use stress is from the tested range. Compare candidate life models, and prefer the one with a mechanism behind it.
Proportional Odds
A proportional odds model acts on the odds rather than on the hazard [Bennett1983reg]. In surpyval it multiplies the odds of survival by time \(x\), \(O(x) = S(x) / F(x)\), by \(e^{\beta' Z}\):
(Writing it for the odds of failure, \(F/S\), is the same model with the sign of \(\beta\) flipped.) Solving for the survival function and differentiating,
What a coefficient means. \(e^{\beta_j}\) is the ratio of the odds of surviving past any given time, per unit of \(z_j\). Because it multiplies the odds of survival, a positive coefficient lengthens life — the opposite sign from the PH and AFT families. Keep this in mind when comparing fits: on the same data a PO model reports coefficients of the opposite sign to a PH model.
Its defining feature is that the covariate effect decays over time. The hazard ratio is
which starts at \(e^{-\beta' Z}\) when \(x\) is small (\(S_0 \approx 1\)) and tends to 1 as \(x \to \infty\) (\(S_0 \to 0\)). Two survival curves under a proportional odds model converge rather than staying a constant multiple apart. This makes it the natural choice when a treatment or covariate matters early but its influence fades, a pattern proportional hazards cannot represent.
Two baselines make it especially natural. With a log-logistic baseline the survival odds are \((x/\alpha)^{-k}\), and multiplying by \(e^{\beta' Z}\) is the same as rescaling \(\alpha\) — so a log-logistic PO model is also an AFT model, the log-logistic’s counterpart of the Weibull’s PH/AFT coincidence. With a logistic baseline the survival odds are \(e^{-(x - \mu)/\sigma}\), and the multiplier shifts \(\mu\) by \(\sigma \beta' Z\) — a location shift of the whole distribution. Pre-built versions (LogisticPO, WeibullPO, …) and the PO(distribution) factory are fitted by the same censored and truncated likelihood as the other parametric families. Proportional odds is not fitted from time-varying data, but its hazard depends only on the time and the current covariate, so a fitted PO model is evaluated exactly along a step covariate path with sf_tvc / Hf_tvc.
Additive Hazards
Where proportional hazards multiplies the baseline hazard, an additive hazards model adds to it:
The covariate shifts the absolute hazard by a constant amount at every age, rather than scaling it. This is often the more natural scale for risk-difference questions (excess deaths per unit time attributable to an exposure), and for reliability settings where hazards from separate mechanisms genuinely add — a component exposed to an extra, age-independent shock process of rate \(\beta' Z\) on top of its own wear-out. Integrating, the cumulative hazard is \(H(x \mid Z) = H_0(x) + x\, \beta' Z\), so \(S(x \mid Z) = S_0(x)\, e^{-x \beta' Z}\): the covariate multiplies survival by an exponential factor. A coefficient is a rate: \(\beta_j = 0.01\) per hour means one extra failure per hundred unit-hours for each unit of \(z_j\), whatever the baseline is doing.
The Lin-Ying estimator. Like Cox, the Lin-Ying form [LinYing1994reg] leaves the baseline hazard unspecified, but unlike Cox it admits a closed-form estimator for \(\beta\) — no iteration and nothing to converge. With \(Y_i(t)\) the at-risk indicator, \(N_i(t)\) the failure counting process and \(\bar Z(t)\) the mean covariate among those at risk,
\(A\) is the spread of the covariates in the risk set integrated over time — how long each covariate configuration was exposed — and \(b\) compares the covariates of those who failed with those at risk. The variance is the Lin-Ying sandwich \(A^{-1} B A^{-1}\), with \(B = \sum_i \int (Z_i - \bar Z)^{\otimes 2} dN_i\), from which the standard errors and Wald p_values follow. The baseline cumulative hazard is a Breslow-type step sum corrected for the covariate drift, \(\hat H_0(t) = \sum_{t_k \le t} d_k / |R_k| - \int_0^t \bar Z(s)' \hat\beta\, ds\); the drift integral is evaluated continuously (\(\bar Z\) is constant between the times at which the risk set changes), so a prediction \(H(t \mid Z) = \hat H_0(t) + t\, \hat\beta' Z\) is the same however the covariates are centred. Because a step function has no rate, the hazard \(h(t \mid Z)\) is reported with a kernel-smoothed baseline (an Epanechnikov kernel over the increments of \(\hat H_0\), bandwidth by a normal-reference rule unless you pass bandwidth=), which is least accurate near the ends of the observed time range. The semi-parametric fitter handles observed and right-censored data.
The parametric version. AH(distribution) and the pre-built WeibullAH, ExponentialAH, … use a parametric baseline, \(h(x \mid Z) = h_0(x; \theta) + \beta' Z\), and are fitted by the censored and truncated likelihood of the proportional hazards section. They give a smooth, extrapolatable version of the same model.
The positivity caveat. Nothing constrains \(h_0(x) + \beta' Z\) to be positive. For a strongly protective covariate the additive hazard can go negative — impossible for a real hazard. The semi-parametric AdditiveHazards predicts with the running maximum of its cumulative-hazard estimate, so its survival stays in \([0, 1]\) and never rises; its fitted \(\hat H_0\) is left as estimated. The parametric models predict with the model as fitted: where \(h_0(x) + \beta' Z\) is negative the cumulative hazard falls, survival exceeds 1 and the density is negative, and each prediction there (sf, Hf, cb, sf_tvc, …) warns once that the values are the model’s, not a distribution’s. The parametric likelihood needs \(\log h\) at every failure, so surpyval’s optimiser treats any parameter value that makes the hazard non-positive at an observed failure as infeasible and stays inside the feasible region. When the data would prefer a negative hazard, the fit therefore returns the best feasible model — pressed against the boundary, with a near-zero hazard at the earliest failures of the protected units and a baseline distorted to compensate — and warns that it has done so; it raises an error only if the optimiser ends at an infeasible point. When effects are strongly protective, the exponential link of proportional hazards, which keeps the hazard positive by construction, is the safer choice.
Semi-Parametric — Buckley-James
Cox leaves the baseline hazard unspecified; Buckley-James [BuckleyJames1979reg] is the accelerated-failure-time counterpart that leaves the error distribution unspecified. It fits an AFT model by iterating between imputing the censored failure times from the current fit (using the Kaplan-Meier residual distribution) and re-estimating the coefficients by least squares on the completed data. The result is a semi-parametric AFT: covariate effects on the log-time scale without committing to a parametric family for the baseline.
In detail, write the model as \(\log T = \gamma' Z + \varepsilon\) with \(\varepsilon\) from an arbitrary distribution. If there were no censoring, least squares of \(y_i = \log x_i\) on \(Z_i\) would be the obvious estimator. A censored \(y_i\) is only a lower bound, so Buckley-James replaces it by its conditional expectation given that it is at least that large:
where \(e_i = y_i - \gamma' Z_i\) are the residuals, \(\hat S\) is their Kaplan-Meier survival function and \(\Delta\hat F(e_k)\) its jumps. The two steps — impute, then refit \(\gamma\) by least squares on the \(\hat y_i\) — alternate until \(\gamma\) stops moving. Some details matter in practice:
the largest residual is treated as a failure (Efron’s tail correction), so the residual distribution is proper and every conditional mean is finite;
the intercept is not identified separately from the error distribution, so only the slopes are estimated and the location lives in the residual distribution;
the iteration can settle into a two-point cycle rather than a fixed point, a known property of the estimator; surpyval detects the cycle and averages it, and warns if the iteration has not converged;
there is no likelihood, so there are no model-based standard errors: uncertainty comes from a percentile bootstrap.
surpyval reports the coefficients as \(\beta = -\gamma\), in the same sign convention as WeibullAFT (positive shortens life), and predicts with the residual distribution directly, \(S(t \mid Z) = \hat S_\varepsilon(\log t + \beta' Z)\) — a step function, like a Kaplan-Meier. The fitter accepts observed and right-censored data with positive times.
Checking a proportional hazards fit
Every proportional hazards model rests on one assumption: that a covariate multiplies the baseline hazard by a constant factor for all time. If that is false — a treatment that helps early but not late, a covariate whose effect drifts — the single coefficient the model reports is a time-average that can be misleading.
The assumption is checked with the Schoenfeld residuals. At each event time the Schoenfeld residual for a covariate is the observed covariate value of the subject who failed minus the risk-weighted mean covariate value over everyone still at risk. If proportional hazards holds, these residuals have no trend in time; if the effect is drifting, they trend. The Grambsch-Therneau test [GrambschTherneau1994reg] formalises this by regressing the scaled Schoenfeld residuals on a transform of time and testing for a non-zero slope, both per covariate and jointly. A small \(p\)-value is evidence against proportional hazards.
The precise definitions. For the unit failing at \(t_k\) the Schoenfeld residual is
(for an Efron fit the mean is averaged over the Efron-reduced risk sets of the tie group, so the residuals are consistent with the likelihood that was maximised). Grambsch and Therneau showed that, if the coefficient is really a function of time \(\beta_j(t)\), the scaled residual \(r^*_k = \hat\beta + d\, V r_k\) — with \(d\) the number of failures and \(V = \mathcal{I}^{-1}\) the coefficient covariance — has expectation approximately \(\beta(t_k)\). Plotting \(r^*_k\) against time is therefore a picture of how the coefficient changes, centred on the fitted constant. For a time transform \(g\) (transform= in check_ph: "km", the default, uses \(g(t) = 1 - \hat S_{KM}(t)\) from a Kaplan-Meier fitted to all the data; "rank" uses average ranks; "identity" and "log" use \(t\) and \(\log t\)), let \(u = \sum_k (g_k - \bar g)\, r_k\) and \(s = \sum_k (g_k - \bar g)^2\). Then
the global and per-covariate statistics of R’s cox.zph and lifelines, which surpyval reproduces. The "km" transform is the usual choice: it spreads the failures evenly and is not dominated by a few long times. The per-covariate tests are screens; with several covariates, look at the global test and at the plots.
When the test rejects, the usual remedies are to stratify on the offending covariate (if it is a nuisance), to let its effect change over time by including an interaction with a function of time as a time-varying covariate, or to move to a family whose effect is not constant in time (AFT, proportional odds).
SurPyval exposes several other residuals for a fitted Cox model, each answering a different question: martingale residuals (observed minus expected events, \(M_i = \delta_i - e^{\hat\beta' Z_i}\bigl(\hat H_0(x_i) - \hat H_0(t_{l,i})\bigr)\), where \(\delta_i = 1\) for a failure; an Efron fit credits a tied failure with only its share of the baseline step at its own time) reveal non-linear covariate functional form when plotted against a covariate; deviance residuals, a symmetrised transform of the martingale residuals, highlight poorly-predicted individuals; score residuals are each observation’s contribution to \(U(\hat\beta)\), and dfbeta residuals (score residuals times \(V\)) approximate how much each observation moves \(\hat\beta\). The Schoenfeld, score and martingale residuals all sum to zero at the maximum of the partial likelihood. All of them respect delayed entry. They follow the tie method of a Breslow or Efron fit; after an 'exact' or 'kalbfleisch-prentice' fit they are computed with the Breslow forms, so under heavy ties they will not sum exactly to zero there. See [TherneauGrambsch2000reg] for a thorough treatment.
Cluster-robust standard errors
The model-based standard errors assume every observation is independent. When the data are clustered — repeated events on the same subject, several failures from one machine, grouped sampling — that assumption is wrong and the naive errors are too small. The Lin-Wei sandwich [LinWei1989reg] (or “robust”) variance corrects for it. Writing \(H\) for the information matrix and \(s_c\) for the sum of a cluster’s score contributions, the robust covariance is
which reduces to the usual variance when there is one observation per cluster and there is no within-cluster correlation. Since the dfbeta residual of an observation is its score residual times \(H^{-1}\), this is the same as summing the dfbeta residuals within each cluster, \(D_c = H^{-1} s_c\), and forming \(\sum_c D_c D_c^{\top}\) — which is how surpyval computes it.
The intuition for why clustering matters: if every observation were accidentally entered twice, a naive analysis would think it had twice the data and shrink the standard errors by \(\sqrt 2\); declaring each pair a cluster tells the sandwich that the copies carry no new information, and the error returns to its correct size. For a start-stop (time-varying-covariate) Cox fit the rows of one subject are correlated by construction, so the robust variance clusters by subject unless told otherwise. Counts n are frequency weights here as everywhere in surpyval: a row with n = 3 is three independent observations (three clusters, unless cluster= groups them), so the robust errors — like the fit, the check_ph ranks and the Buckley-James bootstrap — are exactly those of the data written out one row per observation.
Shared frailty
Cluster-robust errors correct for within-cluster correlation but do not model it. A shared-frailty model does the opposite: it introduces the correlation explicitly through an unobserved random effect [Hougaard2000reg]. Each group \(g\) (a manufacturing lot, a site, a repairable unit) is given a frailty \(u_g\) — a random multiplier shared by every member of the group — acting on the hazard:
with the frailties drawn once per group from a Gamma distribution of mean 1 and variance \(\theta\). The frailty is the survival analogue of a random intercept: it absorbs whatever unmeasured feature makes a whole group fail faster or slower than its covariates predict, and \(\theta\) measures that between-group variability (\(\theta = 0\) recovers ordinary proportional hazards). This is the conditional / random-effects counterpart of the marginal cluster-robust correction above — same within-group correlation, modelled rather than merely accounted for.
Because the frailty multiplies the cumulative hazard, a Gamma frailty integrates out of a group’s likelihood in closed form. Writing \(D_g\) for the number of events in group \(g\) and \(H_g = \sum_{j \in g} e^{\beta' z_j} H_0(t_j)\) for the sum of its members’ cumulative hazards, the group contributes
to the marginal log-likelihood, which is maximised jointly over the baseline parameters, \(\beta\), and \(\theta\). The same conjugacy makes the posterior frailty of an observed group a closed form, \(\hat u_g = (D_g + 1/\theta)/(H_g + 1/\theta)\) — an empirical-Bayes estimate shrunk toward 1, larger for groups that fail early. Standard errors come from the Hessian of the marginal likelihood (exact; numerical where the variance sits at its limit of 0 and the exact one is singular); the Wald interval for \(\theta\) (and for the positive baseline parameters) is formed on the log scale so it stays positive. The closed form requires observed and right-censored data only. Gamma is the only frailty distribution surpyval offers (other choices, such as a log-normal frailty, have no closed form and need numerical integration); the baseline can be any surpyval distribution, via Frailty(distribution).
The distinction between the two curves the model can draw matters. Integrating the frailty out gives the marginal (population-averaged) survival of a unit from an unknown group, \(S(t \mid Z) = (1 + \theta \, e^{\beta' Z} H_0(t))^{-1/\theta}\) — a Laplace transform of the frailty distribution, and always heavier-tailed than the baseline. Conditioning on a value \(u\) gives \(S(t \mid Z, u) = e^{-u \, e^{\beta' Z} H_0(t)}\), used with \(\hat u_g\) to predict a new member of an already-observed group. That is a plug-in of the posterior mean: the group’s frailty is still uncertain after its data (its posterior is a Gamma with shape \(D_g + 1/\theta\) and rate \(H_g + 1/\theta\)), and averaging over that posterior instead gives \(\bigl(1 + \theta\, e^{\beta' Z} H_0(t) / (1 + \theta H_g)\bigr)^{-(D_g + 1/\theta)}\), which has a slightly heavier tail. The two agree closely for a group with many events, whose frailty is well determined. A subtle consequence is that a mixture of groups makes the population hazard bend down over time even when every group’s hazard rises, because the frail groups fail first and leave robust survivors — so an apparent decreasing hazard can be a heterogeneity artifact rather than a real one. Identification requires within-group replication: with a single group, or one observation per group, \(\theta\) is confounded with the baseline shape and cannot be estimated.
The same selection effect changes what the coefficients mean. The marginal hazard is
so the population hazard ratio between two covariate values starts at \(e^{\beta}\) and shrinks towards 1 over time, even though within every group it is exactly \(e^{\beta}\). \(\beta\) in a frailty model is a within-group (conditional) effect, and it is typically larger in magnitude than the coefficient an ordinary PH fit reports on the same data. The cell below shows the bending for a baseline whose hazard increases, with \(\theta = 1\):
from surpyval import Weibull
t = np.linspace(0.01, 30, 300)
h0, H0 = Weibull.hf(t, 10, 1.5), Weibull.Hf(t, 10, 1.5)
for theta in [0.0, 0.5, 1.0]:
plt.plot(t, h0 / (1 + theta * H0), label=f'theta = {theta:g}')
plt.xlabel('t'); plt.ylabel('population (marginal) hazard'); plt.legend()
plt.show()
Every group’s hazard is the rising \(\theta = 0\) curve times its own \(u_g\), yet the population hazard for \(\theta = 1\) rises and then falls. As \(\theta \to 0\) the marginal cumulative hazard \(\log(1 + \theta e^{\beta' Z} H_0)/\theta\) tends continuously to the proportional-hazards one, \(e^{\beta' Z} H_0(t)\), and surpyval’s marginal predictions switch to that limit when \(\theta\) is numerically zero. Because \(\theta\) cannot be negative, data with no between-group heterogeneity push the estimate onto that boundary, where the Wald interval is no longer meaningful (the log-scale interval degenerates to \([0, \infty)\) or to a single point). A \(\hat\theta\) at or very near zero says the grouping explains nothing beyond the covariates, and the frailty fit then reproduces the ordinary proportional-hazards fit — same baseline, coefficients and likelihood: fit and report that model instead.
Stratification
When proportional hazards fails for a nuisance covariate — a study site, a batch, a device generation you would rather not model — the standard remedy is stratification: allow a separate baseline hazard \(h_{0,g}(t)\) for each stratum \(g\) while sharing the coefficients \(\beta\). Because the partial likelihood is summed within strata, risk sets never cross a stratum boundary and the nuisance factor is removed from the comparison without ever estimating its effect. The Cox partial likelihood factorises across strata, so this is a small change to the estimation with a large gain in robustness:
each \(\ell_g\) being the ordinary partial likelihood of stratum \(g\) alone. There is one Breslow baseline per stratum, so a prediction must say which stratum it is for. The price is that the stratifying variable gets no coefficient — you learn nothing about its effect, which is the point when it is a nuisance and a loss when it is not. The residual diagnostics and robust variance above assume a single baseline, so they are not available on a stratified fit.
Time-varying covariates
A covariate can change during a subject’s follow-up — a dose is raised, a treatment begins, a machine is moved to a harsher environment. Such a covariate path is represented in the counting-process (start-stop) format: each subject contributes one row per interval \((x_l, x_r]\) on which its covariate vector is constant, and only the interval ending at the subject’s event carries the terminal status. surpyval writes this with its usual vocabulary — the subject id i, the interval bounds xl / xr, and the censoring flag c (0 event, 1 right-censored interval end). A subject’s rows may not overlap, it may have at most one event, and that event must be on its last interval; gaps between intervals are allowed for the hazard-based families (the subject is simply not at risk in the gap). Equivalently the path can be given as a timeline — one row per covariate change, each value holding until the next row, with the first time the entry and the last row carrying the exit time and status — which surpyval expands into the same intervals.
Fitting. Whether a subject’s episodes can be fitted as if they were independent observations depends on the family. Where the cumulative hazard is additive over disjoint intervals — proportional hazards, additive hazards, and (semi-parametric) Cox — a subject splits exactly into one left-truncated (delayed-entry) observation per constant-covariate interval:
is the interval’s contribution and the subject’s likelihood is the product over its intervals, so the ordinary maximum-likelihood fitter recovers the same estimate from the reshaped rows (the episode-splitting identity). Accelerated failure time is different: the covariate rescales the time axis, so the baseline is evaluated at the subject’s accumulated accelerated age \(\psi = \sum_k e^{\beta' z_k}\,(b_k - a_k)\), and the episode entry ages in that sum depend on \(\beta\). The episodes therefore do not separate into independent rows, and AFT is fitted with a dedicated accumulated-age likelihood that re-accumulates \(\psi\) per subject on each optimiser step. Proportional odds is not fitted from time-varying data.
For AFT the subject’s contribution is
with \(\delta = 1\) for a subject whose last interval ends in failure. The integral runs from time zero, which has a consequence: the accelerated age accumulated before a subject entered observation, or during a gap in its record, depends on covariate values that were never observed. Rather than guess them, surpyval’s AFT time-varying fit requires each subject’s intervals to start at 0 and to be contiguous, and refuses delayed entry and gaps with an error that points to the Cox time-varying fit, whose risk-set likelihood never needs the unobserved history. Information criteria for this fit count subjects, not interval rows.
Evaluation. Given an already-fitted model and a covariate path \(Z(t)\), the survival \(S(t \mid Z(\cdot))\) is exact when the path is piecewise-constant (a step function). For proportional and additive hazards (and Cox) it is the sum of the per-segment cumulative-hazard increments; for AFT it is the baseline evaluated at the accumulated accelerated age. A continuously-varying covariate would break the exactness of these segment sums, so the step-valued requirement is a correctness precondition, not a convenience. Unlike some packages, surpyval is willing to evaluate a model along a future covariate path — because that path is supplied as a plan or hypothesis (mission phases, a duty cycle, a scheduled load), making \(S(t \mid Z(\cdot))\) a well-posed conditional question rather than a claim to know the future.
Two details of the evaluation are worth knowing. For Cox, the baseline hazard jumps at observed event times, and a jump that falls exactly on a covariate change time is weighted by the old covariate — the same (x_l, x_r] convention as the fit, where a unit is still at risk at the end of its interval. And conditional survival — the probability of surviving to \(x\) for a unit already known to have survived to age \(g\) along the same path — is
Identifiability. A time-varying covariate’s effect is estimated from contrast between units at the same time: at each failure, did the unit that failed carry a different covariate value from the others still at risk? If every unit follows the same covariate path on the same clock, the covariate is perfectly confounded with time and its coefficient cannot be separated from the baseline. Staggered starts, different schedules and idle periods are what make the effect measurable.
Validating a survival predictor
Information criteria (AIC, BIC) compare how well models fit the data they were trained on. To judge how well a model predicts, it must be scored on held-out data, and the metrics must account for censoring. Two right-censored-standard measures are used, both handling censoring by inverse-probability-of-censoring weighting (IPCW):
The Brier score \(BS(t)\) is the weighted mean squared error between the predicted survival \(S(t \mid Z)\) and the survival indicator \(\mathbb{1}(T > t)\); the integrated Brier score averages it over a time grid. Lower is better, and 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 \(t\) was assigned a higher risk than one still event-free. 0.5 is chance, 1.0 is perfect.
IPCW works by up-weighting the subjects whose outcome is known at \(t\) by the inverse of their probability of not having been censored yet, estimated by a Kaplan-Meier of the censoring times, so the scored sample stands in for the full one. A subject still event-free at \(t\) gets weight \(1/\hat G(t)\); an event at \(x_i \le t\) gets \(1/\hat G(x_i-)\), the probability of being uncensored just before \(x_i\). An event and a censoring at the same time are recorded as an event, so the censoring Kaplan-Meier takes the event to come first (it is not at risk of being censored at its own time) and a censoring tied with an event does not count against it. The metrics are model-agnostic: they take a matrix of predicted survival probabilities, so parametric, semi-parametric and tree-based predictors can be compared on the same footing (see Comparison Tests and Validation Metrics).
Concordance
The oldest discrimination measure is Harrell’s concordance index [Harrell1982reg]. It asks, over every comparable pair of subjects, whether the model ranked them correctly: the one that failed first should have the higher risk score. A pair is comparable when the earlier time is an observed failure — if the earlier subject was censored we do not know who failed first. Then
with 0.5 for a random ranking and 1 for a perfect one. Ties in time need conventions. A failure and a censoring at the same time form a fully comparable pair, because the censored subject is known to have outlived the failure (so a lower score for the failure counts 0, not 0.5), and two censorings at the same time are not comparable. Two failures at the same time are where the conventions differ. By default surpyval follows Therneau, as R’s survival::concordance and lifelines do (ties="therneau"): neither subject outlived the other, so the pair is not comparable. Harrell’s original definition (ties="harrell") counts the pair, 1 if the scores tie and 0.5 otherwise. The two agree on data without tied failure times; for a Cox model of the lung data on age, sex and ECOG score, whose 28 pairs of tied deaths make the difference, Therneau’s gives 0.6371 (R and lifelines too) and Harrell’s 0.6369. For a proportional hazards model the linear predictor \(\beta' Z\) is a natural risk score. \(C\) measures ranking only — a model can have an excellent \(C\) and badly miscalibrated survival probabilities, which is what the Brier score is for. A four-subject example shows the counting:
from surpyval.metrics import concordance_index as score
x_c = [1.0, 2.0, 3.0, 4.0]
c_c = [0, 1, 0, 0] # the second subject is censored at 2
risk = [0.9, 0.2, 0.5, 0.1] # higher score = expected to fail earlier
# comparable pairs: (1,2) (1,3) (1,4) (3,4); (2,*) are not, 2 was censored
# all four are ranked correctly
score(x_c, c_c, risk)
1.0
Survival trees and forests
All the families above impose a structure — a link function, a linear predictor. A survival tree imposes none: it recursively splits the data on one covariate at a time, at the threshold that best separates the survival of the two halves, and fits a survival model in each leaf [LeBlancCrowley1993reg]. A prediction for a new unit drops it down the tree to a leaf and returns that leaf’s distribution. Trees find interactions and thresholds by themselves (“high temperature matters, but only for the old design”), at the cost of being step functions of the covariates and of being unstable — a slightly different sample can grow a different tree.
surpyval’s trees couple the split rule with the leaf model (kind=):
"non-parametric"scores a candidate split by the standardised log-rank statistic between the two children,\[L = \frac{\sum_j \bigl(d_{j,L} - Y_{j,L}\, d_j / Y_j\bigr)} {\sqrt{\sum_j \frac{Y_{j,L}}{Y_j}\Bigl(1 - \frac{Y_{j,L}}{Y_j}\Bigr)\frac{Y_j - d_j}{Y_j - 1}\, d_j}},\]summed over the pooled distinct times, with \(d\) deaths and \(Y\) numbers at risk (left child and total). The at-risk counts use the same
(entry, exit]convention as everywhere else, so left truncation is handled; the leaves are Nelson-Aalen estimates.A risk-set statistic only exists for observed and right-censored data. For left- and interval-censored data the tree uses the score form of the same test instead [Finkelstein1986reg]: at each node it fits the pooled Turnbull estimate \(\hat S\) once, gives each unit whose event lies in \((L, R]\) the log-rank score
\[c_i = \frac{\hat S(L)\log \hat S(L) - \hat S(R)\log \hat S(R)}{\hat S(L) - \hat S(R)}\](so \(\log \hat S(L)\) if right censored at \(L\), \(1 + \log \hat S(t)\) if observed at \(t\), and the \(L \to 0\) limit if left censored), and scores a split by the standardised sum of the left child’s scores, \(|\sum_{i \in L} c_i - n_L \bar c| / \sqrt{\tfrac{n_L n_R}{n(n-1)}\sum_i (c_i - \bar c)^2}\), with its permutation variance. \(\hat S\) is taken as \(e^{-\hat H}\), where \(\hat H\) is the Nelson-Aalen hazard of the Turnbull estimate’s expected risk sets and events; this keeps every score finite (a Kaplan-Meier type estimate reaches zero at a last observed failure, whose score would be \(-\infty\)). On right-censored data the scores are then exactly the log-rank scores \(\delta_i - \hat H(x_i)\), whose sum over a child is the log-rank numerator, so the two splits agree asymptotically. The leaves are Turnbull estimates. The same split is used for right-truncated data, whose units are seen only because they failed in time and so have no risk sets either. A unit observable only if its event falls in a window \((t_l, t_r]\) contributes the conditional likelihood \(P(L < T \le R)/P(t_l < T \le t_r)\), and its score is the score of that: with \(g(a, b)\) the score above for an event in \((a, b]\),
\[c_i = g\bigl(\max(L, t_l), \min(R, t_r)\bigr) - g(t_l, t_r),\]with \(\hat S\) from the Turnbull estimate fitted with the truncation. Under left truncation the window term is \(\log \hat S(t_l)\), so on observed and right-censored data the scores are the delayed-entry martingale residuals \(\delta_i - [\hat H(x_i) - \hat H(t_{l,i})]\) [FuSimonoff2017reg]. Subtracting the window’s score is what keeps the test honest when the truncation depends on a covariate: a score of a conditional likelihood has mean zero given the window, so a covariate that only delays entry, or brings forward the right-truncation time, is not mistaken for one that changes survival (in simulations, without the window term such a covariate was significant at the 5% level in 39% of data sets with delayed entry and 49% with right truncation; with it, in 0.5% and 3%). Under right truncation the data only identify the distribution of failures before the largest truncation time, and that is the distribution whose proportional-hazards alternative the split tests.
"weibull"(the default) and"exponential"score a split by the gain in the maximised full log-likelihood of a Weibull (or exponential [DavisAnderson1989reg]) model in each child. Because it uses the full likelihood of the proportional hazards section, this works for every kind of censoring and truncation; the leaves are the fitted Weibull or exponential models. On observed and right-censored data each child’s maximum is found directly rather than by an optimiser: the exponential rate is \(\hat\lambda = r / \sum_i n_i x_i\) (\(r\) the failures), and at a fixed Weibull shape \(\beta\) the likelihood is maximised by \(\hat\alpha^\beta = \sum_i n_i x_i^\beta / r\), which leaves the concave one-dimensional profile \(\ell_p(\beta) = r \log\beta - r \log(\sum_i n_i x_i^\beta / r) + (\beta - 1)\sum_{\text{failures}} n_i \log x_i - r\), maximised by Newton steps. Any split, even of pure noise, raises the in-sample likelihood, so these trees keep splitting until their leaves are too small: what a forest of deep trees wants, but not a single tree. A split adds the working model’s \(k\) parameters (1 for the exponential, 2 for the Weibull), andmin_split_gainasks for a gain \(\ell_L + \ell_R - \ell\) above a number, above \(k\) ("aic": the split lowers Akaike’s criterion) or above \(k \log(d) / 2\) ("bic", \(d\) the node’s failures). Neither is a test, since the split is the best of many cuts; the conditional-inference stop below is.
Choosing the split by the best cut over every covariate (selection="greedy", the default) has two known faults. A covariate with many distinct values offers many more cuts than one with a few, so by chance alone its best cut looks better: greedy search prefers continuous covariates whether or not they matter. And it always finds a cut, so the tree has to be stopped by its size. Conditional inference (selection="ctree") separates choosing the covariate from choosing its cut [Hothorn2006reg]. At each node every unit gets a score from the node’s split statistic, computed once: the log-rank scores above for "non-parametric" (with delayed entry \(\delta_i - [\hat H(x_i) - \hat H(t_{l,i})]\), whose sum over a child is again \(O - E\)), and for "exponential" and "weibull" each unit’s contribution to the score \(\partial \ell_i / \partial \theta\) of the working model at the node’s pooled estimate (one score for \(\log \lambda\), two for \((\log \alpha, \log \beta)\)), as in model-based recursive partitioning [Zeileis2008reg]. For a cut leaving \(m\) of \(N\) units on the left, the sum \(T_m\) of the left child’s centred scores has permutation covariance \(\frac{m(N-m)}{N(N-1)}\sum_i (h_i - \bar h)(h_i - \bar h)^\top\), and \(Q_m = T_m^\top \Sigma_m^{+} T_m\) is maximised over the covariate’s allowed cuts: the maximally selected statistic [LausenSchumacher1992reg], which, like the tree, depends on a covariate only through the order of its values. Its p-value accounts for the number of cuts: the standardised statistics at successive cuts form, asymptotically, a Gauss-Markov chain with correlation \(\sqrt{m_{k-1}(N - m_k)/(m_k(N - m_{k-1}))}\) between neighbours, and surpyval computes the probability that the chain’s maximum exceeds the observed one exactly for that chain, by quadrature, rather than by the Hunter-Worsley bound of maxstat [HothornLausen2003reg], which is about twice too large for a continuous covariate. One cut gives the plain \(\chi^2\) p-value, many close cuts not much more than a few independent ones. The covariate with the smallest p-value is chosen, and the node is split, at that covariate’s best cut by the kind’s criterion, only if the p-value times the number of covariates tested (Bonferroni) is below alpha_split (0.05 by default); otherwise the node is a leaf. On data with no effect the tree then usually stays a single leaf, and a covariate with two values competes on equal terms with continuous noise.
A random survival forest [Ishwaran2008reg] averages many trees, each grown on a bootstrap resample of the data and allowed to consider only a random subset of the covariates at each split. The averaging trades the high variance of a single deep tree for a little bias, and usually predicts much better. The forest’s survival curve is the average of the trees’ leaf survival curves (or, optionally, the survival implied by their averaged cumulative hazards), and its risk score for concordance is the leaf cumulative hazard summed over the evaluation times. Because every tree is grown without about a third of the rows, each row can be scored by the trees that never saw it: the forest’s out-of-bag log-likelihood is the mean of those rows’ full likelihoods (density, survival, failure or interval probability, over the truncation probability), an estimate of how well it predicts new data for every kind of censoring, and shuffling one covariate among the out-of-bag rows and measuring the drop gives its permutation importance [Breiman2001reg]. In surpyval both live in surpyval.beta.ml — tested and usable, but with an interface that may still change (see Machine Learning (beta)).
Choosing a model
There is no universally best family, but the questions below settle most cases.
Is the effect a constant multiple of the risk? Fit a Cox model and run the Grambsch-Therneau test. If it passes, proportional hazards is a defensible and very interpretable default, and a Cox fit makes no assumption about the baseline.
Do you need to extrapolate in time, or predict the full distribution? Cox cannot predict beyond the observed times. Use a parametric family, and choose the baseline by likelihood (AIC/BIC) and by comparing the fitted curves with non-parametric estimates.
Is the mechanism “the covariate speeds the clock up”? Use AFT — a stress that accelerates a physical process — or, when the covariate is a controlled stress with a known physical law and you must extrapolate to use conditions, accelerated life. Remember that for a Weibull, AFT and PH are the same model.
Does the effect fade over time? Proportional odds represents exactly that; a PH model forced onto such data reports a time-averaged hazard ratio.
Is the question about excess risk? Additive hazards reports risk differences, in units of failures per unit time.
Are units grouped? For honest standard errors, cluster-robust errors; to model and quantify the between-group variation (and to predict for a known group), a shared frailty; for a nuisance grouping that violates PH, stratification.
Do covariates change during follow-up? The start-stop format, with Cox, PH, AH or AFT.
Is the covariate structure unknown or strongly non-linear? A random survival forest, validated on held-out data, and compared against a simpler model on the same metrics.
Whatever the choice, check it: residuals and the PH test for the assumption, information criteria for the fit, and held-out Brier scores, AUC and concordance for prediction. For worked examples on how to do regression analysis — including checking the proportional hazards assumption, robust and stratified fits, and validating predictions — see the Regression Modelling with SurPyval page.
References
Cox, D.R., 1972. Regression models and life-tables. Journal of the Royal Statistical Society: Series B, 34(2), pp.187-220.
Breslow, N., 1974. Covariance analysis of censored survival data. Biometrics, 30(1), pp.89-99.
Efron, B., 1977. The efficiency of Cox’s likelihood function for censored data. Journal of the American Statistical Association, 72(359), pp.557-565.
DeLong, D.M., Guirguis, G.H. and So, Y.C., 1994. Efficient computation of subset selection probabilities with application to Cox regression. Biometrika, 81(3), pp.607-611.
Gail, M.H., Lubin, J.H. and Rubinstein, L.V., 1981. Likelihood calculations for matched case-control studies and survival studies with tied death times. Biometrika, 68(3), pp.703-707.
Kalbfleisch, J.D. and Prentice, R.L., 2002. The Statistical Analysis of Failure Time Data, 2nd ed. Wiley.
Bennett, S., 1983. Analysis of survival data by the proportional odds model. Statistics in Medicine, 2(2), pp.273-277.
Lin, D.Y. and Ying, Z., 1994. Semiparametric analysis of the additive risk model. Biometrika, 81(1), pp.61-71.
Buckley, J. and James, I., 1979. Linear regression with censored data. Biometrika, 66(3), pp.429-436.
Grambsch, P.M. and Therneau, T.M., 1994. Proportional hazards tests and diagnostics based on weighted residuals. Biometrika, 81(3), pp.515-526.
Therneau, T.M. and Grambsch, P.M., 2000. Modeling Survival Data: Extending the Cox Model. Springer.
Lin, D.Y. and Wei, L.J., 1989. The robust inference for the Cox proportional hazards model. Journal of the American Statistical Association, 84(408), pp.1074-1078.
Hougaard, P., 2000. Analysis of Multivariate Survival Data. Springer.
Harrell, F.E., Califf, R.M., Pryor, D.B., Lee, K.L. and Rosati, R.A., 1982. Evaluating the yield of medical tests. JAMA, 247(18), pp.2543-2546.
LeBlanc, M. and Crowley, J., 1993. Survival trees by goodness of split. Journal of the American Statistical Association, 88(422), pp.457-467.
Finkelstein, D.M., 1986. A proportional hazards model for interval-censored failure time data. Biometrics, 42(4), pp.845-854.
Fu, W. and Simonoff, J.S., 2017. Survival trees for left-truncated and right-censored data, with application to time-varying covariate data. Biostatistics, 18(2), pp.352-369.
Breiman, L., 2001. Random forests. Machine Learning, 45(1), pp.5-32.
Davis, R.B. and Anderson, J.R., 1989. Exponential survival trees. Statistics in Medicine, 8(8), pp.947-961.
Ishwaran, H., Kogalur, U.B., Blackstone, E.H. and Lauer, M.S., 2008. Random survival forests. The Annals of Applied Statistics, 2(3), pp.841-860.
Hothorn, T., Hornik, K. and Zeileis, A., 2006. Unbiased recursive partitioning: a conditional inference framework. Journal of Computational and Graphical Statistics, 15(3), pp.651-674.
Zeileis, A., Hothorn, T. and Hornik, K., 2008. Model-based recursive partitioning. Journal of Computational and Graphical Statistics, 17(2), pp.492-514.
Lausen, B. and Schumacher, M., 1992. Maximally selected rank statistics. Biometrics, 48(1), pp.73-85.
Hothorn, T. and Lausen, B., 2003. On the exact distribution of maximally selected rank statistics. Computational Statistics & Data Analysis, 43(2), pp.121-137.