Custom regression models¶

Like for univariate models, it is possible to create your own custom parametric survival models. Why might you want to do this?

Create new / extend AFT models using known probability distributions
Create a piecewise model using domain knowledge about subjects
Iterate and fit a more accurate parametric model

lifelines has a very simple API to create custom parametric regression models. You only need to define the cumulative hazard function. For example, the cumulative hazard for the constant-hazard regression model looks like:

$$ H(t, x) = \frac{t}{\lambda(x)}\\ \lambda(x) = \exp{(\vec{\beta} \cdot \vec{x}^{\,T})} $$

where $\beta$ are the unknowns we will optimize over.

Below are some example custom models.

In [9]:

from lifelines.fitters import ParametricRegressionFitter
from autograd import numpy as np
from lifelines.datasets import load_rossi
%config InlineBackend.figure_format = 'retina'


class ExponentialAFTFitter(ParametricRegressionFitter):
    
    # this class property is necessary, and should always be a non-empty list of strings. 
    _fitted_parameter_names = ['lambda_']
    
    def _cumulative_hazard(self, params, t, Xs):
        # params is a dictionary that maps unknown parameters to a numpy vector. 
        # Xs is a dictionary that maps unknown parameters to a numpy 2d array 
        beta = params['lambda_']
        X = Xs['lambda_']
        lambda_ = np.exp(np.dot(X, beta))
        return t / lambda_
    

rossi = load_rossi()
rossi['intercept'] = 1.0

# the below variables maps dataframe columns to parameters
regressors = {
    'lambda_': rossi.columns
}

eaf = ExponentialAFTFitter().fit(rossi, 'week', 'arrest', regressors=regressors)
eaf.print_summary()

model	lifelines.ExponentialAFTFitter
duration col	'week'
event col	'arrest'
number of observations	432
number of events observed	114
log-likelihood	-686.37
time fit was run	2020-06-21 17:36:47 UTC

		coef	exp(coef)	se(coef)	coef lower 95%	coef upper 95%	exp(coef) lower 95%	exp(coef) upper 95%	z	p	-log2(p)
lambda_	fin	0.37	1.44	0.19	-0.01	0.74	0.99	2.10	1.92	0.06	4.18
	age	0.06	1.06	0.02	0.01	0.10	1.01	1.10	2.55	0.01	6.52
	race	-0.30	0.74	0.31	-0.91	0.30	0.40	1.35	-0.99	0.32	1.63
	wexp	0.15	1.16	0.21	-0.27	0.56	0.76	1.75	0.69	0.49	1.03
	mar	0.43	1.53	0.38	-0.32	1.17	0.73	3.24	1.12	0.26	1.93
	paro	0.08	1.09	0.20	-0.30	0.47	0.74	1.59	0.42	0.67	0.57
	prio	-0.09	0.92	0.03	-0.14	-0.03	0.87	0.97	-3.03	<0.005	8.65
	intercept	4.05	57.44	0.59	2.90	5.20	18.21	181.15	6.91	<0.005	37.61

AIC	1388.73
log-likelihood ratio test	31.22 on 6 df
-log2(p) of ll-ratio test	15.41

In [10]:

from lifelines.calibration import survival_probability_calibration

fig, ax = plt.subplots(figsize=(8, 5))

survival_probability_calibration(eaf, rossi, 25, ax=ax)

ICI =  0.027226192838563718
E50 =  0.02894862605432391

Out[10]:

(<matplotlib.axes._subplots.AxesSubplot at 0x11af2ea10>,
 0.027226192838563718,
 0.02894862605432391)

Cure models¶

Suppose in our population we have a subpopulation that will never experience the event of interest. Or, for some subjects the event will occur so far in the future that it's essentially at time infinity. In this case, the survival function for an individual should not asymptically approach zero, but some positive value. Models that describe this are sometimes called cure models (i.e. the subject is "cured" of death and hence no longer susceptible) or time-lagged conversion models.

It would be nice to be able to use common survival models and have some "cure" component. Let's suppose that for individuals that will experience the event of interest, their survival distrubtion is a Weibull, denoted $S_W(t)$. For a random selected individual in the population, thier survival curve, $S(t)$, is:

$$ \begin{align*} S(t) = P(T > t) &= P(\text{cured}) P(T > t\;|\;\text{cured}) + P(\text{not cured}) P(T > t\;|\;\text{not cured}) \\ &= p + (1-p) S_W(t) \end{align*} $$

Even though it's in an unconvential form, we can still determine the cumulative hazard (which is the negative logarithm of the survival function):

