|
| 1 | +# -*- coding: utf-8 -*- |
| 2 | +from lifelines.fitters import ParametricRegressionFitter |
| 3 | +from autograd.scipy.special import expit |
| 4 | +from autograd import numpy as np |
| 5 | +from lifelines.utils.safe_exp import safe_exp |
| 6 | + |
| 7 | +exp = safe_exp |
| 8 | +dot = np.dot |
| 9 | + |
| 10 | + |
| 11 | +class CopulaFrailtyWeilbullModel(ParametricRegressionFitter): |
| 12 | + """ |
| 13 | +
|
| 14 | + Assume that subjects have two competing risks that are _not_ independent. |
| 15 | + Assume that their is a frailty term that both risks depend on, and, conditional |
| 16 | + on this term, the risks are IID Weibull models. If the frailty term is |
| 17 | + stable with parameter alpha, then the cumulative hazard can be written down explicitly. |
| 18 | +
|
| 19 | +
|
| 20 | +
|
| 21 | + Reference |
| 22 | + -------------- |
| 23 | + Frees and Valdez, UNDERSTANDING RELATIONSHIPS USING COPULAS |
| 24 | +
|
| 25 | + """ |
| 26 | + |
| 27 | + _fitted_parameter_names = ["lambda1", "rho1", "lambda2", "rho2", "alpha"] |
| 28 | + |
| 29 | + def _cumulative_hazard(self, params, T, Xs): |
| 30 | + lambda1 = exp(dot(Xs["lambda1"], params["lambda1"])) |
| 31 | + lambda2 = exp(dot(Xs["lambda2"], params["lambda2"])) |
| 32 | + rho2 = exp(dot(Xs["rho2"], params["rho2"])) |
| 33 | + rho1 = exp(dot(Xs["rho1"], params["rho1"])) |
| 34 | + alpha = exp(dot(Xs["alpha"], params["alpha"])) |
| 35 | + |
| 36 | + return ((T / lambda1) ** rho1 + (T / lambda2) ** rho2) ** alpha |
| 37 | + |
| 38 | + |
| 39 | +swf = CopulaFrailtyWeilbullModel(penalizer=0.1) |
| 40 | + |
| 41 | +rossi = load_rossi() |
| 42 | +rossi["intercept"] = 1.0 |
| 43 | +rossi["week"] = rossi["week"] / 54.0 |
| 44 | + |
| 45 | +covariates = { |
| 46 | + "lambda1": rossi.columns, |
| 47 | + "lambda2": rossi.columns, |
| 48 | + "rho1": ["intercept"], |
| 49 | + "rho2": ["intercept"], |
| 50 | + "alpha": ["intercept"], |
| 51 | +} |
| 52 | + |
| 53 | +swf.fit(rossi, "week", event_col="arrest", regressors=covariates, timeline=np.linspace(0, 2)) |
| 54 | +swf.print_summary(2) |
0 commit comments