L’Analyse de survie

add.png

\bullet Introduction:

L’Analyse de survie regroupe un ensemble d’outils statistique permettant de construire des indicateurs sur la durée de vie représentée par une variable aléatoire T dans le but de quantifier le délai de survenu d’un évènement d’un instant à l’autre. Une autre façon homogène de positionner le problème est de considérer l’analyse de survie comme une approche visant à étudier une variable réponse Y binaire par rapport à une matrice de P + 1 variables explicatives (T, \mathbf{X})T déjà définie auparavant et \mathbf{X} continues et/ou qualitatives.

Universellement, nous notons S la fonction de survie à modéliser et de formule,

S(t) = P(T > t)

La notion de données censurées aura également un impact important dans une analyse de survie. Une donnée censurée (ou censure) est une unité statistique qui aura été perdue de vue lors de la période d’observation mais dont l’information récoltée avant sa disparition sera prise en compte dans le modèle en étant considérée comme n’ayant ni connu ni pas connu l’évènement à mesurer.

Le principe de l’analyse de survie est principalement de modéliser une courbe de survie. Trois types d’approche existent: non paramétrique, semi-paramétrique et paramétrique.

Enfin, rappelons que le terme de survie n’est pas à prendre au sens stricte du terme puisque que l’analyse de survie s’étend également aux phénomènes d’apparition d’un évènement. La raison étant que l’analyse de survie avait été initialement développée pour mesurer le temps survenu avant un décès et a ensuite été extrapolée à d’autres thématiques de recherche.

\bullet Estimateur de Kaplan-Meier:

Présentation:

Publié en 1958 par Edward L. Kaplan et Paul Meier, l’estimateur de Kaplan-Meier, également connu sous le nom d’estimateur produit-limite, est un estimateur non paramétrique du maximum de vraisemblance visant à estimer la fonction de survie.

L’estimateur:

Pour un temps t fixé, la formule de l’estimateur de Kaplan-Meier est,

\hat{S} (t) = \prod_{t_i \leq t} \frac{n_i - Survie_i}{n_i}

Où,

n_i représente l’effectif restant pour t \geq t_i,

Survie_i = \sharp \lbrace \mbox{Evenement(s) produit(s) a } t_i = t \rbrace.

Gestion des ex-aequos:

Pour un même temps, il peut y avoir plusieurs situations différentes (apparition de l’évènement ou non). Lorsque cette situation se produit, la démarche à suivre est de considérer les cas où l’évènement a eu lieu avant ceux où il n’a pas eu lieu.

Autres estimateurs non paramétriques:

Plusieurs estimateurs existent dans cet objectif tel que celui de Nelson-Arlen, de Breslow, de Harrington-Fleming et de la méthode actuarielle. Il n’y pas vraiment de meilleur estimateur non paramétrique que les autres, disons simplement que l’estimateur de Kaplan-Meier est le plus populaire d’entre eux et donc le plus utilisé.

\bullet La régression de Cox:

Présentation:

La régression de Cox a été publiée en 1975 par David Cox dans l’objectif de modéliser des données de survies. Elle porte également le nom de modèle de Cox ou de modèle à risque proportionnel et est considérée comme une approche semi-paramétrique pour la modélisation de la courbe de survie.

Le modèle:

La régression de Cox consiste à modéliser la probabilité que l’évènement se produise au temps t alors qu’il ne s’était pas produit au temps t-1. Le modèle s’écrit,

\lambda (t, \mathbf{X}_i) = \lambda_0 (t) e ^{\sum_{p = 1} ^P \beta_p X_i ^p}

Avec \lambda_0 (t) le risque instantané de survenu de l’évènement lorsque (X_1, \cdots, X_n) sont nulles. L’Aspect semi-paramétrique de la régression de Cox provient du fait qu’il ne requiert pas l’estimation du paramètre \lambda_0(t), finalement indépendant du temps t.

L’estimation des coefficients \beta:

Posons n_E = \sharp \lbrace \mbox{Evenements produits} \rbrace. Les paramètres du modèle de Cox se déterminent de manière classique par maximisation de la vraisemblance et de formule,