$$ H(t) = -\log{\left(p + (1-p) S_W(t)\right)} $$

In [11]:

from autograd.scipy.special import expit

class CureModel(ParametricRegressionFitter):
    _scipy_fit_method = "SLSQP"
    _scipy_fit_options = {"ftol": 1e-10, "maxiter": 200}

    _fitted_parameter_names = ["lambda_", "beta_", "rho_"]

    def _cumulative_hazard(self, params, T, Xs):
        c = expit(np.dot(Xs["beta_"], params["beta_"]))

        lambda_ = np.exp(np.dot(Xs["lambda_"], params["lambda_"]))
        rho_ = np.exp(np.dot(Xs["rho_"], params["rho_"]))
        sf = np.exp(-(T / lambda_) ** rho_)

        return -np.log((1 - c) + c * sf)


cm = CureModel(penalizer=0.0)

rossi = load_rossi()
rossi["intercept"] = 1.0

covariates = {"lambda_": rossi.columns, "rho_": ["intercept"], "beta_": ['intercept', 'fin']}

cm.fit(rossi, "week", event_col="arrest", regressors=covariates, timeline=np.arange(250)) 
cm.print_summary(2)

model	lifelines.CureModel
duration col	'week'
event col	'arrest'
number of observations	432
number of events observed	114
log-likelihood	-679.51
time fit was run	2020-06-21 17:36:50 UTC

		coef	exp(coef)	se(coef)	coef lower 95%	coef upper 95%	exp(coef) lower 95%	exp(coef) upper 95%	z	p	-log2(p)
beta_	fin	-0.41	0.66	1.09	-2.55	1.73	0.08	5.64	-0.38	0.71	0.50
beta_	intercept	0.70	2.02	0.98	-1.21	2.62	0.30	13.73	0.72	0.47	1.08
lambda_	fin	0.15	1.16	0.38	-0.60	0.90	0.55	2.46	0.39	0.69	0.52
	age	0.03	1.04	0.02	0.00	0.07	1.00	1.07	2.10	0.04	4.82
	race	-0.28	0.75	0.23	-0.73	0.16	0.48	1.18	-1.25	0.21	2.24
	wexp	0.20	1.23	0.18	-0.14	0.55	0.87	1.73	1.16	0.25	2.02
	mar	0.34	1.41	0.26	-0.17	0.86	0.84	2.36	1.30	0.19	2.37
	paro	0.03	1.04	0.15	-0.27	0.34	0.77	1.40	0.23	0.82	0.29
	prio	-0.08	0.93	0.02	-0.12	-0.03	0.89	0.97	-3.23	<0.005	9.64
	intercept	3.76	43.15	0.49	2.80	4.73	16.44	113.22	7.65	<0.005	45.48
rho_	intercept	0.43	1.54	0.12	0.20	0.66	1.22	1.94	3.64	<0.005	11.87

AIC	1381.03
log-likelihood ratio test	34.22 on 9 df
-log2(p) of ll-ratio test	13.58

In [12]:

cm.predict_survival_function(rossi.loc[::100]).plot(figsize=(12,6))

Out[12]:

<matplotlib.axes._subplots.AxesSubplot at 0x11b149f50>

In [13]:

# what's the effect on the survival curve if I vary "age"
fig, ax = plt.subplots(figsize=(12, 6))

cm.plot_covariate_groups(['age'], values=np.arange(20, 50, 5), cmap='coolwarm', ax=ax)

Out[13]:

<matplotlib.axes._subplots.AxesSubplot at 0x11a96b510>

In [14]:

from lifelines.calibration import survival_probability_calibration

fig, ax = plt.subplots(figsize=(8, 5))

survival_probability_calibration(cm, rossi, 25, ax=ax)

ICI =  0.0025660960118482753
E50 =  0.0018552649672565757

Out[14]:

(<matplotlib.axes._subplots.AxesSubplot at 0x11a969d50>,
 0.0025660960118482753,
 0.0018552649672565757)

Spline models¶

See royston_parmar_splines.py and crowther_royston_clements_splines.py in the examples folder: https://github.com/CamDavidsonPilon/lifelines/tree/master/examples

Note in the below model the use of _create_initial_point, and one of the parameters is non-zero initially. This is important as it nudges the model slightly away from the degenerate all-zeros model. Try setting it to 0, and watch the model fail to converge.

In [15]:

