Régression de survie – documentation Lifelines 0.25.8

Ajuster la régression

L’implémentation du modèle Cox dans lignes de vie est sous CoxPHFitter. Nous adaptons le modèle à l’ensemble de données en utilisant fit(). Il a un print_summary() fonction qui imprime une vue tabulaire des coefficients et des statistiques associées.

from lifelines import CoxPHFitter
from lifelines.datasets import load_rossi

rossi = load_rossi()

cph = CoxPHFitter()
cph.fit(rossi, duration_col='week', event_col='arrest')

cph.print_summary()  # access the individual results using cph.summary

"""

             duration col = 'week'
                event col = 'arrest'
   number of observations = 432
number of events observed = 114
   partial log-likelihood = -658.75
         time fit was run = 2019-10-05 14:24:44 UTC

---
       coef  exp(coef)   se(coef)   coef lower 95%   coef upper 95%  exp(coef) lower 95%  exp(coef) upper 95%
fin   -0.38       0.68       0.19            -0.75            -0.00                 0.47                 1.00
age   -0.06       0.94       0.02            -0.10            -0.01                 0.90                 0.99
race   0.31       1.37       0.31            -0.29             0.92                 0.75                 2.50
wexp  -0.15       0.86       0.21            -0.57             0.27                 0.57                 1.30
mar   -0.43       0.65       0.38            -1.18             0.31                 0.31                 1.37
paro  -0.08       0.92       0.20            -0.47             0.30                 0.63                 1.35
prio   0.09       1.10       0.03             0.04             0.15                 1.04                 1.16

         z      p   -log2(p)
fin  -1.98   0.05       4.40
age  -2.61   0.01       6.79
race  1.02   0.31       1.70
wexp -0.71   0.48       1.06
mar  -1.14   0.26       1.97
paro -0.43   0.66       0.59
prio  3.19 <0.005       9.48
---
Concordance = 0.64
Partial AIC = 1331.50
log-likelihood ratio test = 33.27 on 7 df
-log2(p) of ll-ratio test = 15.37
"""

Nouveau dans la v0.25.0, nous pouvons également utiliser des «formules» pour gérer le côté droit du modèle linéaire. Par exemple:

cph.fit(rossi, duration_col='week', event_col='arrest', formula="fin + wexp + age * prio")

est analogue au modèle linéaire avec terme d’interaction:

cph.fit(rossi, duration_col='week', event_col='arrest', formula="fin + wexp + age * prio")
cph.print_summary()

"""

             duration col = 'week'
                event col = 'arrest'
      baseline estimation = breslow
   number of observations = 432
number of events observed = 114
   partial log-likelihood = -659.39
         time fit was run = 2020-07-13 19:30:33 UTC

---
            coef  exp(coef)   se(coef)   coef lower 95%   coef upper 95%  exp(coef) lower 95%  exp(coef) upper 95%
covariate
fin        -0.33       0.72       0.19            -0.70             0.04                 0.49                 1.05
wexp       -0.24       0.79       0.21            -0.65             0.17                 0.52                 1.19
age        -0.03       0.97       0.03            -0.09             0.03                 0.92                 1.03
prio        0.31       1.36       0.17            -0.03             0.64                 0.97                 1.90
age:prio   -0.01       0.99       0.01            -0.02             0.01                 0.98                 1.01

              z    p   -log2(p)
covariate
fin       -1.73 0.08       3.57
wexp      -1.14 0.26       1.97
age       -0.93 0.35       1.51
prio       1.80 0.07       3.80
age:prio  -1.28 0.20       2.32
---
Concordance = 0.64
Partial AIC = 1328.77
log-likelihood ratio test = 31.99 on 5 df
-log2(p) of ll-ratio test = 17.35
"""

Les formules peuvent être utilisées pour créer des interactions, encoder des variables catégorielles, créer des splines de base, etc. Les formules utilisées sont (presque) les mêmes que celles disponibles dans les modèles R et stats.

Interprétation

Pour accéder directement aux coefficients et à l’aléa de référence, vous pouvez utiliser params_ et baseline_hazard_ respectivement. Jetant un œil à ces coefficients pendant un moment, prio (le nombre d’arrestations antérieures) a un coefficient d’environ 0,09. Ainsi, une augmentation d’une unité de prio signifie que le danger de référence augmentera d’un facteur de ( exp {(0,09)} = 1,10 ) – une augmentation d’environ 10%. Rappelons que dans le modèle de risque proportionnel de Cox, un danger plus élevé signifie plus de risque que l’événement se produise. La valeur ( exp {(0,09)} ) s’appelle le Taux de dangerosité, un nom qui sera clair avec un autre exemple.

Considérez le coefficient de mar (que le sujet soit marié ou non). Les valeurs de la colonne sont binaires: 0 ou 1, représentant célibataire ou marié. La valeur du coefficient associé à mar, ( exp {(-. 43)} ), est la valeur du ratio de dangers associé au mariage, c’est-à-dire:

Notez que le côté gauche est une constante (en particulier, il est indépendant du temps, (t )), mais le côté droit a deux facteurs qui peuvent varier avec le temps. La hypothèse de risque proportionnel est-ce que cette relation est vraie. Autrement dit, les dangers peuvent changer avec le temps, mais leur rapport entre les niveaux reste constant. Plus tard, nous traiterons de vérifier cette hypothèse. Cependant, en réalité, il est très courant que le rapport de risque change au cours de la durée de l’étude. Le rapport de risque a alors l’interprétation d’une sorte de moyenne pondérée des rapports de risque spécifiques à la période. En conséquence, le rapport de risque peut dépendre de manière critique de la durée du suivi.

Convergence

L’ajustement du modèle de Cox aux données implique l’utilisation de méthodes itératives. lignes de vie demande des efforts supplémentaires pour aider à la convergence, soyez donc attentif à tout avertissement qui apparaît. La correction des avertissements aidera généralement la convergence et diminuera le nombre d’étapes itératives requises. Si vous souhaitez voir plus d’informations lors du montage, il y a un show_progress paramètre dans fit() fonction. Pour plus d’aide, voir Problèmes de convergence dans le modèle de risque proportionnel de Cox.

Après l’ajustement, la valeur de la log-vraisemblance maximale disponible en utilisant log_likelihood_. La matrice de variance des coefficients est disponible sous variance_matrix_.

Qualité de l’ajustement

Après l’ajustement, vous voudrez peut-être savoir dans quelle mesure votre modèle correspondait aux données. Quelques méthodes que l’auteur a trouvées utiles sont de

Pénalités et régression éparse

Il est également possible d’ajouter un terme pénalisant à la régression de Cox. On peut les utiliser pour i) stabiliser les coefficients, ii) réduire les estimations à 0, iii) encourager un point de vue bayésien et iv) créer des coefficients clairsemés. Tous les modèles de régression, y compris le modèle de Cox, incluent à la fois une pénalité L1 et L2:

Pour utiliser ceci dans lignes de vie, à la fois le penalizer et l1_ratio peut être spécifié dans la création de classe:

from lifelines import CoxPHFitter
from lifelines.datasets import load_rossi

rossi = load_rossi()

cph = CoxPHFitter(penalizer=0.1, l1_ratio=1.0) # sparse solutions,
cph.fit(rossi, 'week', 'arrest')
cph.print_summary()

Au lieu d’un flotteur, un tableau peut être fourni avec la même taille que le nombre de paramètres pénalisés. Les valeurs du tableau sont des coefficients de pénalité spécifiques pour chaque covariable. Ceci est utile pour une structure de covariable plus compliquée. Quelques exemples:

1. vous avez beaucoup de facteurs de confusion que vous souhaitez pénaliser, mais pas le (s) traitement (s) principal (s).

from lifelines import CoxPHFitter
from lifelines.datasets import load_rossi

rossi = load_rossi()

# variable `fin` is the treatment of interest so don't penalize it at all
penalty = np.array([0, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5])

cph = CoxPHFitter(penalizer=penalty)
cph.fit(rossi, 'week', 'arrest')
cph.print_summary()

2. vous devez fusionner les catégories ensemble.
3. vous souhaitez implémenter un solution très rare.

En savoir plus sur les sanctions et leur mise en œuvre sur notre blog de développement.

Tracer les coefficients

Avec un modèle ajusté, une autre façon de visualiser les coefficients et leurs plages est d’utiliser le plot méthode.

from lifelines.datasets import load_rossi
from lifelines import CoxPHFitter

rossi = load_rossi()
cph = CoxPHFitter()
cph.fit(rossi, duration_col='week', event_col='arrest')

cph.plot()

Tracer l’effet de la variation d’une covariable

Après l’ajustement, nous pouvons tracer à quoi ressemblent les courbes de survie en faisant varier une seule covariable tout en
tenant tout le reste égal. Ceci est utile pour comprendre l’impact d’une covariable, étant donné le modèle. Pour ce faire, nous utilisons le plot_partial_effects_on_outcome() et donnez-lui la covariable d’intérêt, ainsi que les valeurs à afficher.

from lifelines.datasets import load_rossi
from lifelines import CoxPHFitter

rossi = load_rossi()
cph = CoxPHFitter()
cph.fit(rossi, duration_col='week', event_col='arrest')

cph.plot_partial_effects_on_outcome(covariates='prio', values=[0, 2, 4, 6, 8, 10], cmap='coolwarm')

S’il existe des entités dérivées dans votre jeu de données, par exemple, supposons que vous ayez inclus prio et prio**2 dans votre ensemble de données. Cela n’a pas de sens de simplement varier year et part year**2 fixé. Vous devrez spécifier manuellement les valeurs prises par les covariables dans un tableau ou une liste N-d (où N est le nombre de covariables à modifier.)