L(\beta) = \prod_{i = 1} ^{n_E} \frac{e ^{\sum_{p = 1} ^P \beta_p X_i ^p}}{\sum_{i' \in [i,n_E]} e ^{\sum_{p = 1} ^P \beta_p X_i ^p}}

La résolution se fait, par exemple, via l’algorithme de Newton-Raphson,

– Etape intiale \beta = \beta_0 = (0, \cdots, 0)

– Etape 1: Calcul de la dérivée partiel de la vraisemblance,

U(\beta) = \frac{\partial}{\partial \beta} L(\beta) = (\frac{\partial}{\partial \beta_1} L(\beta), \cdots, \frac{\partial}{\partial \beta_P} L(\beta))

\Rightarrow U(\beta) = (\sum_{i = 1} ^{n_E} [X_i ^1 - \frac{\sum_{i' \in [i,n]} X_i ^1 e ^{\sum_{p = 1} ^P \beta_p X_{i'} ^p}}{\sum_{i' \in [i,n]} e ^{\sum_{p = 1} ^P \beta_p X_{i'} ^p}}], \cdots, \sum_{i = 1} ^{n_E} [X_i ^P - \frac{\sum_{i' \in [i,n]} X_i ^P e ^{\sum_{p = 1} ^P \beta_p X_{i'} ^p}}{\sum_{i' \in [i,n]} e ^{\sum_{p = 1} ^P \beta_p X_{i'} ^p}}])

– Etape 2: Calcul de la matrice jacobienne des dérivées secondes de la vraisemblance: I(\beta). Les éléments de la diagonale sont de formules,

\frac{\partial}{\partial \beta_p ^2} L(\beta) = - \frac{\sum_{i' \in [i,n]} e ^{\sum_{p = 1} ^P \beta_p X_{i'} ^p} \sum_{i' \in [i,n]} X_{i'} ^p e ^{\sum_{p = 1} ^P \beta_p X_{i'} ^p} - (\sum_{i' \in [i,n]} X_{i'} ^p e ^{\sum_{p = 1} ^P \beta_p X_{i'} ^p}) ^2}{(\sum_{i' \in [i,n]} e ^{\sum_{p = 1} ^P \beta_p X_{i'} ^p}) ^2}

, et ceux non diagonaux de formules,

\frac{\partial}{\partial \beta_{p_1} \partial \beta_{p_2}} L(\beta) = - \frac{\sum_{i' \in [i,n]} e ^{\sum_{p = 1} ^P \beta_p X_{i'} ^p} \sum_{i' \in [i,n]} X_{i'} ^{p_1} X_{i'} ^{p_2} e ^{\sum_{p = 1} ^P \beta_p X_{i'} ^p} - \sum_{i' \in [i,n]} X_{i'} ^{p_1} e ^{\sum_{p = 1} ^P \beta_p X_{i'} ^p} \sum_{i' \in [i,n]} X_{i'} ^{p_2} e ^{\sum_{p = 1} ^P \beta_p X_{i'} ^p}}{(\sum_{i' \in [i,n]} e ^{\sum_{p = 1} ^P \beta_p X_{i'} ^p}) ^2}

– Etape 3: a l’itération k+1,

\beta ^{k+1} = \beta - U(\beta) \cdot I(\beta) ^{-1}

– Etape finale: Si || \beta ^{k+1} - \beta ^k || > \epsilon alors nous reprenons les étapes 1, 2, 3 avec \beta = \beta ^{k+1}. Sinon, l’algorithme s’arrête et nous avons notre vecteur de coefficients finaux.

Les hypothèses de validité:

Le modèle de Cox repose sur l’hypothèse forte des risques proportionnels et nécessite qu’elle soit vérifiée avant application. Pour cela il s’agit de vérifier que le temps n’a aucune influence bénéfique, nocif ou nul sur les différentes variables de \mathbf{X}.

Le risque relatif ajusté ou hazard ratio (RR):

Du modèle de Cox, il est possible de déduire le risque relatif associé à chaque variable. Deux formats sont à considérer,

– si Y est binaire ou continue alors,

RR_p = e ^{\beta_p}, p \in [1,P]

– si Y est ordinale, à valeurs dans [a, b], alors,

