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

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$.

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 seront utilisées pour tester la validité de l’hypothèse de constance des rapports de risque pour calculer les résidus de Schoenfeld… Le grand sujet (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 strictement 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$ (jour, mois, années…), le risque d’observer l’évènement est $e^b$ fois plus important/plus faible pour $X=1$ que pour $X=0$.

Covariable mesurée (quantitativement) (mais 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 la même lecture que lors 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 (ie 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.
Les résultats apparaissent, sans connaissance médicale très poussée, conforme à ce qu’on pourrait attendre.

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 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.1.3.3 SAS

Le modèle est estimé avec la proc 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 mais à faire un jour). En post-estimation, les valeurs estimées du risque pourront présenter des valeurs aberrantes , si on dévie trop de cette constance, en particulier en obtenant des négatives des taux de risque.
Important: analyser cette hypothèse revient à introduire une interaction entre les rapports de risque 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, leur croisement non tardif impliquera nécessairement un problème sur cette hypothèse. Il est donc important de les analyser en amont d’un modèle.

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 (Sas, Stata, Python et Survival v2). T.Therneau avec la Version 3 de package survival a substitué - assez brutalement - un test GLS⁵ au test OLS. Et cela pose de gros problème, en particulier en sciences sociales. Ceci a été montré de manière très convaincante par Shawna K. Metzger.

On doit également 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. Ou tout simplement bien regarder les données qu’on a en face des yeux…. par exemple pour la variable surgery, il y a t-il besoin de faire un test? Regardez bien.

Principe du test: consiste à regarder la corrélation entre les résidus de Schoenfeld obtenus directement avec la fonction de score de la vraisemblance 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.
OLS: variance des estimateurs du modèles
GLS: calcule à chaque des variances à chaque évènement. Possible instabilité lorsque le nombre de personnes restent soumises au risque est faible. La moyenne de ces variances est égale à la variance des estimateurs du modèle.

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 dans le temps, 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_0=100$) avec 50 hommes et 50 femmes. Si l’hypothèse PH est 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.

En revanche, 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. Et on obtiendrait au final un rapport de risque moyen égal à 1.

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-values peuvent présenter des valeurs très faibles alors que ce n’est pas le cas pour les covariables prises une à une. Rester si possible le cas à un test à un degré de liberté, variable par variable.
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,5)$. 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)$

TRES IMPORTANT: Pourquoi la version exacte (GLS) du test pose problème

Il sagit de la version d’origine du test (1994) mais implémenté seulement en 2020 dans la version 3 du package survival de R. Les autres logiciels ne l’implémente pas.
S.Metzger a effectué une décomposition complète de sa statistique dans sa version exacte lien pour les plus acharné.e.s et l’a comparé à la version appelé “approximation” (je préfère “version Ols”). Résultats:
- Un terme résiduel se situant au numérateur de la statistique, qui traduit le degré de corrélation entre les covariables est présent. Ce qui signifie:
- Qu’en présence de corrélations entre covariables, même limitées, la statistique du test augmente artificiellement, et produira tout aussi artificiellement mais mécaniquement des pvalues très faibles. Le test exact s’écarte alors fortement du test OLS, alors qu’il ne devrait être juste qu’un peu plus précis. Si les travaux de S.Metzger s’appuient bien sur des simulations, ses résultats sont vérifiés avec des données réelles qu’elle a également confrontée. C’est également ici le cas avec les données du support comme celles des TP. Pour la base utilisée dans ce support, c’est la corrélation entre les variables surgery et year qui explique les forte différences pour la deuxième variable entre les deux versions du test.
Cela signifie également que le test, qui consiste à synthétiser l’approche par l’intéraction (voir plus bas), donnera des résultats contradictoires avec cette dernière.
Dans R, et seulement dans R:
- T.Therneau, proprétaire du package survival et visiblement discutant du travail de S.Metzger, n’a toujours pas remis le test OLS dans le package, ce qui est quand même problématique au niveau reproductibilité et réplicabilité dans le temps et dans l’espace des logiciels. C’est facheux.
- Ne pas paniquer, j’ai récupéré le test OLS et on peut le charger et l’utiliser très facilement. Voir la section programmation dédiée à R pour le détail de la démarche (ou directement ce lien).