rossi['prio**2'] = rossi['prio'] ** 2

cph.fit(rossi, 'week', 'arrest')

cph.plot_partial_effects_on_outcome(
    covariates=['prio', 'prio**2'],
    values=[
        [0, 0],
        [1, 1],
        [2, 4],
        [3, 9],
        [8, 64],
    ],
    cmap='coolwarm')

Cependant, si vous avez utilisé le formula kwarg in fit, toutes les transformations nécessaires seront effectuées en interne pour vous.

cph.fit(rossi, 'week', 'arrest', formula="prio + I(prio**2)")

cph.plot_partial_effects_on_outcome(
    covariates=['prio'],
    values=[0, 1, 2, 3, 8],
    cmap='coolwarm')

Cette fonctionnalité est également utile pour analyser les variables catégorielles:

cph.plot_partial_effects_on_outcome(
    covariates=["a_categorical_variable"]
    values=["A", "B", ...],
    plot_baseline=False)

Vérification de l’hypothèse des risques proportionnels

Pour faire des inférences appropriées, nous devons nous demander si notre modèle de Cox est approprié pour notre ensemble de données. Rappelons ci-dessus que lorsque nous utilisons le modèle de Cox, nous appliquons implicitement l’hypothèse de risque proportionnel. Nous devrions nous demander si notre ensemble de données obéit à cette hypothèse?

CoxPHFitter a un check_assumptions() méthode qui produira des violations de l’hypothèse de risque proportionnel. Pour obtenir un didacticiel sur la résolution des violations, consultez Test des hypothèses de risque proportionnel. Les suggestions consistent à rechercher des moyens de stratifier une colonne (voir la documentation ci-dessous), ou utilisez un modèle variant dans le temps.

Stratification

Parfois, une ou plusieurs covariables peuvent ne pas obéir à l’hypothèse du risque proportionnel. Dans ce cas, nous pouvons permettre que la ou les covariables soient toujours incluses dans le modèle sans en estimer l’effet. C’est ce qu’on appelle la stratification. À un niveau élevé, considérez cela comme une division du jeu de données en m ensembles de données plus petits, partitionnés par les valeurs uniques de la ou des covariables de stratification. Chaque jeu de données a son propre risque de référence (la partie non paramétrique du modèle), mais ils partagent tous les paramètres de régression (la partie paramétrique du modèle). Puisque les covariables sont les mêmes dans chaque ensemble de données, il n’y a pas de paramètre de régression pour les covariables stratifiées sur, par conséquent, elles n’apparaîtront pas dans la sortie. Cependant il y aura m dangers de base sous baseline_cumulative_hazard_.

Pour spécifier les variables à utiliser dans la stratification, nous les définissons dans l’appel à fit():

from lifelines.datasets import load_rossi
from lifelines import CoxPHFitter
rossi = load_rossi()

cph = CoxPHFitter()
cph.fit(rossi, 'week', event_col='arrest', strata=['wexp'])
cph.print_summary()

"""

             duration col = 'week'
                event col = 'arrest'
                   strata = ['wexp']
      baseline estimation = breslow
   number of observations = 432
number of events observed = 114
   partial log-likelihood = -580.89
         time fit was run = 2020-08-09 21:25:37 UTC

---
            coef  exp(coef)   se(coef)   coef lower 95%   coef upper 95%  exp(coef) lower 95%  exp(coef) upper 95%
covariate
fin        -0.38       0.68       0.19            -0.76            -0.01                 0.47                 0.99
age        -0.06       0.94       0.02            -0.10            -0.01                 0.90                 0.99
race        0.31       1.36       0.31            -0.30             0.91                 0.74                 2.49
mar        -0.45       0.64       0.38            -1.20             0.29                 0.30                 1.34
paro       -0.08       0.92       0.20            -0.47             0.30                 0.63                 1.35
prio        0.09       1.09       0.03             0.03             0.15                 1.04                 1.16
              z      p   -log2(p)
covariate
fin       -1.99   0.05       4.42
age       -2.64   0.01       6.91
race       1.00   0.32       1.65
mar       -1.19   0.23       2.09
paro      -0.42   0.67       0.57
prio       3.16 <0.005       9.33
---
Concordance = 0.61
Partial AIC = 1173.77
log-likelihood ratio test = 23.77 on 6 df
-log2(p) of ll-ratio test = 10.77

"""

cph.baseline_survival_.shape
# (49, 2)
cph.baseline_cumulative_hazard_.plot(drawstyle="steps")

Poids et erreurs robustes