RR_p = e ^{\beta_p (b - a)}, p \in [1,P]

Tests associés au modèle de Cox:

Trois principaux tests existent pour juger de la fiabilité des coefficients,

– Le test du rapport de vraisemblance,

\chi_V ^2 = 2 [log L(\beta) - log(\beta_0)]

– Le test de Wald,

\chi_W ^2 = (\beta - \beta_0) I(\beta) (\beta - \beta_0)

– Le test du Score,

\chi_S ^2 = (U(\beta)) ^T (I(\beta_0)) ^{-1} U(\beta_0)

Ces trois tests suivent une loi du \chi ^2 à P degrés de liberté et l’hypothèse H_0 est \beta = \beta_0 soit que le modèle considéré avec l’information auxiliaire \mathbf{X} est aussi performant que sans.

Gestion des ex-aequos:

Dans le cas de survenue de plusieurs évènements pour un même temps, le calcul des dérivées partielles de la vraisemblance doit être modifié. Plusieurs méthodes existent, parmi les plus populaires nous avons la méthode de Breslow et celle d’Efron.

\bullet Approches paramétriques:

Une approche purement paramétrique est également possible. Il s’agit de partir du constat que la courbe de survie est identifiable à l’une des nombreuses loi de probabilité dont les principales sont,

– La loi Exponentielle,

S(t | \theta) = e ^{-\frac{t}{\theta}}

– La loi de Weibull,

S(t | \theta, \alpha) = e ^{-(\frac{t}{\theta}) ^{\alpha}}

– La loi Gamma,

S(t | \theta, \alpha) = \frac{\frac{\alpha^{\theta}}{\Gamma(\theta)} t ^{\theta - 1} e ^{\alpha t}}{1 - \frac{1}{\Gamma(\theta)} \int_0 ^{\alpha t} u ^{\theta - 1} e ^{-u} du}

\bullet Le test du log-rank:

Présentation:

Publié par Nathan Mantel et David Cox en 1966, le test du log-rank (connu aussi sous le nom de test de Mantel-Cox) est une approche non paramétrique permettant de comparer les distributions de deux séries de survie.

Le test:

Le calcul de la statistique de test peut se résumer en trois étapes itératives. Soit T l’ensemble des temps considéré pour les deux séries considérées.

– Sortie des temps pour lesquelles l’évènement que nous cherchons à mesurer n’est pas survenu au sein des deux séries d’intérêt,

T_S = T_S (serie_1) \cup T_S (serie_2)

– Pour chacun des temps retenus dans T_S = (T_S ^1, \cdots, T_S ^t, \cdots), construire le tableau,

add2

– Calcul des effectifs observés (O_t), théoriques (E_t) et de la variance (V_t) associé au tableau ci-dessus,

O_t = C

E_t = (A + B) \frac{C + D}{A + B + C + D}

V_t = \frac{(A + B) (C + D) (A + C) [(A + B + C + D) - (A + C)]}{(A + B + C + D) ^2 [(A + B + C + D) - 1]}

A noter que si V_t \rightarrow \infty alors .</span></p> <p style="text-align:justify;"><span style="color:#808080;">– Une fois les différents temps parcourus, calculer la statistique de test du log-rank,</span></p> <p style="text-align:center;"><span style="color:#808080;">[latex]Z = |\frac{(\sum_t O_t - \sum_t E_t) ^2}{\sum_t V_t}|

Elle suit une loi du \chi ^2 à un degré de liberté et l’hypothèse H_0 est: « les deux courbes de survie sont différentes ».

Autres tests pour la comparaison de deux courbes de survie:

D’autres tests existent dans cet objectif, parmi les plus populaires nous avons le test du log-rank pondéré, le test de Gehan et le test de Peto et Prentice.

Tous partagent la même condition que le test du log-rank, les courbes de survies ne doivent pas se croiser sous peine d’occasionner une forte chute de la puissance statistique.

\bullet Exemple:

Soit l’échantillon suivant,

add

, où la variable « Temps » est le moment où nous observons la production de l’évènement mesuré (« Survie = 0« ) ou non (« Survie = 1« ).

