8 Le modèle de Cox

On peut ignorer la partie sur l’estimation du modèle. On retiendra tout de même qu’il est déconseillé d’utiliser la méthode dite exacte pour la correction de la vraisemblance, qui ne peut matériellement fonctionner qu’avec un nombre très limité d’évènements observés simultanément. Ce qui est plutôt rare avec des données à durées discrètes ou groupées, très fréquentes dans les sciences sociales.

8.1 Le modèle semi-paramétrique de Cox

8.1.1 La vraisemblance partielle et estimation des paramètres

On se situe dans une situation où la durée est mesurée sur une échelle strictement continue. Il ne peut donc y avoir qu’un seul évènement observé en $t_i$ (idem pour la censure).

On peut représenter le processus aléatoire d’une analyse de survie en présence de censure à droite, avec l’équation de vraisemblance suivante:

\[L_i=f(t_i)^{\delta_i}S(t_i)^{1-\delta_i}\]

$f(t_i)$ est la valeur de la fonction de densité en $t_i$
$S(t_i)$ est la valeur de la fonction de survie en $t_i$
$\delta_i=1$ si l’évènement est observé: $L_i=f(t_i)$
$\delta_i=0$ si l’observation est censurée: $L_i=S(t_i)$

Vraisemblance partielle de Cox

Comme $f(t_i)=h(t_i)\times S(t_i)$ ¹, on obtient: $L_i=[h(t_i)S(t_i)]^{\delta_i}S(t_i)^{1-\delta_i} = h(t_i)^{\delta_i}S(t_i)$.

Pour $i=1,2,.....,n$, la vraisemblance s’ecrit donc: $L_i=\prod_{i=1}^{n}h(t_i)^{\delta_i}S(t_i)$.

On peut réécrire cette vraisemblance en la multipliant et en la divisant par: $\sum_{j\in R_i}h(t_i)$, où $j\in R_i$ est l’ensemble des observations soumises au risque en $t_i$.

\[L=\prod_{i=1}^{n}\left[h(t_i)\frac{\sum_{j\in R}h(t_i)}{\sum_{j\in R}h(t_i)}\right]^{\delta_i}S(t_i)= \prod_{i=1}^{n}\left[\frac{h(t_i)}{\sum_{j\in R_i}h(t_i)}\right]^{\delta_i}\sum_{j\in R_i}h(t_i)^{\delta_i}S(t_i)\]

La vraisemblance partielle retient seulement le premier terme de la vraisemblance, soit:

\[PL=\prod_{i=1}^{n}\left[\frac{h(t_i)}{\sum_{j\in R}h(t_i)}\right]^{\delta_i}\]

Une fois remplacée la valeur de $h(t_i)$ par son expression en tant que modèle à risques proportionnels, la vraisemblance partielle ne dépendra plus de la durée. Mais elle va dépendre de l’ordre d’arrivée des évènements, c’est à dire leur rang.

Remarque: pour les observations censurées($\delta_i=0$), $PL=1$. Toutefois, ces censures à droite entrent dans l’expression $\sum_{j\in R}h(t_i)$ tant qu’elles sont soumises au risque.