On ne présentera donc ici que la version OLS (idem Stata, Python, Sas .

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-Python-Sas

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")

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

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’.

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.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), c’est le principe même du modèle. Difficilement soutenable pour l’analyse des effets cliniques, il faut reconnaître que ceci peut être envisagé dans d’autres domaines, comme en sciences sociales si on ne veut pas aller plus loin. 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. La lecture des résultats ne sera pas affecté. 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.

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. Si on souhaite corriger la non proportionnalité avec une intéraction, il est inutile de complexifier le modèle pour des variables introduites comme simples contrôles… qui ne sont là que pour jouer le rôle de contrôle. Et justement la vérification de la proportionnalité peut améner à s’intérroger sur la bonne spécification du modèle.

En effet, 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: Le risque de base pour chaque groupe est alors estimé en amont.

Intéraction

Introduire une interaction avec la durée.
Cela permet d’enrichir le modèle au niveau de l’interprétation. Valable si peu de covariables présentent des problèmes d’uniformité dans le temps 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, cela facilité la lecture, mais c’est une hypothèse très forte…. car au final 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 à durée discrète/groupée (voir section suivante). 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 qui sont par définition des probabilités conditionnelles, resteront toujours dans les bonnes 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 seront présentés par des graphiques, se commentent assez difficilement. Et dont l’utilisation est devenue je crois plus que rare.

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é très prédictive, moins riche en interprétation.

Se reporter à la définition des grandeurs dans la section Théorie↩︎
Dans un entretien en 1994, Cox a précisé qu’il avait cette manipulation pour voir, sans idée préconcue↩︎
$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↩︎
Moindres carrés généralisés↩︎
Se reporter à leur ouvrage Modeling Survival Data: Extending the Cox Model (2001)↩︎