class SplineFitter:
    # this is also available in lifelines.fitters.mixins - it's reproduced here for example's sake.
    
    @staticmethod
    def relu(x):
        return np.maximum(0, x)

    def basis(self, x, knot, min_knot, max_knot):
        lambda_ = (max_knot - knot) / (max_knot - min_knot)
        return self.relu(x - knot) ** 3 - (
            lambda_ * self.relu(x - min_knot) ** 3 + (1 - lambda_) * self.relu(x - max_knot) ** 3
        )

class PHSplineFitter(SplineFitter, ParametricRegressionFitter):
    """
    Proportional Hazard model with baseline modelled as a spline

    References
    ------------
    Royston, P., & Parmar, M. K. B. (2002). Flexible parametric proportional-hazards and proportional-odds models for censored survival data, with application to prognostic modelling and estimation of treatment effects. Statistics in Medicine, 21(15), 2175–2197. doi:10.1002/sim.1203 
    """

    _fitted_parameter_names = ["beta_", "phi0_", "phi1_", "phi2_"]

    def _create_initial_point(self, Ts, E, entries, weights, Xs):
        return {
             "beta_": np.zeros(len(Xs.mappings["beta_"])), 
             "phi0_": np.array([0.0]), 
             "phi1_": np.array([0.1]),
             "phi2_": np.array([0.0])
            }
        
    
    def set_knots(self, T, E):
        self.knots = np.percentile(np.log(T[E.astype(bool).values]), np.linspace(5, 95, 3))

    def _pre_fit_model(self, Ts, E, df):
        # this function runs before the model is fit and can be used to set data-determined values (like knots)
        self.set_knots(Ts[0], E)
        
    def _cumulative_hazard(self, params, T, Xs):
        lT = np.log(T)
        return np.exp(
            np.dot(Xs["beta_"], params["beta_"]) 
            + params["phi0_"]
            + (params["phi1_"]) * lT
            + params["phi2_"] * self.basis(lT, self.knots[1], self.knots[0], self.knots[-1])
        )
    
    
rossi = load_rossi()
rossi["intercept"] = 1.0

covariates = {"beta_": rossi.columns.difference(['intercept', 'arrest', 'week']), 
              "phi0_": ["intercept"], 
              "phi1_": ["intercept"], 
              "phi2_": ["intercept"], 
             }

phf = PHSplineFitter(penalizer=0.0)
phf.fit(rossi, "week", "arrest", regressors=covariates) 
phf.print_summary(2)

model	lifelines.PHSplineFitter
duration col	'week'
event col	'arrest'
number of observations	432
number of events observed	114
log-likelihood	-679.86
time fit was run	2020-06-21 17:36:55 UTC

		coef	exp(coef)	se(coef)	coef lower 95%	coef upper 95%	exp(coef) lower 95%	exp(coef) upper 95%	z	p	-log2(p)
beta_	fin	-0.38	0.68	0.19	-0.76	-0.01	0.47	0.99	-1.99	0.05	4.43
	age	-0.06	0.94	0.02	-0.10	-0.01	0.90	0.99	-2.60	0.01	6.73
	race	0.32	1.37	0.31	-0.29	0.92	0.75	2.51	1.03	0.31	1.71
	wexp	-0.15	0.86	0.21	-0.57	0.27	0.57	1.31	-0.71	0.48	1.06
	mar	-0.44	0.65	0.38	-1.19	0.31	0.31	1.37	-1.14	0.25	1.99
	paro	-0.08	0.92	0.20	-0.47	0.30	0.63	1.35	-0.42	0.67	0.57
	prio	0.09	1.10	0.03	0.04	0.15	1.04	1.16	3.21	<0.005	9.56
phi0_	intercept	-5.81	0.00	0.98	-7.72	-3.89	0.00	0.02	-5.94	<0.005	28.42
phi1_	intercept	1.48	4.41	0.27	0.96	2.01	2.61	7.43	5.56	<0.005	25.13
phi2_	intercept	0.05	1.05	0.14	-0.23	0.33	0.79	1.40	0.34	0.73	0.45

AIC	1379.71
log-likelihood ratio test	33.31 on 8 df
-log2(p) of ll-ratio test	14.17

In [16]:

from lifelines.calibration import survival_probability_calibration

fig, ax = plt.subplots(figsize=(8, 5))

survival_probability_calibration(phf, rossi, 25, ax=ax)

ICI =  0.005080507423341433
E50 =  0.00440577738638781

Out[16]:

(<matplotlib.axes._subplots.AxesSubplot at 0x11b1e7590>,
 0.005080507423341433,
 0.00440577738638781)