En remplaçant donc $h(t_i)$ par l’expression $h_0(t)e^{X_i^{'}b}$:

\[PL=\prod_{i=1}^{n}\left[\frac{h_0(t)e^{X_{i}^{'}b}}{\sum_{j\in R_i}h_0(t)e^{X_{j}^{'}b}}\right]^{\delta_i} =\prod_{i=1}^{n}\left[\frac{e^{X_i^{'}b}}{\sum_{j\in R_i}e^{X_{j}^{'}b}}\right]^{\delta_i}\]

L’expression $\frac{e^{Xb}}{\sum_{j\in R}e^{Xb}}$ est donc bien une probabilité, et la vraisemblance partielle est donc bien un produit de probabilités. Pour un individu ayant connu l’évènement, la contribution à la vraisemblance partielle est la probabilité que l’individu observe l’évènement en $t_i$ sachant qu’un évènement (et un seul) s’est produit.

Si $\delta_i = 0$: $PL_i = 1$
Si $\delta_i = 1$: $PL_i =\frac{e^{X_i^{'}b}}{\sum_{j\in R_i}e^{X_{j}^{'}b}}$

Condition nécessaire: pas d’évènement simultané: en présence d’évènements mesurés simultanément, l’estimation de la vraisemblance doit faire l’objet d’une correction.

Correction de la vraisemblance avec des évènements simultanés:

La méthode dite exacte: Comme il ne doit pas y avoir d’évènement simultané, on va introduire à la vraisemblance partielle toutes les permutations possibles des évènements observés au même moment. Bien qu’en $t_i$ on observe au même moment l’évènement pour 2 observations (A,B) une métrique temporelle plus précise permettrait de savoir si A s’est produit avant B ou B s’est produit avant A (2 permutations). Comme le nombre de permutations est calculé à l’aide d’une factorielle ², avec 3 évènements mesurés simultanément, on obtient 6 permutations ($3\times2\times1$). Problème: le nombre de permutations pour chaque $t_i$ peut devenir très vite particulièrement élevé. Par exemple pour 10 évènements simultanés, le nombre de permutations est égal à 3,628,800. Le temps de calcul devient extrêmement long, et ce type de correction totalement inopérant.
La méthode dite de Breslow: il s’agit d’une approximation de la méthode exacte permettant de ne pas avoir à intégrer chaque permutation. Cette approximation est utilisée par défaut par les logiciels Sas et Stata.
La méthode dite d’Efron: elle corrige l’approximation de Breslow, et est jugée plus proche de la méthode exacte. C’est la méthode utilisée par défaut avec le logiciel R, et elle est disponible avec les autres applications.

8.1.2 Estimation des paramètres

On utilise la méthode habituelle, à savoir la maximisation de la log-vraisemblance (ici partielle).

Conditions de premier ordre: calcul des équations de score à partir des dérivées partielles. Solution: $\frac{\partial log(PL)}{\partial{b_k}}=0$. On ne peut pas obtenir de solution numérique directe.

Remarque: les équations de score sont utilisées pour tester la validité de l’hypothèse de constance des rapports de risque pour calculer les résidus de Schoenfeld (voir plus loin).

Conditions de second ordre: calcul des dérivées secondes qui permettent d’obtenir la matrice d’information de Fisher et la matrice des variances-covariances des paramètres.
Comme il n’y a pas de solution numérique directe, on utilise un algorithme d’optimisation (ex: Newton-Raphson) à partir des équations de score et de la matrice d’information de Fisher.

Eléments de calcul

En logarithme (sans évènement simultané), la vraisemblance partielle s’ecrit:

\[pl(b)=\sum_{i=1}^n\delta_i\left(log(e^{X_{i}^{'}b})-log\sum_{j\in R_i}e^{X_{j}^{'}b}\right)\]

\[pl(b)=\sum_{i=1}^n\delta_i\left(X_{i}^{'}b-log\sum_{j\in R_i}e^{X_{j}^{'}b}\right)\]

Calcul de l’équation de score pour une covariable $X_k$:

\[\frac{\partial pl(b)}{\partial{b_k}}=\sum_{i=1}^n\delta_i\left(X_{ik}-\sum_{j\in R_i}X_{jk}\frac{e^{X_{j}^{'}b}}{\sum_{j\in R_i}e^{X_{j}^{'}b}}\right)\]

Comme $\frac{e^{X_{j}b}}{\sum_{j\in R}e^{X_{j}b}}$ est une probabilité, et $\sum_{j\in R}X_{ik}\times p_i$ est l’espérance (la moyenne) $E(X_k)$ d’avoir la caractéristique $X_k$ lorsqu’un évènement a été observé. Au final:

\[\frac{\partial lp(b)}{\partial{b_k}}= \sum_{i=1}^n\delta_i\left(X_{ik} - E(X_{j\in R_i,k}) \right)\]

Cette expression va permettre d’analyser le respect ou non de l’hypothèse de risques proportionnels via les résidus de Schoenfeld.

8.1.3 Lecture des résultats

Comme il s’agit d’un modèle à risque proportionnel, les rapports de risques sont constants pendant toute la période d’observation. Il s’agit d’une propriété de l’estimation.

Covariable binaire (indicatrice) $X=(0,1)$: $RR=\frac{h(t\ |\ X=1)}{h(t\ |\ X=0)}=e^b$.
A chaque moment de la durée $t$, le risque d’observer l’évènement est $e^b$ fois plus important/plus faible pour $X=1$ que pour $X=0$.

Covariable quantitative (fixe dans le temps)

$RR=\frac{h(t\ |\ X=a+c)}{h(t\ |\ X=a)}=e^{c \times{b}}$. On prendra pour illustrer une variable type âge au début de l’exposition au risque (a) et un delta de comparaison avec un âge inférieur c.
Si $c=1$ (résultat de l’estimation): A un âge donnée, le risque de connaitre l’évènement est $e^b$ fois inférieur/supérieur à celui d’une personne qui a un an de moins.

Exemple pour les insuffisances cardiaques

Correction de la vraisemblance: méthode d’Efron
Nombre d’observations: 103
Nombre de décès: 75
Log-Vraisemblance: -289.30639

Cox: log Hazard Ratio (Risks Ratio)
Variables	logRR	Std.Err	z	$P>\|z\|$	95% IC
year	-0.119	0.0673	-1.78	0.076	-0.2516;+0.0124
age	+0.0296	0.0135	2.19	0.029	+0.0031;+0.0561
surgery	-0.9873	0.4363	-2.26	0.024	-1.8424;-0.1323

Cox: Hazard Ratio (Risks Ratio)
Variables	RR	Std.Err	z	$P>\|z\|$	95%CI
year	0.8872	0.0597	-1.78	0.076	0.7775; 1.0124
age	1.0300	0.0139	2.19	0.029	1.0031; 1.0577
surgery	0.3726	0.1625	-2.26	0.024	0.1584; 0.8761

On retrouve les des tests non paramétriques pour l’opération, à savoir qu’un pontage réduit les risques journaliers de décès pendant la période d’observation (augmente la durée de survie).
De la même manière, plus on entre à un âge élevé dans la liste d’attente plus le risque de décès augmente. La variable year, qui traduit des progrès en médecine, renvoie à une réduction plutôt modérée du risque journalier de décès durant l’attente d’une greffe.

R-Stata-Sas-Python

Le modèle est estimé avec la fonction coxph de la librairie survival. Hors options, la syntaxe est identiques aux fonctions survfit et survdif.

8.1.3.1 Stata

Le modèle est estimé avec la commande stcox.

8.1.3.2 SAS

Le modèle est estimé avec la proc phreg.

8.1.3.3 Python

Avec la librairie lifelines, le modèle est estimé avec la fonction CoxPHFitter. Avec la librairie statmodels, il est estimé avec la fonction smf.phreg.

8.2 Analyse de la constance des rapports de risque

Les rapports de risque (RR) estimés par le modèle sont contraints à être constant sur toute la période d’observation. C’est une hypothèse forte.
Le respect de cette hypothèse doit être analysé, en particulier pour le modèle de Cox où la baseline du risque est habituellement estimée à l’aide de ces rapports (par exemple la méthode dite de Breslow, non traitée). En post-estimation, les valeurs estimées du risque pourront présenter des valeurs aberrantes si on dévie trop de constance, en particulier en obtenant des négatives des taux de risque.
Analyser cette hypothèse revient à introduire une interaction entre les rapports et la durée ou plutôt précisément une fonction de la durée).
Plusieurs méthodes disponibles, on traitera celles basées sur les résidus de Schoenfeld, et l’introduction directe d’une intéraction entre une fonction la durée et les covariables du modèle. Cette dernière fait également office de méthode de correction lorsque la violation de l’hypothèse est jugée trop importante ou problématique du point de vue des résultats obtenus.
Si on regarde les courbes de Kaplan-Meier, leurs croisement non tardif impliquera nécessairement un problème sur cette hypothèse.

8.2.1 Test de Grambsch-Therneau sur les résidus de Schoenfeld

Ce test a été proposé par P.Grambsch et T.Therneau ³ dans un cadre à durée strictement continue. Il repose originellement sur une régression linéaire estimé avec les moindres carrés généralisés (GLS) correction de l’autocorrélation des erreurs avec des sér) . Dans un premier temps pour des raisons plutôt pratiques (informatique), le test a une version moindres carrés ordinaires (OLS). Jusqu’en 2020, tous les logiciels ne proposaient que le test OLS. T.Therneau avec la V3 de package survival a substitué - assez brutalement - le test GLS au test OLS. Si les résultats sont proches dans le cadre d’une durée continue et que le test GLS peut être considéré comme un test exact, cela devient problématique dans une situation de durée discrète/groupée ⁴. Le test OLS reste, à mon sens, la méthode à privilégier dans le cas discret.

Il est également important de souligner que pour P.Grambsch et T.Therneau ⁵ n’est qu’un moyen parmi d’autres d’analyser une violation de l’hypothèse de proportionnalité. Ce n’est pas the solution (comme tout autre test au passage). Le croisement des courbes de séjours peut-être suffisant pour alerter sur cette violation.

Principe du test: consiste à regarder la corrélation entre les résidus de Schoenfeld obtenus directement avec la fonction de score de la vraisemblanc partielle de Cox et une fonction de la durée.

Principe de calcul des résidus

Les résidus bruts sont directement calculés à partir des équations de scores [voir section estimation].
Ils ne sont calculés que pour les observations qui ont connues l’évènement, au moment où un évènement s’est produit.
La somme des résidus pour chaque covariable est égale à 0. Il s’agit de la propriété de l’équation de score à l’équilibre.
On utilise généralement les résidus standardisés (remis à l’échelle / scaled) - par leur variance -. C’est la mesure de cette variance qui distingue le test OLS du test GLS.

Pour une observation dont l’évènement s’est produit en $t_i$, le résidu brut de Schoenfeld pour la covariable $X_k$, après estimation du modèle, est égal à:

\[rs_{ik}=X_{ik}- \sum_{j\in R_i}X_{jk}\frac{e^{X_{j}^{'}b}}{\sum_{j\in R_i}e^{X_{j}^{'}b}}= X_{ik} - E(X_{j\in R_i})\]

Ce résidu est formellement la contribution d’une observation ou d’un moment d’évènement au score. Il se lit comme la différence entre la valeur observée d’une covariable et sa valeur espérée au moment où l’évènement s’est produit.
Si la constance des rapports de risque varie peu les résidus ne doivent pas suivre une tendance précise localement ou globalement, à la hausse ou à la baisse.