Les observations peuvent également être accompagnées de poids. Ces poids peuvent être des valeurs entières représentant une observation courante ou des valeurs flottantes représentant certains poids d’échantillonnage (ex: poids de probabilité inverses). dans le fit() méthode, un kwarg est présent pour spécifier quelle colonne du DataFrame doit être utilisée comme poids, ex: CoxPHFitter(df, 'T', 'E', weights_col='weights').

Lorsque vous utilisez des poids d’échantillonnage, il est correct de modifier également les calculs d’erreur standard. Cela se fait en allumant le robust drapeau dans fit(). Intérieurement, CoxPHFitter utilisera l’estimateur sandwich pour calculer les erreurs.

import pandas as pd
from lifelines import CoxPHFitter

df = pd.DataFrame({
    'T': [5, 3, 9, 8, 7, 4, 4, 3, 2, 5, 6, 7],
    'E': [1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0],
    'weights': [1.1, 0.5, 2.0, 1.6, 1.2, 4.3, 1.4, 4.5, 3.0, 3.2, 0.4, 6.2],
    'month': [10, 3, 9, 8, 7, 4, 4, 3, 2, 5, 6, 7],
    'age': [4, 3, 9, 8, 7, 4, 4, 3, 2, 5, 6, 7],
})

cph = CoxPHFitter()
cph.fit(df, 'T', 'E', weights_col='weights', robust=True)
cph.print_summary()

Voir plus d’exemples dans Ajout de poids aux observations dans un modèle de Cox.

Clusters et corrélations

Les groupes de sujets liés sont une autre propriété de votre ensemble de données. Cela peut être causé par:

• un seul individu ayant plusieurs occurrences, et par conséquent apparaissant dans l’ensemble de données plus d’une fois.
• les sujets qui partagent certains biens communs, comme les membres de la même famille ou qui sont appariés sur des scores de propension.

Nous appelons ces sujets groupés «clusters» et supposons qu’ils sont désignés par une colonne dans le DataFrame (exemple ci-dessous). Lorsque vous utilisez un cluster, les estimations ponctuelles du modèle ne changent pas, mais les erreurs standard augmenteront. Un argument intuitif pour cela est que 100 observations sur 100 individus fournissent plus d’informations que 100 observations sur 10 individus (ou grappes).

from lifelines import CoxPHFitter

df = pd.DataFrame({
    'T': [5, 3, 9, 8, 7, 4, 4, 3, 2, 5, 6, 7],
    'E': [1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0],
    'month': [10, 3, 9, 8, 7, 4, 4, 3, 2, 5, 6, 7],
    'age': [4, 3, 9, 8, 7, 4, 4, 3, 2, 5, 6, 7],
    'id': [1, 1, 1, 1, 2, 3, 3, 4, 4, 5, 6, 7]
})

cph = CoxPHFitter()
cph.fit(df, 'T', 'E', cluster_col='id')
cph.print_summary()

Pour plus d’exemples, consultez Corrélations entre les sujets dans un modèle de Cox.

Résidus

Après avoir ajusté un modèle de Cox, nous pouvons regarder en arrière et calculer des résidus de modèle importants. Ces résidus peuvent nous renseigner sur les non-linéarités non capturées, les violations des risques proportionnels et nous aider à répondre à d’autres questions de modélisation utiles. Voir Évaluation de l’ajustement du modèle de Cox à l’aide de résidus.

Modélisation de l’aléa de base et de la survie avec des modèles paramétriques

Normalement, le modèle Cox est semi-paramétrique, ce qui signifie que son danger de base, (h_0

Dans lignes de vie, il existe une option pour s’adapter à une ligne de base paramétrique avec 1) des splines cubiques ou 2) des dangers constants par morceaux. Les splines cubiques sont très flexibles et peuvent capturer les données sous-jacentes presque aussi bien que les méthodes non paramétriques, et avec beaucoup plus d’efficacité.

from lifelines.datasets import load_rossi
from lifelines import CoxPHFitter

rossi = load_rossi()

cph_spline = CoxPHFitter(baseline_estimation_method="spline", n_baseline_knots=5)
cph_spline.fit(rossi, 'week', event_col='arrest')

Pour accéder à l’aléa de base et à la survie de base, on peut utiliser baseline_hazard_ et baseline_survival_ respectivement. Une bonne chose à propos des modèles paramétriques est que nous pouvons également interpoler la survie / les dangers de base, voir baseline_hazard_at_times() et baseline_survival_at_times().

Ci-dessous, nous comparons les survies de base non paramétriques et entièrement paramétriques:

cph_semi = CoxPHFitter().fit(rossi, 'week', event_col='arrest')
cph_piecewise = CoxPHFitter(baseline_estimation_method="piecewise", breakpoints=[20, 35]).fit(rossi, 'week', event_col='arrest')

ax = cph_spline.baseline_cumulative_hazard_.plot()
cph_semi.baseline_cumulative_hazard_.plot(ax=ax, drawstyle="steps-post")
cph_piecewise.baseline_cumulative_hazard_.plot(ax=ax)