Estimateur de Kaplan-Meier:

Pour le groupe A,

– Nous avons pour « Temps = 59« ,

\hat{S} (t = 59) = \frac{6 - 1}{6} = 0.83333333333

– Puis pour « Temps \in [59, 115]« ,

\hat{S} (t \in [59, 115]) = \frac{5- 1}{5} \times \hat{S} (t = 59) = 0.8 \times 0.83333333333 = 0.666666667

– Puis pour « Temps \in [59, 115, 156]« ,

\hat{S} (t \in [59, 115, 156]) = \frac{4 - 1}{4} \times \hat{S} (t \in [59,115]) = 0.75 \times 0.666666667 = 0.5

– Puis pour « Temps \in [59, 115, 156, 431]« ,

\hat{S} (t \in [59, 115, 156, 431]) = \frac{3 - 1}{3} \times \hat{S} (t \in [59,115,156]) = 0.333333333

– Puis pour « Temps \in [59, 115, 156, 431, 448]« ,

\hat{S} (t \in [59, 115, 156, 431, 448]) = \frac{2 - 0}{2} \times \hat{S} (t \in [59, 115, 156, 431]) = 0.333333333

– Puis pour « Temps \in [59, 115, 156, 431, 448, 477]« ,

\hat{S} (t \in [59, 115, 156, 431, 448, 477]) = \frac{1 - 0}{1} \times \hat{S} (t \in [59, 115, 156, 431, 448])

= 0.333333333

En procédant aux mêmes calculs pour le groupe B, nous pouvons alors les deux courbes de Kaplan-Meier suivantes,

add

Test du Log-Rank:

La courbe ci-dessus semble laisser penser qu’il y a de grandes différences entre le groupe A et le groupe B. Vérifions le au sens statistique du terme.

Dans un premier temps, relevons les temps pour lesquelles l’évènement n’a pas lieu (« survie = 1« ). Nous avons pour les groupes A et B:

(T_S)_A = (59,115,156,431) et (T_S)_B = (464,475,563)

Quatre cas pour le groupe A et trois cas pour le groupe B. Nous obtenons alors les temps fusionnés suivants,

T_S = (59,115,156,431,464,475,563)

Procédons itérativement pour l’ensemble des temps considérés (évènement produit ou non),

– pour T = 59,

add2

\Rightarrow O_{59} = 0

\Rightarrow E_{59} = (1 + 0) \times \frac{0 + 4}{1 + 5 + 0 + 4} = 1 \times \frac{4}{10} = 0.4