Pourquoi?

Par l’exemple, sans censure à droite et en ne considérant que les résidus bruts: Avec un rapport de risque strictement égal à 1 en début d’exposition, une population soumise au risque $R_i=100$ avec 50 hommes et 50 femmes. Si l’hypothèse PH (strictement) respectée, lorsqu’il reste 90 personnes soumises au risque, on devrait avoir 45 hommes et 45 femmes. Avec $R_i=50$, 25 hommes et 25 femmes,…….avec $R_i=10$, 5 hommes et 5 femmes.
Au final l’espérance d’avoir la caractéristique $X$ est toujours égal à 0.5 et les résidus bruts prendront toujours la valeur -.5 si $X=0$ et .5 si $X=1$. En faisant une simple régression linéaire entre les résidus, qui alternent ces deux valeurs, et $t$, le coefficient estimé sera en toute logique très proche de 0.

De manière encore plus simple, cette proportionnalité avec un risque ratio égal à 1 suggère qu’au cours de la durée d’observation, on observe une succession d’un même nombre d’hommes et de femmes qui connaissent l’évènement. Si tous les hommes ou presque avaient observés l’évènement plutôt en début d’éxposition et si toutes les femmes ou presque avaient observé l’évènement plutôt en fin d’exposition, l’hypothèse de proportionnalité pourraient fortement remise en cause.

On trouvera des éléments de calcul du test OLS ici

Avertissement

Test omnibus: Ne pas l’utiliser bien qu’il figure généralement en bas des output. Il n’a pas d’interprétation directe, et les p-value peuvent présenter des valeurs très faibles alors que ce n’est pas le cas pour les covariables prises une à une. Rester comme c’est souvent le cas à un test à un degré de liberté.
Transformations de la durée: n’importe quelle fonction de la durée peut être utilisée pour réaliser le test. On retient généralement les fonctions suivantes: $g(t)=t$ (« identity »), $g(t)=log(t)$, $g(t)=KM(t)$ ou $g(t)=1- S(t$) où $S(t$) est l’estimateur de Kaplan-Meier. Enfin une transformation appelée « rank » est utilisée seulement pour les durées strictement continue ou suffisamment dispersées . Par exemple $t=(0.1,0.5,1,2.6,3)$ donne une transformation $t=(1,2,3,4)$. A savoir : $g(t)=t rend le test relativement sensible aux évènements tardifs lorsque la population restant soumise est peu nombreuse (outliers).
Par défaut Stata, Sas, Python: $g(t)=t$
Par défaut R: $g(t)= 1 - S(t)$

Pour des raisons de reproductibilité dans l’espace des logiciels et dans le temps pour les différentes versions du package survival de R, on ne présente ici que la version OLS.

Test OLS avec $g(t)=t$

Test OLS Grambsch-Therneau avec $g(t)=t$
Variables	chi2	df	P>Chi2
year	0.80	1	0.3720
age	1.61	1	0.2043
surgery	5.54	1	0.0186

Ici l’hypothèse de proportionnalité des risques est questionnable pour la variable surgery. Le risque ratio pourrait ne pas constant dans le temps. Ce n’est pas du tout étonnant, le premier décès pour les personnes opérées d’un pontage n’est observé qu’au bout de 165 jours. Au final, un test était-il bien nécessaire pour arriver à ce constat ???????

Test OLS avec $g(t)=1- S(t)$

Test Grambsch-Therneau avec $g(t)=1- S(t)$
Variables	chi2	df	P>Chi2
year	1.96	1	0.162
age	1.15	1	0.284
surgery	3.96	1	0.046

R-Stata-Sas-Python

Attention seulement version GLS du test depuis le V3 de survival.

Après avoir créer un objet à l’estimation du modèle de Cox, on utilise la fonction cox.zph. Cette fonction utilise par défaut $g(t)=1-S(t)$ où $S(t)$ sont les estimateurs de la courbe de Kaplan-Meier. On peut modifier cette fonction. Il est préférable de conserver cette fonction par défaut.
Test OLS: j’ai récupéré le programme du test antérieur, renommé cox.zphold. On peut le charger simplement, et il est facilement exécutable. Pour le charger: source("https://raw.githubusercontent.com/mthevenin/analyse_duree/main/cox.zphold/cox.zphold.R")

8.2.1.1 Stata

Le test (OLS) est obtenu avec la commande estat phtest, d. Par défaut Stata utilise $g(t)=t$. On peut modifier cette fonction.

8.2.1.2 SAS

Le test (OLS) est disponible depuis quelques années avec l’argument zph sur la ligne proc lifetest. Par défaut SAS utilise $g(t)=t$. On peut modifier cette fonction.

8.2.1.3 Python

Le test (OLS) est donné avec la fonction proportional_hazard_test de la librairie lifelines. La fonction utilise par défaut $g(t)=t$, mais on peut afficher les résultats pour toutes les transformations de $t$ disponibles avec l’option time_transform=’all’.

8.2.2 Intéraction avec la durée

Petit retour sur l’estimation du modèle

Pour estimer le modèle de Cox, les données sont dans un premier temps splitées aux moment où au moins un évènement a été observé.