Modélisation de la survie de base avec des splines vs non paramétrique.

lignes de vie » Le modèle spline Cox peut également utiliser presque toutes les options non paramétriques, notamment: couches, pénalisant, chronologie, formule, etc.

Modèles de temps de panne accéléré

Supposons que nous ayons deux populations, A et B, avec des fonctions de survie différentes, (Sam)) et (S_B

Cela peut être interprété comme un ralentissement ou une accélération du déplacement le long de la fonction de survie. Un exemple classique de ceci est que les chiens vieillissent à 7 fois le taux des humains, c.-à-d. ( lambda = frac {1} {7} ). Ce modèle a d’autres propriétés intéressantes: le temps de survie moyen de la population B est ({ lambda} ) fois la durée moyenne de survie de la population A. De même avec le médian temps de survie.

Plus généralement, nous pouvons modéliser la ( lambda ) en fonction des covariables disponibles, soit:

Ce modèle peut accélérer ou ralentir les temps d’échec en fonction des covariables des sujets. Une autre caractéristique intéressante est la facilité d’interprétation des coefficients: une augmentation unitaire de (x_i ) signifie que le temps de survie moyen / médian change d’un facteur de ( exp (b_i) ).

Ensuite, nous choisissons une forme paramétrique pour la fonction de survie, (S

Nous appelons ces modèles de temps de défaillance accéléré, souvent raccourcis aux modèles AFT. En utilisant lignes de vie, nous pouvons adapter ce modèle (et l’inconnu ( rho ) paramètre aussi).

Le modèle Weibull AFT

Le modèle Weibull AFT est implémenté sous WeibullAFTFitter. L’API de la classe est similaire aux autres modèles de régression de lignes de vie. Après l’ajustement, les coefficients sont accessibles via params_ ou summary, ou bien imprimée en utilisant print_summary().

from lifelines import WeibullAFTFitter
from lifelines.datasets import load_rossi

rossi = load_rossi()

aft = WeibullAFTFitter()
aft.fit(rossi, duration_col='week', event_col='arrest')

aft.print_summary(3)  # access the results using aft.summary

"""

      duration col = 'week'
         event col = 'arrest'
number of subjects = 432
  number of events = 114
    log-likelihood = -679.917
  time fit was run = 2019-02-20 17:47:19 UTC

---
                     coef  exp(coef)  se(coef)      z       p  -log2(p)  lower 0.95  upper 0.95
lambda_ fin         0.272      1.313     0.138  1.973   0.049     4.365       0.002       0.543
        age         0.041      1.042     0.016  2.544   0.011     6.512       0.009       0.072
        race       -0.225      0.799     0.220 -1.021   0.307     1.703      -0.656       0.207
        wexp        0.107      1.112     0.152  0.703   0.482     1.053      -0.190       0.404
        mar         0.311      1.365     0.273  1.139   0.255     1.973      -0.224       0.847
        paro        0.059      1.061     0.140  0.421   0.674     0.570      -0.215       0.333
        prio       -0.066      0.936     0.021 -3.143   0.002     9.224      -0.107      -0.025
        Intercept  3.990     54.062     0.419  9.521 <0.0005    68.979       3.169       4.812
rho_    Intercept  0.339      1.404     0.089  3.809 <0.0005    12.808       0.165       0.514
---
Concordance = 0.640
AIC = 1377.833
log-likelihood ratio test = 33.416 on 7 df
-log2(p) of ll-ratio test = 15.462
"""

D’en haut, on peut voir que prio, qui est le nombre d’incarcérations antérieures, a un coefficient négatif élevé. Cela signifie que chaque addition d’incarcérations modifie le temps de survie moyen / médian d’un sujet de ( exp (-0,066) = 0,936 ), une diminution d’environ 7% du temps de survie moyen / médian. Quelle est la durée de survie moyenne / médiane?

print(aft.median_survival_time_)
print(aft.mean_survival_time_)

# 100.325
# 118.67

Que fait le rho_    _intercept ligne moyenne dans le tableau ci-dessus? En interne, nous modélisons le log du rho_ paramètre, donc la valeur de ( rho ) est l’exponentielle de la valeur, donc en cas supérieur, c’est ( hat { rho} = exp0.339 = 1.404 ). Cela nous amène au point suivant – la modélisation ( rho ) avec des covariables également:

Modélisation des paramètres auxiliaires

Dans le modèle ci-dessus, nous avons laissé le paramètre ( rho ) comme un seul inconnu. Nous pouvons également choisir de modéliser ce paramètre. Pourquoi voudrions-nous faire cela? Cela peut aider dans la prédiction de survie pour permettre l’hétérogénéité dans le ( rho ) paramètre. Le modèle n’est plus un modèle AFT, mais nous pouvons toujours récupérer et comprendre l’influence du changement d’une covariable en regardant son graphique de résultat (voir la section ci-dessous). Modeler ( rho ), nous utilisons le ancillary argument de mot-clé dans l’appel à fit(). Il existe quatre options valides:

1. False ou None: ne modélisez pas explicitement le rho_ paramètre (sauf pour son interception).
2. un Pandas DataFrame. Cette option utilisera les colonnes du Pandas DataFrame comme covariables dans la régression pour rho_. Ce DataFrame peut être un égal ou un sous-ensemble de l’ensemble de données d’origine utilisé pour la modélisation lambda_, ou il pourrait s’agir d’un ensemble de données totalement différent.
3. True. En passant True réutilisera en interne l’ensemble de données utilisé pour modéliser lambda_.
4. Une formule de type R.

aft = WeibullAFTFitter()

aft.fit(rossi, duration_col='week', event_col='arrest', ancillary=False)
# identical to aft.fit(rossi, duration_col='week', event_col='arrest', ancillary=None)


aft.fit(rossi, duration_col='week', event_col='arrest', ancillary=some_df)


aft.fit(rossi, duration_col='week', event_col='arrest', ancillary=True)
# identical to aft.fit(rossi, duration_col='week', event_col='arrest', ancillary=rossi)
# identical to aft.fit(rossi, duration_col='week', event_col='arrest', ancillary="fin + age + race + wexp + mar + paro + prio")

aft.print_summary()

"""

      duration col = 'week'
         event col = 'arrest'
number of subjects = 432
  number of events = 114
    log-likelihood = -669.40
  time fit was run = 2019-02-20 17:42:55 UTC

---
                    coef  exp(coef)  se(coef)     z      p  -log2(p)  lower 0.95  upper 0.95
lambda_ fin         0.24       1.28      0.15  1.60   0.11      3.18       -0.06        0.55
        age         0.10       1.10      0.03  3.43 <0.005     10.69        0.04        0.16
        race        0.07       1.07      0.19  0.36   0.72      0.48       -0.30        0.44
        wexp       -0.34       0.71      0.15 -2.22   0.03      5.26       -0.64       -0.04
        mar         0.26       1.30      0.30  0.86   0.39      1.35       -0.33        0.85
        paro        0.09       1.10      0.15  0.61   0.54      0.88       -0.21        0.39
        prio       -0.08       0.92      0.02 -4.24 <0.005     15.46       -0.12       -0.04
        Intercept  2.68      14.65      0.60  4.50 <0.005     17.14        1.51        3.85
rho_    fin        -0.01       0.99      0.15 -0.09   0.92      0.11       -0.31        0.29
        age        -0.05       0.95      0.02 -3.10 <0.005      9.01       -0.08       -0.02
        race       -0.46       0.63      0.25 -1.79   0.07      3.77       -0.95        0.04
        wexp        0.56       1.74      0.17  3.32 <0.005     10.13        0.23        0.88
        mar         0.10       1.10      0.27  0.36   0.72      0.47       -0.44        0.63
        paro        0.02       1.02      0.16  0.12   0.90      0.15       -0.29        0.33
        prio        0.03       1.03      0.02  1.44   0.15      2.73       -0.01        0.08
        Intercept  1.48       4.41      0.41  3.60 <0.005     11.62        0.68        2.29
---
Concordance = 0.63
Log-likelihood ratio test = 54.45 on 14 df, -log2(p)=19.83
"""

Traçage

L’API de traçage est la même que dans CoxPHFitter. Nous pouvons afficher toutes les covariables dans un graphique forestier:

from matplotlib import pyplot as plt

wft = WeibullAFTFitter().fit(rossi, 'week', 'arrest', ancillary=True)
wft.plot()

Nous pouvons observer l’influence d’une variable dans le modèle en traçant la résultat (c.-à-d. survie) de changer la variable. Ceci est fait en utilisant plot_partial_effects_on_outcome(), et c’est aussi un bon moment pour observer les effets de la modélisation rho_ vs le garder fixe. Ci-dessous, nous ajustons le modèle Weibull au même ensemble de données deux fois, mais dans le premier modèle, nous modélisons rho_ et dans le deuxième modèle, nous ne le faisons pas. Nous quand varier le prio (qui est le nombre d’arrestations antérieures) et observez comment la survie change.

fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(10, 4))

times = np.arange(0, 100)
wft_model_rho = WeibullAFTFitter().fit(rossi, 'week', 'arrest', ancillary=True, timeline=times)
wft_model_rho.plot_partial_effects_on_outcome('prio', range(0, 16, 3), cmap='coolwarm', ax=ax[0])
ax[0].set_title("Modelling rho_")

wft_not_model_rho = WeibullAFTFitter().fit(rossi, 'week', 'arrest', ancillary=False, timeline=times)
wft_not_model_rho.plot_partial_effects_on_outcome('prio', range(0, 16, 3), cmap='coolwarm', ax=ax[1])
ax[1].set_title("Not modelling rho_");

En comparant quelques-unes de ces fonctions de survie côte à côte, on peut voir que la modélisation rho_ produit un ensemble plus flexible (diversifié) de fonctions de survie.

fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(7, 4))