--- #title: "**Cox**" --- # 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 ### 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$. 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^[Dans un entretien en 1994, Cox a précisé qu'il avait cette manipulation pour voir, sans idée préconcue], soit: ::: {.box_img} $$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}$: ::: {.box_img} $$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 seront utilisées pour tester la validité de l'hypothèse de constance des rapports de risque pour calculer les **résidus de Schoenfeld**... Le grand sujet (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: ::: {.box_img} $$\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 strictement 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$ (jour, mois, années...), le risque d'observer l'évènement est $e^b$ fois plus important/plus faible pour $X=1$ que pour $X=0$. ***Covariable mesurée (quantitativement)*** (mais 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 ::: {.box_img} | 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} ::: ::: {.box_img} | 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 la même lecture que lors 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 (ie 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. Les résultats apparaissent, sans connaissance médicale très poussée, conforme à ce qu'on pourrait attendre. **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`**. #### 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`. #### SAS {{< fa solid cross>}} Le modèle est estimé avec la **`proc 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 mais à faire un jour). En post-estimation, les valeurs estimées du risque pourront présenter des valeurs aberrantes , si on dévie trop de cette constance, en particulier en obtenant des négatives des taux de risque. * Important: analyser cette hypothèse revient à introduire une interaction entre les rapports de risque 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, leur croisement non tardif impliquera nécessairement un problème sur cette hypothèse. Il est donc important de les analyser en amont d'un modèle. ### Test de Grambsch-Therneau sur les résidus de Schoenfeld {#test} 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 (Sas, Stata, Python et Survival v2). T.Therneau avec la Version 3 de package `survival` a substitué - assez brutalement - un test GLS^[Moindres carrés généralisés] au test OLS. Et cela pose de gros problème, en particulier en sciences sociales. Ceci a été montré de manière très convaincante par [Shawna K. Metzger](https://www.cambridge.org/core/journals/political-analysis/article/implementation-matters-evaluating-the-proportional-hazard-tests-performance/05E2657A64B18FB1E25C2BE1F7D92C3B). On doit également 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. Ou tout simplement bien regarder les données qu'on a en face des yeux.... par exemple pour la variable *surgery*, il y a t-il besoin de faire un test? Regardez bien. ***Principe du test***: consiste à regarder la corrélation entre les **résidus de Schoenfeld** obtenus directement avec la fonction de score de la vraisemblance 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. * OLS: variance des estimateurs du modèles * GLS: calcule à chaque des variances à chaque évènement. Possible instabilité lorsque le nombre de personnes restent soumises au risque est faible. La moyenne de ces variances est égale à la variance des estimateurs du modèle. 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 à: ::: {.box_img} $$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 dans le temps, 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_0=100$) avec 50 hommes et 50 femmes. Si l'hypothèse PH est 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. En revanche, 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. Et on obtiendrait au final un rapport de risque moyen égal à 1. 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-values peuvent présenter des valeurs très faibles alors que ce n'est pas le cas pour les covariables prises une à une. Rester si possible le cas à un test à un degré de liberté, variable par variable. - **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,5)$. 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)$ ::: ::: {.box_img} ::: {.callout-important } #### TRES IMPORTANT: Pourquoi la version exacte (GLS) du test pose problème * Il sagit de la version d'origine du test (1994) mais implémenté seulement en 2020 dans la version 3 du package `survival` de R. Les autres logiciels ne l'implémente pas. * S.Metzger a effectué une décomposition complète de sa statistique dans sa version exacte [lien pour les plus acharné.e.s](https://www.cambridge.org/core/journals/political-analysis/article/implementation-matters-evaluating-the-proportional-hazard-tests-performance/05E2657A64B18FB1E25C2BE1F7D92C3B) et l'a comparé à la version appelé "approximation" (je préfère "version Ols"). Résultats: * Un terme *résiduel* se situant au numérateur de la statistique, qui traduit le degré de corrélation entre les covariables est présent. Ce qui signifie: * Qu'en présence de corrélations entre covariables, même limitées, la statistique du test augmente artificiellement, et produira tout aussi artificiellement mais mécaniquement des pvalues très faibles. Le test exact s'écarte alors fortement du test OLS, alors qu'il ne devrait être juste qu'un peu plus précis. Si les travaux de S.Metzger s'appuient bien sur des simulations, ses résultats sont vérifiés avec des données réelles qu'elle a également confrontée. C'est également ici le cas avec les données du support comme celles des TP. Pour la base utilisée dans ce support, c'est la corrélation entre les variables *surgery* et *year* qui explique les forte différences pour la deuxième variable entre les deux versions du test. * Cela signifie également que le test, qui consiste à synthétiser l'approche par l'intéraction (voir plus bas), donnera des résultats contradictoires avec cette dernière. * Dans R, et seulement dans R: * T.Therneau, proprétaire du package `survival` et visiblement discutant du travail de S.Metzger, n'a toujours pas remis le test OLS dans le package, ce qui est quand même problématique au niveau reproductibilité et réplicabilité dans le temps et dans l'espace des logiciels. C'est facheux. * Ne pas paniquer, j'ai récupéré le test OLS et on peut le charger et l'utiliser très facilement. Voir la section programmation dédiée à R pour le détail de la démarche (ou directement ce [lien](https://mthevenin.github.io/assistoolsms/R/assist/posts/coxzph/coxzph.html)). ::: ::: On ne présentera donc ici que la version OLS (idem Stata, Python, Sas {{<fa solid cross >}}. **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)$ <br> **R-Stata-Python-Sas** ::: 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. #### 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’`. #### SAS {{<fa solid cross>}} 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. ::: ### 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), c'est le principe même du modèle. Difficilement soutenable pour l’analyse des effets cliniques, il faut reconnaître que ceci peut être envisagé dans d’autres domaines, comme en sciences sociales si on ne veut pas aller plus loin. 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. La lecture des résultats ne sera pas affecté. 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. 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. Si on souhaite corriger la non proportionnalité avec une intéraction, il est inutile de complexifier le modèle pour des variables introduites comme simples contrôles... qui ne sont là que pour jouer le rôle de contrôle. Et justement la vérification de la proportionnalité peut améner à s'intérroger sur la bonne spécification du modèle. En effet, 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: Le risque de base pour chaque groupe est alors estimé en amont. ***Intéraction*** Introduire une interaction avec la durée. Cela permet d'enrichir le modèle au niveau de l’interprétation. Valable si peu de covariables présentent des problèmes d'uniformité dans le temps 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, cela facilité la lecture, mais c'est une hypothèse très forte.... car au final 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 à durée discrète/groupée (voir section suivante). 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 qui sont par définition des probabilités conditionnelles, resteront toujours dans les bonnes 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 seront présentés par des graphiques, se commentent assez difficilement. Et dont l'utilisation est devenue je crois plus que rare. **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é très 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