Sur l’application, avec 2 individus avec la covariable age (rappel: il s’agit de l’âge en $t_0$:

Base spittées sur les intervals d’évènement
id	age	died	$t_0$	$t$
2	51	0	0	1
2	51	0	1	2
2	51	0	2	3
2	51	0	3	5
2	51	1	5	6
3	54	0	0	1
3	54	0	1	2
3	54	0	2	3
3	54	0	3	5
3	54	0	5	6
3	54	0	6	8
3	54	0	8	9
3	54	0	9	12
3	54	1	12	16

Les bornes des intervalles $[t_0;t]$ présentent des valeurs seulement lorsqu’un évènement s’est produit (principe de la vraisemblance partielle). Il n’y a donc pas de valeurs pour $t$ et $t_0$ en $t=4$ pour $id=(2,3)$et $t=7,10,11,13,14,15$ pour $id=3$.
Les deux individus observent l’évènement en $t=6$ pour $id=2$, et en $t=16$ pour $id=3$. Avant ce moment la valeur de la variable évènement/censure (ici $d$) prend toujours la valeur 0, et prend la valeur 1 le jour du décès.

Sur cette base splitée aux moments d’évènement (n=3573), on pourra vérifier facilement que les résultats obtenus par le modèle de Cox sont identiques à ceux obtenus précédemment.

Introduction d’une intéraction avec une fonction de la durée

On a une variable de durée (on prendra $g(t)=t$) qui sera croisée avec la variable surgery.

Le modèle s’écrit:

\[h(t | X,t) = h_0(t)e^{b_1age + b_2year + b_3 surgery + b_4 (surgery\times t)}\]

Le modèle avec cette intéraction donne les résultats suivants:

Modèle de Cox avec une intéraction entre une fonction de la durée et la variable *surgery
Variable	$e^b$	Std.err	z	P>\|z\|	95% IC
year	0.884	0.059	-1.84	0.066	0.776 ; 1.008
age	1.029	0.014	+2.15	0.032	1.003 ; 1.057
$surgery(t_{0+})$	0.173	0.117	-2.60	0.009	0.046 ; 0.649
$surgery\times t$	1.002	0.001	+2.02	0.043	1.000 ; 1.004

On retrouve donc un résultat proche de celui obtenu à partir du test OLS sur les résidus de Schoenfeld pour la variable surgery. Et c’est normal. Avec $g(t)=t$, il a le mérite de pouvoir être interprété directement. Ce qui ne veut pas dire qu’il s’agit de la meilleure solution.

Donc, malgré une hypothèse plutôt forte sur la forme fonctionnelle de l’intéraction, et dans les faits surement pas pertinente, on peut dire que chaque jour le rapport des risques entre personnes opérées et personnes non opérées augmente de +0.2%. Pour plus précis, étant à l’origine <1, l’écart se modère. L’effet de l’opération sur la survie des individus s’estompe donc avec le temps.

A noter

Le modèle n’est plus un modèle à risque proportionnel. La variable surgery n’est plus une variable fixe mais une variable tronquée dynamique qui prend la valeur de $t$ pour les personnes qui ont été opérées d’un pontage avant leur entrée dans le registre de greffe.

Si $surgery = 0$

id	died	$t_0$	$t$
2	0	0	1
2	0	1	2
2	0	2	3
2	0	3	5
2	1	5	6

Si $surgery = 1$ (jusqu’à $t=6$ car aucun décès précoce pour ce groupe)

id	surgery	died	$t_0$	$t$	surgery*t
40	1	0	0	1	1
40	1	0	1	2	2
40	1	0	2	3	3
40	1	0	3	5	5
40	1	1	5	6	6

Exemple pour une variable quantitative (age)

id	age	died	$t_0$	$t$	age*t
2	51	0	0	1	51
2	51	0	1	2	102
2	51	0	2	3	153
2	51	0	3	5	255
2	51	1	5	6	306

L’altération des rapports de risque dépend de la forme fonctionnelle de l’intéraction choisie. Ici la variation dans la durée du rapport des risque est constante, ce qui est une hypothèse assez forte. On a, en quelques sorte, réintroduit une hypothèse de proportionnalité, ici sur le degré d’altération des écarts de risques dans le temps, qui devient lui même strictement constant.

8.2.3 Que faire ?

Ne rien faire

On interprète le risque ratio comme un ratio moyen pendant la durée d’observation (P.Allison). Difficilement soutenable pour l’analyse des effets cliniques, elle peut être envisagée dans d’autres domaines. Attention au nombre de variables qui ne respectent pas l’hypothèse, l’estimation de la baseline du risque pourrait être sensiblement affectée si l’analyse a des visée prédictives. Il convient tout de même lors de l’interprétation, de préciser les variables qui seront analysées sous cette forme très « moyenne » sur la période d’observation.

On peut également adapter cette stratégie du « ne rien faire » selon sens de l’altération des rapports de risque. Si aux cours du temps des écarts de risque, s’accentuent à la hausse comme à la baisse, on peut conserver cet estimateur moyen. Mais si cette non proportionnalité conduit à un changement du sens des rapport de risque je suis moins convaincu de la pertinence de cette stratégie. Encore une fois, et il faut le rappeler, l’estimation des courbes de survie doit permette d’anticiper ce dernier cas de figure.

Il faut également tenir compte de l’intérêt portée par les variables qui présentent un problème par rapport à l’hypothèse. Il n’est peut-être pas nécessaire de complexifier le modèle pour des variables introduites comme simples contrôles.

Mais plus problématique [important]… On sait qu’une des causes du non respect de l’hypothèse peut provenir d’effets de sélection liées à des variables omises ou non observables. En analyse de durée ce problème prend le nom de frailty (fragilité) lorsque cette non homogénéité n’est pas observable. Des estimations, plus complexes, sont possibles dans ce cas, et sont en mesure malgré leur interprétation plutôt difficile de régler le problème. Il convient donc de bien spécifier le modèle au niveau des variables de contrôle observables et disponibles.

Modèle de Cox stratifié

Utiliser la méthode dite de « Cox stratifiée » (non traitée). Utile si l’objectif est de présenter des fonctions de survie prédites ajustées, et si une seule covariable (binaire) présente un problème. Les HR ne seront pas estimés pour la variable qui ne respecte pas l’hypothèse.

Intéraction

Introduire une interaction avec la durée. Cela peut permettre en plus d’enrichir le modèle au niveau de l’interprétation. Valable si peu de covariables présentent des problèmes de stabilité des rapports de risque, dans l’idéal une seule variable. Attention tout de même à la forme de la fonction, dans l’exemple on a contraint l’effet d’interaction à être strictement linéaire, ce qui est une hypothèse plutôt forte…. on introduit de nouveau une contrainte de proportionnalité dans le modèle.

Modèles alternatifs

Utiliser un modèle alternatif: modèles paramétriques à risques proportionnels si la distribution du risque s’ajuste bien, le modèle paramétrique « flexible » de Parmar-Royston ou un modèle à temps discret. Pour la dernière solution, on peut également corriger la non proportionnalité avec l’introduction d’une intéraction. Si on ne le fait pas, les risques prédits, par définition des probabilités conditionnelles, resteront toujours dans les bornes contrairement au modèle de Cox.
Utiliser un modèle non paramétrique additif dit d’Aalen ou une de ses variantes (non traité). Mais ces modèles, dont les résultats sont présentés par des graphiques, se commentent assez difficilement.

Forêt aléatoire
Autre méthode : les forêts aléatoires. L.Breiman a dès le départ proposé une estimation des modèles de survie par cette méthode. Par définition, pas sensible à l’hypothèse PH. Mais cela reste des méthodes à finalité prédictive, moins riche en interprétation.

Se reporter à la définition des grandeurs dans la section Théorie↩︎
$n! = (n)\times(n-1)\times(n-2)\times....\times3\times2\times1$↩︎
Il s’agit bien de la personne qui maintient le package survival dans R↩︎
Pour les personnes utilisant R, je donne un moyen pour récupérer et exécuter le test OLS sous R↩︎
Se reporter à leur ouvrage Modeling Survival Data: Extending the Cox Model (2001)↩︎

--- #title: "**Le modèle de Cox**" filters: - lightbox lightbox: auto --- # **Le modèle de Cox** On peut ignorer la partie sur l’estimation du modèle. On retiendra tout de même qu'il est déconseillé d'utiliser la méthode dite *exacte* pour la correction de la vraisemblance, qui ne peut matériellement fonctionner qu’avec un nombre très limité d’évènements observés simultanément. Ce qui est plutôt rare avec des données à durées discrètes ou groupées, très fréquentes dans les sciences sociales. ## Le modèle semi-paramétrique de Cox ### La vraisemblance partielle et estimation des paramètres On se situe dans une situation où la durée est mesurée sur une échelle strictement continue. Il ne peut donc y avoir qu'un seul évènement observé en $t_i$ (idem pour la censure). On peut représenter le processus aléatoire d'une analyse de survie en présence de censure à droite, avec l'équation de vraisemblance suivante: $$L_i=f(t_i)^{\delta_i}S(t_i)^{1-\delta_i}$$ * $f(t_i)$ est la valeur de la fonction de densité en $t_i$ * $S(t_i)$ est la valeur de la fonction de survie en $t_i$ * $\delta_i=1$ si l'évènement est observé: $L_i=f(t_i)$ * $\delta_i=0$ si l'observation est censurée: $L_i=S(t_i)$ **Vraisemblance partielle de Cox** Comme $f(t_i)=h(t_i)\times S(t_i)$ ^[Se reporter à la définition des grandeurs dans la section *Théorie*], on obtient: $L_i=[h(t_i)S(t_i)]^{\delta_i}S(t_i)^{1-\delta_i} = h(t_i)^{\delta_i}S(t_i)$. Pour $i=1,2,.....,n$, la vraisemblance s'ecrit donc: $L_i=\prod_{i=1}^{n}h(t_i)^{\delta_i}S(t_i)$. On peut réécrire cette vraisemblance en la multipliant et en la divisant par: $\sum_{j\in R_i}h(t_i)$, où $j\in R_i$ est l'ensemble des observations soumises au risque en $t_i$. $$L=\prod_{i=1}^{n}\left[h(t_i)\frac{\sum_{j\in R}h(t_i)}{\sum_{j\in R}h(t_i)}\right]^{\delta_i}S(t_i)= \prod_{i=1}^{n}\left[\frac{h(t_i)}{\sum_{j\in R_i}h(t_i)}\right]^{\delta_i}\sum_{j\in R_i}h(t_i)^{\delta_i}S(t_i)$$ La vraisemblance partielle retient seulement le premier terme de la vraisemblance, soit: $$PL=\prod_{i=1}^{n}\left[\frac{h(t_i)}{\sum_{j\in R}h(t_i)}\right]^{\delta_i}$$ Une fois remplacée la valeur de $h(t_i)$ par son expression en tant que modèle à risques proportionnels, la vraisemblance partielle ne dépendra plus de la durée. **Mais elle va dépendre de l'ordre d'arrivée des évènements, c'est à dire leur rang**. *Remarque*: pour les observations censurées($\delta_i=0$), $PL=1$. Toutefois, ces censures à droite entrent dans l'expression $\sum_{j\in R}h(t_i)$ tant qu'elles sont soumises au risque. En remplaçant donc $h(t_i)$ par l'expression $h_0(t)e^{X_i^{'}b}$: $$PL=\prod_{i=1}^{n}\left[\frac{h_0(t)e^{X_{i}^{'}b}}{\sum_{j\in R_i}h_0(t)e^{X_{j}^{'}b}}\right]^{\delta_i} =\prod_{i=1}^{n}\left[\frac{e^{X_i^{'}b}}{\sum_{j\in R_i}e^{X_{j}^{'}b}}\right]^{\delta_i}$$ L'expression $\frac{e^{Xb}}{\sum_{j\in R}e^{Xb}}$ est donc bien une probabilité, et la vraisemblance partielle est donc bien un produit de probabilités. Pour un individu ayant connu l'évènement, la contribution à la vraisemblance partielle est **la probabilité que l'individu observe l'évènement en $t_i$ sachant qu'un évènement (et un seul) s'est produit**. - Si $\delta_i = 0$: $PL_i = 1$ - Si $\delta_i = 1$: $PL_i =\frac{e^{X_i^{'}b}}{\sum_{j\in R_i}e^{X_{j}^{'}b}}$ *Condition nécessaire: pas d'évènement simultané*: en présence d'évènements mesurés simultanément, l’estimation de la vraisemblance doit faire l'objet d'une correction. *Correction de la vraisemblance avec des évènements simultanés*: * La ***méthode dite exacte***: Comme il ne doit pas y avoir d’évènement simultané, on va introduire à la vraisemblance partielle toutes les permutations possibles des évènements observés au même moment. Bien qu'en $t_i$ on observe au même moment l’évènement pour 2 observations (A,B) une métrique temporelle plus précise permettrait de savoir si A s'est produit avant B ou B s'est produit avant A (2 permutations). Comme le nombre de permutations est calculé à l'aide d'une factorielle ^[$n! = (n)\times(n-1)\times(n-2)\times....\times3\times2\times1$], avec 3 évènements mesurés simultanément, on obtient 6 permutations ($3\times2\times1$). Problème: le nombre de permutations pour chaque $t_i$ peut devenir très vite particulièrement élevé. Par exemple pour 10 évènements simultanés, le nombre de permutations est égal à 3,628,800. Le temps de calcul devient extrêmement long, et ce type de correction totalement inopérant. * La ***méthode dite de Breslow***: il s'agit d'une approximation de la méthode exacte permettant de ne pas avoir à intégrer chaque permutation. *Cette approximation est utilisée par défaut par les logiciels Sas et Stata*. * La ***méthode dite d'Efron***: elle corrige l'approximation de Breslow, et est jugée plus proche de la méthode exacte. *C'est la méthode utilisée par défaut avec le logiciel R*, et elle est disponible avec les autres applications. ### Estimation des paramètres On utilise la méthode habituelle, à savoir la maximisation de la log-vraisemblance (ici partielle). * Conditions de premier ordre: calcul des équations de score à partir des dérivées partielles. Solution: $\frac{\partial log(PL)}{\partial{b_k}}=0$. On ne peut pas obtenir de solution numérique directe. Remarque: les équations de score sont utilisées pour tester la validité de l'hypothèse de constance des rapports de risque pour calculer les **résidus de Schoenfeld** (voir plus loin). * Conditions de second ordre: calcul des dérivées secondes qui permettent d'obtenir la matrice d'information de Fisher et la matrice des variances-covariances des paramètres. * Comme il n'y a pas de solution numérique directe, on utilise un algorithme d'optimisation (ex: Newton-Raphson) à partir des équations de score et de la matrice d'information de Fisher. **Eléments de calcul** En logarithme (sans évènement simultané), la vraisemblance partielle s'ecrit: $$pl(b)=\sum_{i=1}^n\delta_i\left(log(e^{X_{i}^{'}b})-log\sum_{j\in R_i}e^{X_{j}^{'}b}\right)$$ $$pl(b)=\sum_{i=1}^n\delta_i\left(X_{i}^{'}b-log\sum_{j\in R_i}e^{X_{j}^{'}b}\right)$$ Calcul de l'équation de score pour une covariable $X_k$: $$\frac{\partial pl(b)}{\partial{b_k}}=\sum_{i=1}^n\delta_i\left(X_{ik}-\sum_{j\in R_i}X_{jk}\frac{e^{X_{j}^{'}b}}{\sum_{j\in R_i}e^{X_{j}^{'}b}}\right)$$ Comme $\frac{e^{X_{j}b}}{\sum_{j\in R}e^{X_{j}b}}$ est une probabilité, et $\sum_{j\in R}X_{ik}\times p_i$ est l'espérance (la moyenne) $E(X_k)$ d'avoir la caractéristique $X_k$ lorsqu'un évènement a été observé. Au final: $$\frac{\partial lp(b)}{\partial{b_k}}= \sum_{i=1}^n\delta_i\left(X_{ik} - E(X_{j\in R_i,k}) \right)$$ Cette expression va permettre d'analyser le respect ou non de l’hypothèse de risques proportionnels via les *résidus de Schoenfeld*. ### Lecture des résultats Comme il s'agit d'un modèle à risque proportionnel, **les rapports de risques sont constants pendant toute la période d'observation**. Il s'agit d'une **propriété de l'estimation**. ***Covariable binaire (indicatrice)*** $X=(0,1)$: $RR=\frac{h(t\ |\ X=1)}{h(t\ |\ X=0)}=e^b$. A chaque moment de la durée $t$, le risque d'observer l'évènement est $e^b$ fois plus important/plus faible pour $X=1$ que pour $X=0$. ***Covariable quantitative*** (fixe dans le temps) $RR=\frac{h(t\ |\ X=a+c)}{h(t\ |\ X=a)}=e^{c \times{b}}$. On prendra pour illustrer une variable type âge au début de l’exposition au risque (a) et un delta de comparaison avec un âge inférieur c. Si $c=1$ (résultat de l'estimation): A un âge donnée, le risque de connaitre l'évènement est $e^b$ fois inférieur/supérieur à celui d'une personne qui a un an de moins. **Exemple pour les insuffisances cardiaques** * Correction de la vraisemblance: méthode d'Efron * Nombre d'observations: 103 * Nombre de décès: 75 * Log-Vraisemblance: -289.30639 | Variables | logRR | Std.Err | z | $P>|z|$| 95% IC | |-----------|---------|---------|-------|--------|-----------------| | year | -0.119 | 0.0673 | -1.78 | 0.076 | -0.2516;+0.0124 | | age | +0.0296 | 0.0135 | 2.19 | 0.029 | +0.0031;+0.0561 | | surgery | -0.9873 | 0.4363 | -2.26 | 0.024 | -1.8424;-0.1323 | : Cox: log Hazard Ratio (Risks Ratio) {.striped .hover} | Variables | RR | Std.Err | z | $P>|z|$| 95%CI | |-----------|---------|---------|-------|--------|------------------| | year | 0.8872 | 0.0597 |-1.78 | 0.076 | 0.7775; 1.0124 | | age | 1.0300 | 0.0139 | 2.19 | 0.029 | 1.0031; 1.0577 | | surgery | 0.3726 | 0.1625 |-2.26 | 0.024 | 0.1584; 0.8761 | : Cox: Hazard Ratio (Risks Ratio) {.striped .hover} On retrouve les des tests non paramétriques pour l'opération, à savoir qu’un pontage réduit les risques journaliers de décès pendant la période d’observation (augmente la durée de survie). De la même manière, plus on entre à un âge élevé dans la liste d'attente plus le risque de décès augmente. La variable *year*, qui traduit des progrès en médecine, renvoie à une réduction plutôt modérée du risque journalier de décès durant l'attente d'une greffe. **R-Stata-Sas-Python** ::: panel-tabset ### R Le modèle est estimé avec la fonction **`coxph`** de la librairie **`survival`**. Hors options, la syntaxe est identiques aux fonctions `survfit` et `survdif`. #### Stata Le modèle est estimé avec la commande **`stcox`**. #### SAS Le modèle est estimé avec la **`proc phreg`**. #### Python Avec la librairie `lifelines`, le modèle est estimé avec la fonction **`CoxPHFitter`**. Avec la librairie `statmodels`, il est estimé avec la fonction `smf.phreg`. ::: ## Analyse de la constance des rapports de risque * Les rapports de risque (RR) estimés par le modèle sont contraints à être constant sur toute la période d’observation. C’est une hypothèse forte. * Le respect de cette hypothèse doit être analysé, en particulier pour le modèle de Cox où la baseline du risque est habituellement estimée à l’aide de ces rapports (par exemple la méthode dite de Breslow, non traitée). En post-estimation, les valeurs estimées du risque pourront présenter des valeurs aberrantes si on dévie trop de constance, en particulier en obtenant des négatives des taux de risque. * Analyser cette hypothèse revient à introduire une interaction entre les rapports et la durée ou plutôt précisément une fonction de la durée). * Plusieurs méthodes disponibles, on traitera celles basées sur les **résidus de Schoenfeld**, et l'introduction directe d'une intéraction entre une fonction la durée et les covariables du modèle. Cette dernière fait également office de méthode de correction lorsque la violation de l'hypothèse est jugée trop importante ou problématique du point de vue des résultats obtenus. * Si on regarde les courbes de Kaplan-Meier, leurs croisement non tardif impliquera nécessairement un problème sur cette hypothèse. ### Test de Grambsch-Therneau sur les résidus de Schoenfeld Ce test a été proposé par P.Grambsch et T.Therneau ^[Il s'agit bien de la personne qui maintient le package **`survival`** dans R] dans un cadre à durée strictement continue. Il repose originellement sur une régression linéaire estimé avec les moindres carrés généralisés (GLS) correction de l'autocorrélation des erreurs avec des sér) . Dans un premier temps pour des raisons plutôt pratiques (informatique), le test a une version moindres carrés ordinaires (OLS). Jusqu'en 2020, tous les logiciels ne proposaient que le test OLS. T.Therneau avec la V3 de package `survival` a substitué - assez brutalement - le test GLS au test OLS. Si les résultats sont proches dans le cadre d'une durée continue et que le test GLS peut être considéré comme un test *exact*, cela devient problématique dans une situation de durée discrète/groupée ^[Pour les personnes utilisant R, je donne un moyen pour récupérer et exécuter le test OLS sous R]. Le test OLS reste, à mon sens, la méthode à privilégier dans le cas discret. Il est également important de souligner que pour P.Grambsch et T.Therneau ^[Se reporter à leur ouvrage *Modeling Survival Data: Extending the Cox Model* (2001)] n'est qu'un moyen parmi d'autres d'analyser une violation de l'hypothèse de proportionnalité. Ce n'est pas *the solution* (comme tout autre test au passage). Le croisement des courbes de séjours peut-être suffisant pour alerter sur cette violation. ***Principe du test***: consiste à regarder la corrélation entre les **résidus de Schoenfeld** obtenus directement avec la fonction de score de la vraisemblanc partielle de Cox et une fonction de la durée. ***Principe de calcul des résidus*** * Les résidus *bruts* sont directement calculés à partir des équations de scores [voir section estimation]. * Ils ne sont calculés que pour les observations qui ont connues l’évènement, au moment où un évènement s’est produit. * La somme des résidus pour chaque covariable est égale à 0. Il s'agit de la propriété de l’équation de score à l’équilibre. * On utilise généralement les *résidus standardisés* (*remis à l'échelle* / scaled) - par leur variance -. C'est la mesure de cette variance qui distingue le test OLS du test GLS. Pour une observation dont l’évènement s’est produit en $t_i$, le *résidu brut de Schoenfeld* pour la covariable $X_k$, après estimation du modèle, est égal à: $$rs_{ik}=X_{ik}- \sum_{j\in R_i}X_{jk}\frac{e^{X_{j}^{'}b}}{\sum_{j\in R_i}e^{X_{j}^{'}b}}= X_{ik} - E(X_{j\in R_i})$$ * Ce résidu est formellement la contribution d'une observation ou d'un moment d'évènement au score. Il se lit comme la différence entre la valeur observée d'une covariable et sa valeur espérée au moment où l'évènement s'est produit. * Si la constance des rapports de risque varie peu les résidus ne doivent pas suivre une tendance précise localement ou globalement, à la hausse ou à la baisse. Pourquoi? Par l'exemple, sans censure à droite et en ne considérant que les résidus bruts: Avec un rapport de risque strictement égal à 1 en début d'exposition, une population soumise au risque $R_i=100$ avec 50 hommes et 50 femmes. Si l'hypothèse PH (strictement) respectée, lorsqu'il reste 90 personnes soumises au risque, on devrait avoir 45 hommes et 45 femmes. Avec $R_i=50$, 25 hommes et 25 femmes,.......avec $R_i=10$, 5 hommes et 5 femmes. Au final l'espérance d'avoir la caractéristique $X$ est toujours égal à 0.5 et les résidus bruts prendront toujours la valeur -.5 si $X=0$ et .5 si $X=1$. En faisant une simple régression linéaire entre les résidus, qui alternent ces deux valeurs, et $t$, le coefficient estimé sera en toute logique très proche de 0. De manière encore plus simple, cette proportionnalité avec un risque ratio égal à 1 suggère qu'au cours de la durée d'observation, on observe une succession d'un même nombre d'hommes et de femmes qui connaissent l'évènement. Si tous les hommes ou presque avaient observés l'évènement plutôt en début d'éxposition et si toutes les femmes ou presque avaient observé l'évènement plutôt en fin d'exposition, l'hypothèse de proportionnalité pourraient fortement remise en cause. On trouvera des éléments de calcul du test OLS [ici](https://mthevenin.github.io/analyse_duree/annexes/residusch.html) ::: callout-warning - **Test omnibus**: Ne pas l'utiliser bien qu'il figure généralement en bas des output. Il n'a pas d'interprétation directe, et les p-value peuvent présenter des valeurs très faibles alors que ce n'est pas le cas pour les covariables prises une à une. Rester comme c'est souvent le cas à un test à un degré de liberté. - **Transformations de la durée**: n’importe quelle fonction de la durée peut être utilisée pour réaliser le test. On retient généralement les fonctions suivantes: $g(t)=t$ (« identity »), $g(t)=log(t)$, $g(t)=KM(t)$ ou $g(t)=1- S(t$) où $S(t$) est l’estimateur de Kaplan-Meier. Enfin une transformation appelée « rank » est utilisée seulement pour les durées strictement continue ou suffisamment dispersées . Par exemple $t=(0.1,0.5,1,2.6,3)$ donne une transformation $t=(1,2,3,4)$. A savoir : $g(t)=t rend le test relativement sensible aux évènements tardifs lorsque la population restant soumise est peu nombreuse (outliers). * Par défaut Stata, Sas, Python: $g(t)=t$ * Par défaut R: $g(t)= 1 - S(t)$ ::: Pour des raisons de reproductibilité dans l'espace des logiciels et dans le temps pour les différentes versions du package `survival` de R, on ne présente ici que la version OLS. **Test OLS avec $g(t)=t$** | Variables | chi2 | df | P>Chi2 | |-----------|---------|-------|--------| | year | 0.80 | 1 | 0.3720 | | age | 1.61 | 1 | 0.2043 | | surgery | 5.54 | 1 | 0.0186 | : Test OLS Grambsch-Therneau avec $g(t)=t$ Ici l'hypothèse de proportionnalité des risques est questionnable pour la variable *surgery*. Le risque ratio pourrait ne pas constant dans le temps. Ce n'est pas du tout étonnant, le premier décès pour les personnes opérées d'un pontage n'est observé qu'au bout de 165 jours. Au final, un test était-il bien nécessaire pour arriver à ce constat ??????? **Test OLS avec $g(t)=1- S(t)$** | Variables | chi2 | df | P>Chi2 | |-----------|---------|-------|--------| | year | 1.96 | 1 | 0.162 | | age | 1.15 | 1 | 0.284 | | surgery | 3.96 | 1 | 0.046 | : Test Grambsch-Therneau avec $g(t)=1- S(t)$ **R-Stata-Sas-Python** ::: panel-tabset ### R ***Attention seulement version GLS du test depuis le V3 de survival***. * Après avoir créer un objet à l’estimation du modèle de Cox, on utilise la fonction **`cox.zph`**. Cette fonction utilise par défaut $g(t)=1-S(t)$ où $S(t)$ sont les estimateurs de la courbe de Kaplan-Meier. On peut modifier cette fonction. Il est préférable de conserver cette fonction par défaut. * Test OLS: j’ai récupéré le programme du test antérieur, renommé **`cox.zphold`**. On peut le charger simplement, et il est facilement exécutable. Pour le charger: `source("https://raw.githubusercontent.com/mthevenin/analyse_duree/main/cox.zphold/cox.zphold.R")` #### Stata Le test (OLS) est obtenu avec la commande **`estat phtest, d`**. Par défaut Stata utilise $g(t)=t$. On peut modifier cette fonction. #### SAS Le test (OLS) est disponible depuis quelques années avec l’argument **`zph`** sur la ligne `proc lifetest`. Par défaut SAS utilise $g(t)=t$. On peut modifier cette fonction. #### Python Le test (OLS) est donné avec la fonction **`proportional_hazard_test`** de la librairie `lifelines`. La fonction utilise par défaut $g(t)=t$, mais on peut afficher les résultats pour toutes les transformations de $t$ disponibles avec l’option `time_transform=’all’`. ::: ### Intéraction avec la durée ***Petit retour sur l'estimation du modèle*** Pour estimer le modèle de Cox, les données sont dans un premier temps splitées aux moment où au moins un évènement a été observé. Sur l'application, avec 2 individus avec la covariable *age* (rappel: il s'agit de l'âge en $t_0$: | id | age | died |$t_0$| $t$ | |-----|----------|------|-----|-----| | 2 | 51 | 0 | 0 | 1 | | 2 | 51 | 0 | 1 | 2 | | 2 | 51 | 0 | 2 | 3 | | 2 | 51 | 0 | 3 | 5 | | 2 | 51 | 1 | 5 | 6 | | 3 | 54 | 0 | 0 | 1 | | 3 | 54 | 0 | 1 | 2 | | 3 | 54 | 0 | 2 | 3 | | 3 | 54 | 0 | 3 | 5 | | 3 | 54 | 0 | 5 | 6 | | 3 | 54 | 0 | 6 | 8 | | 3 | 54 | 0 | 8 | 9 | | 3 | 54 | 0 | 9 | 12 | | 3 | 54 | 1 | 12 | 16 | : Base spittées sur les intervals d'évènement Les bornes des intervalles $[t_0;t]$ présentent des valeurs seulement lorsqu'un évènement s'est produit (principe de la vraisemblance partielle). Il n'y a donc pas de valeurs pour $t$ et $t_0$ en $t=4$ pour $id=(2,3)$et $t=7,10,11,13,14,15$ pour $id=3$. Les deux individus observent l'évènement en $t=6$ pour $id=2$, et en $t=16$ pour $id=3$. Avant ce moment la valeur de la variable évènement/censure (ici $d$) prend toujours la valeur 0, et prend la valeur 1 le jour du décès. Sur cette base *splitée* aux moments d'évènement (n=3573), on pourra vérifier facilement que les résultats obtenus par le modèle de Cox sont identiques à ceux obtenus précédemment. *Introduction d'une intéraction avec une fonction de la durée* On a une variable de durée (on prendra $g(t)=t$) qui sera croisée avec la variable *surgery*. Le modèle s'écrit: $$h(t | X,t) = h_0(t)e^{b_1age + b_2year + b_3 surgery + b_4 (surgery\times t)}$$ Le modèle avec cette intéraction donne les résultats suivants: | Variable | $e^b$ | Std.err | z | P\>\|z\| | 95% IC | |-------------------|-------|---------|-------|----------|---------------| | year | 0.884 | 0.059 | -1.84 | 0.066 | 0.776 ; 1.008 | | age | 1.029 | 0.014 | +2.15 | 0.032 | 1.003 ; 1.057 | | $surgery(t_{0+})$ | 0.173 | 0.117 | -2.60 | 0.009 | 0.046 ; 0.649 | | $surgery\times t$ | 1.002 | 0.001 | +2.02 | 0.043 | 1.000 ; 1.004 | : Modèle de Cox avec une intéraction entre une fonction de la durée et la variable *surgery On retrouve donc un résultat proche de celui obtenu à partir du test OLS sur les résidus de Schoenfeld pour la variable *surgery*. Et c'est normal. Avec $g(t)=t$, il a le mérite de pouvoir être interprété directement. Ce qui ne veut pas dire qu'il s'agit de la meilleure solution. Donc, malgré une hypothèse plutôt forte sur la forme fonctionnelle de l'intéraction, et dans les faits surement pas pertinente, on peut dire que chaque jour le rapport des risques entre personnes opérées et personnes non opérées augmente de +0.2%. Pour plus précis, étant à l'origine <1, l'écart se modère. L'effet de l'opération sur la survie des individus s'estompe donc avec le temps. **A noter** - Le modèle n'est plus un modèle à risque proportionnel. La variable *surgery* n'est plus une variable **fixe** mais une variable tronquée dynamique qui prend la valeur de $t$ pour les personnes qui ont été opérées d'un pontage avant leur entrée dans le registre de greffe. Si $surgery = 0$ | id | surgery | died |$t_0$| $t$ | surgery*t | |-----|----------|------|-----|-----|------------| | 2 | 0 | 0 | 0 | 1 | 0 | | 2 | 0 | 0 | 1 | 2 | 0 | | 2 | 0 | 0 | 2 | 3 | 0 | | 2 | 0 | 0 | 3 | 5 | 0 | | 2 | 0 | 1 | 5 | 6 | 0 | Si $surgery = 1$ (jusqu'à $t=6$ car aucun décès précoce pour ce groupe) | id | surgery | died |$t_0$| $t$ | surgery*t | |-----|----------|------|-----|-----|-----------| | 40 | 1 | 0 | 0 | 1 | 1 | | 40 | 1 | 0 | 1 | 2 | 2 | | 40 | 1 | 0 | 2 | 3 | 3 | | 40 | 1 | 0 | 3 | 5 | 5 | | 40 | 1 | 1 | 5 | 6 | 6 | Exemple pour une variable quantitative (*age*) | id | age | died |$t_0$| $t$ | age*t | |-----|----------|------|-----|-----|-----------| | 2 | 51 | 0 | 0 | 1 | 51 | | 2 | 51 | 0 | 1 | 2 | 102 | | 2 | 51 | 0 | 2 | 3 | 153 | | 2 | 51 | 0 | 3 | 5 | 255 | | 2 | 51 | 1 | 5 | 6 | 306 | - L'altération des rapports de risque dépend de la forme fonctionnelle de l'intéraction choisie. Ici la variation dans la durée du rapport des risque est constante, ce qui est une hypothèse assez forte. On a, en quelques sorte, réintroduit une hypothèse de proportionnalité, ici sur le degré d'altération des écarts de risques dans le temps, qui devient lui même strictement constant. ### Que faire ? ***Ne rien faire*** On interprète le risque ratio comme un ratio moyen pendant la durée d’observation (P.Allison). Difficilement soutenable pour l’analyse des effets cliniques, elle peut être envisagée dans d’autres domaines. Attention au nombre de variables qui ne respectent pas l’hypothèse, l’estimation de la baseline du risque pourrait être sensiblement affectée si l'analyse a des visée prédictives. Il convient tout de même lors de l’interprétation, de préciser les variables qui seront analysées sous cette forme très « moyenne » sur la période d’observation. On peut également adapter cette stratégie du « ne rien faire » selon sens de l’altération des rapports de risque. Si aux cours du temps des écarts de risque, s’accentuent à la hausse comme à la baisse, on peut conserver cet estimateur moyen. Mais si cette non proportionnalité conduit à un changement du sens des rapport de risque je suis moins convaincu de la pertinence de cette stratégie. Encore une fois, et il faut le rappeler, l'estimation des courbes de survie doit permette d'anticiper ce dernier cas de figure. Il faut également tenir compte de l’intérêt portée par les variables qui présentent un problème par rapport à l’hypothèse. Il n’est peut-être pas nécessaire de complexifier le modèle pour des variables introduites comme simples contrôles. Mais plus problématique [**important**]… On sait qu’une des causes du non respect de l’hypothèse peut provenir d’effets de sélection liées à des variables omises ou non observables. En analyse de durée ce problème prend le nom de ***frailty*** (fragilité) lorsque cette non homogénéité n’est pas observable. Des estimations, plus complexes, sont possibles dans ce cas, et sont en mesure malgré leur interprétation plutôt difficile de régler le problème. Il convient donc de bien spécifier le modèle au niveau des variables de contrôle observables et disponibles. ***Modèle de Cox stratifié*** Utiliser la méthode dite de « Cox stratifiée » (non traitée). Utile si l’objectif est de présenter des fonctions de survie prédites ajustées, et si une seule covariable (binaire) présente un problème. Les HR ne seront pas estimés pour la variable qui ne respecte pas l'hypothèse. ***Intéraction*** Introduire une interaction avec la durée. Cela peut permettre en plus d’enrichir le modèle au niveau de l’interprétation. Valable si peu de covariables présentent des problèmes de stabilité des rapports de risque, dans l'idéal une seule variable. Attention tout de même à la forme de la fonction, dans l’exemple on a contraint l’effet d’interaction à être strictement linéaire, ce qui est une hypothèse plutôt forte.... on introduit de nouveau une contrainte de proportionnalité dans le modèle. ***Modèles alternatifs*** Utiliser un modèle alternatif: modèles paramétriques à risques proportionnels si la distribution du risque s’ajuste bien, le modèle paramétrique « flexible » de Parmar-Royston ou un modèle à temps discret. Pour la dernière solution, on peut également corriger la non proportionnalité avec l'introduction d'une intéraction. Si on ne le fait pas, les risques prédits, par définition des probabilités conditionnelles, resteront toujours dans les bornes contrairement au modèle de Cox. Utiliser un modèle non paramétrique additif dit d’Aalen ou une de ses variantes (non traité). Mais ces modèles, dont les résultats sont présentés par des graphiques, se commentent assez difficilement. **Forêt aléatoire** Autre méthode : les forêts aléatoires. L.Breiman a dès le départ proposé une estimation des modèles de survie par cette méthode. Par définition, pas sensible à l’hypothèse PH. Mais cela reste des méthodes à finalité prédictive, moins riche en interprétation.

Base spittées sur les intervals d’évènement
id	age	died	\(t_0\)	\(t\)
2	51	0	0	1
2	51	0	1	2
2	51	0	2	3
2	51	0	3	5
2	51	1	5	6
3	54	0	0	1
3	54	0	1	2
3	54	0	2	3
3	54	0	3	5
3	54	0	5	6
3	54	0	6	8
3	54	0	8	9
3	54	0	9	12
3	54	1	12	16