# modeling rho == solid line
wft_model_rho.plot_partial_effects_on_outcome('prio', range(0, 16, 5), cmap='coolwarm', ax=ax, lw=2, plot_baseline=False)

# not modeling rho == dashed line
wft_not_model_rho.plot_partial_effects_on_outcome('prio', range(0, 16, 5), cmap='coolwarm', ax=ax, ls='--', lw=2, plot_baseline=False)

ax.get_legend().remove()

Pour en savoir plus et voir d’autres exemples d’extensions à, consultez la documentation de plot_partial_effects_on_outcome()

Prédiction

Étant donné un nouveau sujet, nous aimerions poser des questions sur leur survie future. Quand sont-ils susceptibles de vivre l’événement? À quoi ressemble leur fonction de survie? La WeibullAFTFitter est en mesure d’y répondre. Si nous avons modélisé les covariables auxiliaires, nous devons également les inclure:

X = rossi.loc[:10]

aft.predict_cumulative_hazard(X, ancillary=X)
aft.predict_survival_function(X, ancillary=X)
aft.predict_median(X, ancillary=X)
aft.predict_percentile(X, p=0.9, ancillary=X)
aft.predict_expectation(X, ancillary=X)

Il existe deux hyper-paramètres qui peuvent être utilisés pour obtenir un meilleur score au test. Ceux-ci sont penalizer et l1_ratio dans l’appel à WeibullAFTFitter. Le pénaliseur est similaire à celui de scikit-learn ElasticNet modèle, voir leur docs. (Pourtant, lignes de vie acceptera également un tableau pour la valeur de pénalité personnalisée par variable, voir la documentation Cox ci-dessus)