\Rightarrow V_{59} = \frac{(1 + 5) \times (0 + 4) \times (1 + 0) \times ((1 + 5 + 0 + 4) - (1 + 0))}{(1 + 5 + 0 + 4) ^2 \times ((1 + 5 + 0 + 4) - 1} = \frac{6 \times 4 \times 1 \times 9}{10 ^2 \times 9} = \frac{216}{900} = 0.24

– pour T = 115,

add2

\Rightarrow O_{115} = 0

\Rightarrow E_{115} = (1 + 0) \times \frac{0 + 4}{1 + 4 + 0 + 4} = 1 \times \frac{4}{9} = 0.444444444

\Rightarrow V_{115} = \frac{(1 + 4) \times (0 + 4) \times (1 + 0) \times ((1 + 4 + 0 + 4) - (1 + 0))}{(1 + 4 + 0 + 4) ^2 \times ((1 + 4 + 0 + 4) - 1} = \frac{5 \times 4 \times 1 \times 8}{9 ^2 \times 8} = \frac{160}{648} = 0.2469136

\cdots

– pour T = 563,

add2

\Rightarrow O_{563} = 1

\Rightarrow E_{563} = (0 + 1) \times \frac{1 + 0}{0 + 0 + 1 + 0} = 1 \times \frac{1}{1} = 1

\Rightarrow V_{563} = \frac{(0 + 0) \times (1 + 0) \times (0 + 1) \times ((0 + 0 + 1 + 0) - (0 + 1))}{(0 + 0 + 1 + 0) ^2 \times ((0 + 0 + 1 + 0) - 1} = \frac{0 \times 1 \times 1 \times 0}{1 ^2 \times 0} = 0

Nous pouvons alors calculer la statistique du test du Log-Rank,

O = \sum_{t \in [59,115,156,431,464,475,563]} O_t =0 + 0 + 0 + 0 + 1 + 1 + 1 = 3

E = \sum_{t \in [59,115,156,431,464,475,563]} E_t

= 0.4 + 0.444444444 + 0.5 + 0.5 + 0.75 + 0.6666667 + 1

= 4.261111

V = \sum_{t \in [59,115,156,431,464,475,563]} V_t

= 0.24 + 0.2469136 + 0.25 + 0.25 + 0.1875 + 0.2222222 + 0

= 1.138737

Par conséquent,

|Z| = |\frac{(3 - 4.261111) ^2}{1.396636}| = 1.138737

, que nous reportons à la loi du \chi ^2 pour un degré de liberté et nous obtenons p = 0.2926337. Nous ne pouvons rejeter H_0 et donc nous en concluons que les deux courbes ne sont pas identiques.

Régression de Cox:

Passons maintenant à la modélisation de l’évènement « Survie » pour les groupes A et B fusionnés. Nous initialisons l’algorithme de Newton-Raphson à \beta ^0 = 0.

Nous avons pour la première itération,

U(\beta ^0) = \sum_{i = 1} ^{10} Survie_i - \sum_{i = 1} ^{10} \frac{\sum_{i' \in [i,10]} Survie_{i'} e ^{0 \times Survie_{i'}}}{\sum_{i' \in [i,10]} e ^{0 \times Survie_{i'}}}

= 7 - (\frac{7}{10} + \frac{6}{9} + \cdots + \frac{1}{1})

= 7 - 6.746429

= 0.253571

I(\beta ^0) = - \sum_{i = 1} ^{10} \frac{\sum_{i' \in [i,10]} e ^{0 \times Survie_{i'}}\times \sum_{i' \in [i,10]} Survie_{i'} ^2 e ^{0 \times Survie_{i'}} - (\sum_{i' \in [i,10]} Survie_{i'} e ^{0 \times Survie_{i'}}) ^2}{(\sum_{i' \in [1,10]} e ^{0 \times Survie_{i'}}) ^2}

= - \frac{70 - 49}{100} + \frac{54 - 36}{81} + \cdots + \frac{1 - 1}{1}

= - 2.03344

Nous avons alors,

\beta ^1 = 0 - \frac{0.253571}{- 2.03344} = 0.1247005

Pour la seconde itération et en considérant \beta ^1 à la place de \beta ^0, nous obtenons,

U(\beta ^1) = 0.004441066

I(\beta ^1) = -1.960438

\beta ^2 = 0.1269661

Pour la troisième itération et en considérant \beta ^2 à la place de \beta ^1, nous obtenons,

U(\beta ^2) = 0.000001609696

I(\beta ^2) = -1.959017

\beta ^3 = 0.1269669

Pour la quatrième itération et en considérant \beta ^3 à la place de \beta ^2, nous obtenons,

U(\beta ^2) = 0.0000000000002113865

I(\beta ^2) = -1.959016

\beta ^4 = 0.1269669

Nous avons \beta ^4 = \beta ^3, nous en déduisons que l’agorithme de Newton-Raphson a convergé. Le coefficient associé à la variable « Survie » est donc \beta = 0.1269669.

\bullet Application informatique:

– Procédure SAS: https://support.sas.com/documentation/cdl/en/statug/63347/HTML/default/viewer.htm#statug_lifetest_sect004.htm

– Package R: https://stat.ethz.ch/R-manual/R-devel/library/survival/html/survfit.html

\bullet Bibliographie: 

– Comprendre et utiliser les statistiques dans les sciences de la vie de Bruno Falissard

– Nonparametric estimation from incomplete observations de Edward L. Kaplan et Paul Meier

– Evaluation of survival data and two new rank order de Nathan Mantel

– Asymptotically Efficient Rank Invariant Test Procedures de Richard Peto

– Linear Rank Tests in Survival Analysis de David Harrington

– Introduction à l’analyse de survie de Philippe Saint Pierre

– Partial Likelihood de David R. Cox