aft_with_elastic_penalty = WeibullAFTFitter(penalizer=1e-4, l1_ratio=1.0)
aft_with_elastic_penalty.fit(rossi, 'week', 'arrest')
aft_with_elastic_penalty.predict_median(rossi)

aft_with_elastic_penalty.print_summary(columns=['coef', 'exp(coef)'])

"""

             duration col = 'week'
                event col = 'arrest'
                penalizer = 0.0001
   number of observations = 432
number of events observed = 114
           log-likelihood = -679.97
         time fit was run = 2020-08-09 15:04:35 UTC

---
                    coef  exp(coef)
param   covariate
lambda_ age         0.04       1.04
        fin         0.27       1.31
        mar         0.31       1.36
        paro        0.06       1.06
        prio       -0.07       0.94
        race       -0.22       0.80
        wexp        0.11       1.11
        Intercept   3.99      54.11
rho_    Intercept   0.34       1.40
---
Concordance = 0.64
AIC = 1377.93
log-likelihood ratio test = 33.31 on 7 df
-log2(p) of ll-ratio test = 15.40

"""

Les modèles AFT log-normal et log-logistique

Il y a aussi les LogNormalAFTFitter et LogLogisticAFTFitter modèles, qui au lieu de supposer que la distribution du temps de survie est de Weibull, nous supposons qu’elle est respectivement Log-Normal ou Log-Logistique. Ils ont des API identiques aux WeibullAFTFitter, mais les noms des paramètres sont différents.

from lifelines import LogLogisticAFTFitter
from lifelines import LogNormalAFTFitter

llf = LogLogisticAFTFitter().fit(rossi, 'week', 'arrest')
lnf = LogNormalAFTFitter().fit(rossi, 'week', 'arrest')

Plus de modèles AFT: modèle CRC et modèle gamma généralisé

Pour une flexibilité et lisse modèle paramétrique, il y a le GeneralizedGammaRegressionFitter. Ce modèle est en fait une généralisation de tous les modèles AFT ci-dessus (c’est-à-dire que les valeurs spécifiques de ses paramètres représentent un autre modèle) – voir la documentation pour les valeurs de paramètres spécifiques. L’API est cependant légèrement différente et ressemble plus à la façon dont les modèles de régression personnalisés sont construits (voir la section suivante sur Modèles de régression personnalisés).

from lifelines import GeneralizedGammaRegressionFitter
from lifelines.datasets import load_rossi

df = load_rossi()
df['Intercept'] = 1.

# this will regress df against all 3 parameters
ggf = GeneralizedGammaRegressionFitter(penalizer=1.).fit(df, 'week', 'arrest')
ggf.print_summary()

# If we want fine control over the parameters <-> covariates.
# The values in the dict become can be formulas, or column names in lists:
regressors = {
    'mu_': rossi.columns.difference(['arrest', 'week']),
    'sigma_': ["age", "Intercept"],
    'lambda_': 'age + 1',
}

ggf = GeneralizedGammaRegressionFitter(penalizer=0.0001).fit(df, 'week', 'arrest', regressors=regressors)
ggf.print_summary()

De même, il existe le modèle CRC qui utilise des splines pour modéliser le temps. Voir un article de blog à ce sujet Ici.

Les modèles de régression exponentielle par morceaux

Une autre classe de modèles paramétriques implique une modélisation plus flexible de la fonction de risque. La PiecewiseExponentialRegressionFitter peut modéliser les sauts dans l’aléa (pensez: les différences de «survie-de-rester à l’école» entre les élèves de 1re, 2e année, 3e année et 4e année) et les valeurs constantes entre les sauts. La possibilité de spécifier quand ces sauts se produisent, appelés points d’arrêt, offre aux modélisateurs une grande flexibilité. Un exemple d’application impliquant le désabonnement des clients est disponible dans ce carnet.

AIC et sélection de modèles pour les modèles paramétriques

Souvent, vous ne savez pas a priori quel modèle paramétrique utiliser. Chaque modèle a des hypothèses intégrées (pas encore implémentées dans lignes de vie), mais une méthode rapide et efficace consiste à comparer les AIC pour chaque modèle ajusté. (Dans ce cas, le nombre de paramètres pour chaque modèle est le même, il s’agit donc vraiment de comparer la log-vraisemblance). Le modèle avec le plus petit AIC fait le meilleur travail d’adaptation aux données avec un minimum de degrés de liberté.

from lifelines import LogLogisticAFTFitter, WeibullAFTFitter, LogNormalAFTFitter
from lifelines.datasets import load_rossi

rossi = load_rossi()

llf = LogLogisticAFTFitter().fit(rossi, 'week', 'arrest')
lnf = LogNormalAFTFitter().fit(rossi, 'week', 'arrest')
wf = WeibullAFTFitter().fit(rossi, 'week', 'arrest')

print(llf.AIC_)  # 1377.877
print(lnf.AIC_)  # 1384.469
print(wf.AIC_)   # 1377.833, slightly the best model.


# with some heterogeneity in the ancillary parameters
ancillary = rossi[['prio']]
llf = LogLogisticAFTFitter().fit(rossi, 'week', 'arrest', ancillary=ancillary)
lnf = LogNormalAFTFitter().fit(rossi, 'week', 'arrest', ancillary=ancillary)
wf = WeibullAFTFitter().fit(rossi, 'week', 'arrest', ancillary=ancillary)

print(llf.AIC_) # 1377.89, the best model here, but not the overall best.
print(lnf.AIC_) # 1380.79
print(wf.AIC_)  # 1379.21

Données censurées à gauche, à droite et à intervalles

Les modèles paramétriques ont des API qui gèrent également les données censurées à gauche et à intervalle. L’API pour eux est différente de l’API pour l’ajustement aux données censurées à droite. Voici un exemple avec des données censurées par intervalle.

from lifelines.datasets import load_diabetes

df = load_diabetes()
df['gender'] = df['gender'] == 'male'

print(df.head())
"""
   left  right  gender
1    24     27    True
2    22     22   False
3    37     39    True
4    20     20    True
5     1     16    True
"""

wf = WeibullAFTFitter().fit_interval_censoring(df, lower_bound_col='left', upper_bound_col='right')
wf.print_summary()

"""

          lower bound col = 'left'
          upper bound col = 'right'
                event col = 'E_lifelines_added'
   number of observations = 731
number of events observed = 595
           log-likelihood = -2027.20
         time fit was run = 2020-08-09 15:05:09 UTC

---
                    coef  exp(coef)   se(coef)   coef lower 95%   coef upper 95%  exp(coef) lower 95%  exp(coef) upper 95%
param   covariate
lambda_ gender      0.05       1.05       0.03            -0.01             0.10                 0.99                 1.10
        Intercept   2.91      18.32       0.02             2.86             2.95                17.53                19.14
rho_    Intercept   1.04       2.83       0.03             0.98             1.09                 2.67                 2.99
                       z      p   -log2(p)
param   covariate
lambda_ gender      1.66   0.10       3.38
        Intercept 130.15 <0.005        inf
rho_    Intercept  36.91 <0.005     988.46
---
AIC = 4060.39
log-likelihood ratio test = 2.74 on 1 df
-log2(p) of ll-ratio test = 3.35
"""

Un autre exemple d’utilisation de lignes de vie pour les données censurées par intervalle est situé Ici.

Modèles de régression paramétrique personnalisés

lignes de vie a une syntaxe très générale pour créer vos propres modèles de régression paramétrique. Si vous souhaitez créer vos propres modèles personnalisés, consultez la documentation Modèles de régression personnalisés.