Analyses de survie et modèles de Cox

Analyses de survie et modèle de Cox

Table des matières

Introduction

Dans un précédent article, j’ai abordé les courbes de survie de Kaplan-Meier, ainsi que le test du log rank pour comparer deux courbes de survie. Comme expliqué dans cet article, le test du log rank est un test non paramétrique univarié, c’est-à-dire qui ne prend pas en compte les covariables. Or, dans de nombreuses situations, il est nécessaire de prendre en compte des covariables dans nos analyses de survie. C’est là que le modèle de Cox intervient. Ce modèle semi-paramétrique, est similaire à la régression linéaire multiple ou à la régression logistique multiple, mais il est spécialement conçu pour les données de survie.

Présentation du modèle de Cox

Le modèle de Cox est un modèle semi-paramétrique qui permet de modéliser la relation entre le temps de survie et des covariables. Autrement dit, le modèle de Cox permet d’évaluer le lien entre la survie et une, ou plusieurs variables explicatives (qualitative ou quantitative).

L’équation du modèle de Cox est :

\[ h(t|X) = h_0(t) \exp(\beta_1X_1 + \beta_2X_2 + … + \beta_pX_p)\] 

Où :

  • h(t|X) est la fonction de risque instantané de l’événement à l’instant t pour un individu avec les co-variables X
  • h_0(t) est la fonction de risque de base, c’est le risque instantané d’événement à l’instant t pour un individu avec toutes les co-variables égales à 0.
  • et où  : \[ exp(\beta_1X_1 + \beta_2X_2 + … + \beta_pX_p)\] est l’exponentielle des coefficients de régression, interprétée comme l’effet de la covariable Xi sous la forme d’un rapport de risque (Hazard Ratio). C’est un facteur multiplicatif constant qui modifie le risque de base en fonction des covariables Xi

Remarque : Le modèle de Cox ne fait aucune hypothèse sur la forme de la fonction de risque de base $h_0(t)$, c’est pourquoi on parle de modèle semi-paramétrique.

Les coefficients du modèle de Cox

Les coefficients estimés du modèle de Cox, notés beta_i, mesurent l’effet des covariables sur le risque de survenue de l’événement.

Ces coefficients correspondent aux logarithmes des Hazard Ratio (logHR). Le hazard Ratio est le ratio de risque instantané de survenue de l’événement entre deux groupes.

Pour le comprendre, envisageons un modèle de Cox simple, avec une seule variable explicative X. Supposons que nous observons deux individus avec des modalités X_1 pour l’un (par exemple « Hommes ») et X_2 pour l’autre (par exemple « Femmes »). Les risques instantanés de ces deux individus sont donnés par : 

\[  h(t|X_1) = h_0(t) \exp(\beta_1X_1)\] 

et :

\[  h(t|X_1) = h_0(t) \exp(\beta_1X_2)\]

Le rapport de risque (Hazard Ratio) entre ces deux individus est donné par : 

\[ \frac{h(t|X_1)}{h(t|X_2)} = \frac{h_0(t) \exp(\beta_1X_1)}{h_0(t) \exp(\beta_1X_2)} = \exp(\beta_1(X_1 – X_2))\]

En utilisant la transformation logarithmique, nous obtenons :

\[ \log\left(\frac{h(t|X_1)}{h(t|X_2)}\right) = \log\left(\frac{h_0(t) \exp(\beta_1 X_1)}{h_0(t) \exp(\beta_1 X_2)}\right) = \beta_1 (X_1 – X_2)\]

Il est important de noter que le hazard ratio est indépendant du temps. Cela signifie que l’effet des co-variables sur le risque est multiplicatif et constant au cours du temps. Cela illustre le principe des « risques proportionnels ».

Par exemple, si le Hazard Ratio associé à une covariable est de 2, alors le risque de l’événement pour une personne avec cette covariable sera toujours deux fois plus élevé que pour une personne sans cette covariable, et ce à n’importe quel moment de l’étude.

Les tests associés au coefficient

Comme pour n’importe quel modèle de régression, l’égalité à 0 des coefficients estimés est évaluée par un test statistique. La statistique du test, et la p-value associée, sont reportés dans les résultats du modèle de Cox.

Le test de Wald est utilisé pour tester l’hypothèse nulle selon laquelle un coefficient (logHR) est égal à zéro (c’est-à-dire qu’il n’a pas d’effet significatif). La statistique z est calculée comme le rapport du coefficient estimé $\hat{\beta}$ à son erreur standard (ES) :

\[ z = \frac{\hat{\beta}}{ES} \]


Les valeurs de z suivent une distribution normale standard sous l’hypothèse nulle. Par conséquent, une valeur absolue élevée de $z$ indique que le coefficient est significativement différent de zéro.

Si nous raisonnons en termes de Hazard Ration (HR) plutôt qu’en logHR, alors le test de Wald permet de tester l’hypothèse nulle selon laquelle le hazard ratio associé à une covariable est égal à 1.

Une p-value inférieure à un seuil alpha (généralement 0,05) suggère que le Hazard ratio significativement différent de 1, ce qui implique que la covariable a un effet significatif sur le risque de survenue de l’événement.

Interprétation des Hazard Ratios

Le Hazard Ratio (HR) est une mesure clé dans l’analyse de survie avec le modèle de Cox. Il permet d’évaluer l’effet des covariables sur le risque de survenue d’un événement.

  • HR = 1 : Lorsque le Hazard Ratio est égal à 1, cela signifie que la covariable n’a aucun effet sur le risque de survenue de l’événement. Autrement dit, la présence ou l’absence de cette covariable n’affecte pas le taux auquel l’événement se produit.
  • HR < 1 : Lorsque le Hazard Ratio est significativement inférieur à 1, cela signifie que la covariable est associée à une réduction du risque de survenue de l’événement. Par exemple, si une étude examine l’effet de l’activité physique régulière sur la survie (par rapport à l’absence d’activité physique régulière), un HR de 0,75 pour l’activité physique signifierait que l’activité physique régulière est associée à une réduction de 25% du risque de décès.
  • HR > 1 : Lorsque le Hazard Ratio est significativement supérieur à 1, cela indique que la covariable est associée à une augmentation du risque de survenue de l’événement. Par exemple, un HR de 2 signifie que le risque de survenue de l’événement est doublé pour les individus présentant cette covariable (être un homme, par exemple) par rapport à ceux qui ne la présentent pas (être une femme, par exemple).

L’interprétation des Hazard Ratios doit toujours être faite dans le contexte de l’étude et en tenant compte de la signification statistique des résultats. Un Hazard Ratio proche de 1 peut ne pas être cliniquement significatif même s’il est statistiquement significatif, et vice versa.

Les conditions d'applications

Pour que le modèle de Cox soit valide et que ses résultats soient interprétables, deux conditions principale doivent être respectées :

  • Risques Proportionnels : L’hypothèse des risques proportionnels est centrale dans le modèle de Cox. Elle stipule que le rapport des risques (Hazard Ratios) entre les groupes comparés est constant au cours du temps. Autrement dit, les effets des co variables sur le risque de survenue de l’événement sont multiplicatifs et ne varient pas dans le temps. Cette hypothèse peut être vérifiée graphiquement ou à l’aide de tests statistiques spécifiques, comme le test de Schoenfeld.
  • Linéarité des Covariables: Les covariables doivent avoir un effet linéaire sur le logarithme du risque. Cela signifie que chaque unité de changement dans une covariable est associée à une augmentation ou une diminution proportionnelle du risque de l’événement. Cette condition peut être vérifiée en examinant les résidus et en utilisant des transformations des covariables ou l’inclusion de termes polynomiaux si nécessaire.

Mise en application avec R

Dans cette partie, nous allons appliquer le modèle de Cox sur le jeu de données colon du package survival. Le jeu de données colonprovient d’un des premiers essais cliniques réussis de chimiothérapie adjuvante pour le cancer du côlon.

Les données incluent des informations sur la survie des patients, la récidive de la maladie et les décès. Chaque patient a deux enregistrements : l’un pour la récidive et l’autre pour le décès. Nous utiliserons uniquement les données concernant les décès.

La variable « rx » dans le jeu de données « colon » représente le traitement reçu par les patients. Cette variable est catégorique et comporte trois modalités différentes :

  • Observation : Aucun traitement adjuvant.
  • Levamisole : Traitement avec le levamisole, un composé à faible toxicité précédemment utilisé pour traiter les infestations de vers chez les animaux.
  • Levamisole + 5-FU : Traitement combiné avec le levamisole et le 5-FU (fluorouracile), un agent chimiothérapeutique modérément toxique.
library(survival)
library(tidyverse)

# filtre des données concernant les décès (et non la récidive)
colon2 <- colon %>% 
  filter(etype == 2) %>% 
  na.omit() 

Ajustement du Modèle de Cox

Le modèle de Cox est ajusté en utilisant la fonction coxph() du package survival. Dans ce modèle, nous allons examiner l’effet du traitement (rx) sur la survie des patients, en tenant compte des covariables suivantes : l’âge, le sexe, la perforation de la tumeur, la différenciation de la tumeur et le type de chirurgie.

La fonction Surv(time, status) crée une variable de survie en utilisant les colonnes time (temps de survie) et status (indicateur de statut de survie, par exemple, vivant ou décédé) du jeu de données colon.

fit <-coxph(Surv(time, status) ~ rx + age+sex+perfor+surg, data = colon2) 

Résultats

Coefficients

La fonction summary()permet d’obtenir les coefficients estimés du modèle de Cox, ainsi que leurs erreurs standards, les statistiques du test de Wald et les p-values associées. Il est également important de noter que les exponentielles des coefficients sont présentées. Ces valeurs correspondent aux Hazard Ratios (HR) associés à chaque covariable, ainsi que leurs intervalles de confiance à 95%.

# Obtenir les coefficients
summary(fit)

Call:
coxph(formula = Surv(time, status) ~ rx + age + sex + perfor + 
    surg, data = colon2)

  n= 888, number of events= 430 

               coef exp(coef)  se(coef)      z Pr(>|z|)   
rxLev     -0.055783  0.945744  0.113373 -0.492  0.62270   
rxLev+5FU -0.377600  0.685504  0.121373 -3.111  0.00186 **
age        0.003460  1.003466  0.004156  0.833  0.40511   
sex       -0.024593  0.975707  0.096643 -0.254  0.79913   
perfor     0.164588  1.178907  0.263106  0.626  0.53161   
surg       0.196714  1.217396  0.105480  1.865  0.06219 . 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

          exp(coef) exp(-coef) lower .95 upper .95
rxLev        0.9457     1.0574    0.7573    1.1811
rxLev+5FU    0.6855     1.4588    0.5404    0.8696
age          1.0035     0.9965    0.9953    1.0117
sex          0.9757     1.0249    0.8073    1.1792
perfor       1.1789     0.8482    0.7039    1.9744
surg         1.2174     0.8214    0.9900    1.4970

Concordance= 0.552  (se = 0.014 )
Likelihood ratio test= 15.89  on 6 df,   p=0.01
Wald test            = 15.48  on 6 df,   p=0.02
Score (logrank) test = 15.61  on 6 df,   p=0.02 

La fonction parameters() permet d’obtenir une présentation plus synthétique.

library(parameters)
parameters(fit, exponentiate = TRUE)

Parameter    | Coefficient |       SE |       95% CI |     z |     p
--------------------------------------------------------------------
rx [Lev]     |        0.95 |     0.11 | [0.76, 1.18] | -0.49 | 0.623
rx [Lev+5FU] |        0.69 |     0.08 | [0.54, 0.87] | -3.11 | 0.002
age          |        1.00 | 4.17e-03 | [1.00, 1.01] |  0.83 | 0.405
sex          |        0.98 |     0.09 | [0.81, 1.18] | -0.25 | 0.799
perfor       |        1.18 |     0.31 | [0.70, 1.97] |  0.63 | 0.532
surg         |        1.22 |     0.13 | [0.99, 1.50] |  1.86 | 0.062

Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed
  using a Wald z-distribution approximation. 

Ici, les résultats mettent en évidence un effet significatif du traitement (rx) sur la survie des patients. Les patients recevant le traitement Levamisole + 5-FU ont un risque de décès significativement plus faible que les patients qui ne recoivent aucun traitement (Obs), avec un Hazard Ratio (HR) de 0,70 (IC 95% : 0,55 – 0,88, p =0,002). Cela signifie que le risque de décès est réduit de 30% pour les patients recevant le traitement Levamisole + 5-FU par rapport à ceux ne recevant aucun traitement, après ajustement pour les autres covariables.

Cependant, pour avoir confiance dans ces résultats, il est important de vérifier si les conditions d’application du modèle de Cox sont respectées.

Evaluation des conditions d’application

Risques Proportionnels

Pour vérifier si l’hypothèse des risques proportionnels est respectée, nous pouvons utiliser le test de Schoenfeld. Ce test vérifie si les résidus partiels sont indépendants du temps. La fonction cox.zph() du package survival permet de réaliser ce test pour chacune des covariables incluses dans le modèle.

cox.zph(fit)

         chisq df    p
rx     2.05794  2 0.36
age    0.77111  1 0.38
sex    0.14544  1 0.70
perfor 0.75925  1 0.38
surg   0.00919  1 0.92
GLOBAL 3.71006  6 0.72 

Ici, aucune p-value n’est inférieure à 0,05, ainsi rien ne permet de conclure que l’hypothèse des risques proportionnels n’est pas respectée pour aucune des covariables incluses dans le modèle.

Nous pouvons également visualiser graphiquement les résidus partiels en fonction du temps pour chaque covariable en utilisant la fonction plot() sur les résultats de cox.zph().

par(mfrow=c(2,3)) # division de la fen^être graphique en 2 lignes et 3 colonnes
plot(cox.zph(fit)) 
Evaluer graphiquement l'hypothèse des risques proportionnels

Si l’hypothèse des risques proportionnels est respectée, les résidus partiels ne devraient pas montrer de tendance significative en fonction du temps.Ils devraient donc être globalement plats, comme c’est le cas ici. Si les résidus partiels montrent une tendance significative en fonction du temps, cela suggère que l’hypothèse des risques proportionnels n’est pas respectée. Dans cette situtation, il est possible de stratifier les covariables qui ne respectent pas l’hypothèse des risques proportionnels, ou de modéliser explicitement leur effet en fonction du temps.

L’hypothèse de linéarité

Cette hypothèse consiste à évaluer le lien linéaire entre les covariables numériques et le logarithme du risque.

Elle est généralement évaluer de manière visuelle en traçant les résidus de martingale du modèle nul, en fonction des covariables. Les résidus de martingale sont une forme de résidus utilisée dans les modèles de Cox. Ils mesurent la différence entre le nombre d’événements observés et le nombre d’événements prédits par le modèle pour chaque individu. Si la relation entre les covariables et le log-risque est linéaire, les résidus de martingale devraient être répartis aléatoirement autour de zéro lorsque tracés contre les covariables. Si une tendance non linéaire est observée, cela suggère une violation de l’hypothèse de linéarité.

# ajustement du modèle nul
result.0.coxph <- coxph(Surv(time, status) ~ 1, data = colon2)

# stockage des résidus de Martinagle
rr.0 <- residuals(result.0.coxph, type = "martingale")

# Tracé des résidus de Martingale en fonction de l'âge
plot(rr.0~colon2$age 
Evaluer graphiquement l'hypothèse de linéarité dans les analyses de survie avec les résidus de martingale

Ici, l’hypothèse de log-linéarite ne semble pas être rejetée.

Remarque: En principe il est possible d’employer la fonction ggfunctionnal() du package survminer pour visualiser les résidus de martingale en fonction des covariables. Cependant, elle me renvoie une erreur liée aux dimensions que je n’ai pas réussi à résoudre.

Plot des résultats

Les résultats du modèles de Cox peuvent être représentés comme ceci :

library(GGally)
ggcoef_model(fit, exponentiate = TRUE)
 
Représenter graphiquement les coefficients du modèle de Cox

Courbes de survie

Il est également possible de tracer les courbes de survie ajustées pour chaque groupe de traitement en utilisant la fonction survfit() du package survival, puis de les visualiser à l’aide de ggsurvfit() du package ggsurvfit.

library(survival)
library(ggsurvfit)
coxph(Surv(time, status) ~ age + sex + perfor  + surg + strata(rx), data = colon2) %>% 
  survfit2() %>%
  ggsurvfit() +
  add_confidence_interval() +
  add_risktable() +
  scale_y_continuous(limits = c(0, 1)) 
Afficher les courbes de survie ajustées

Comparaisons

Il semble possible de comparer les groupes de traitements en utilisant la fonction emmeans() du package emmeans. Cette fonction permet de calculer les moyennes estimées pour chaque groupe, ainsi que les comparaisons entre les groupes.

À noter : c’est la première fois que j’utilise cette fonction, et il s’agit donc d’un test préliminaire. Soyez prudent dans l’interprétation de ces résultats.

library(emmeans)
emmeans(fit, specs=pairwise~ rx)

$emmeans
 rx        emmean    SE  df asymp.LCL asymp.UCL
 Obs      0.37529 0.293 Inf    -0.198     0.949
 Lev      0.31951 0.312 Inf    -0.292     0.931
 Lev+5FU -0.00231 0.314 Inf    -0.617     0.613

Results are averaged over the levels of: sex, perfor, surg 
Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

$contrasts
 contrast        estimate    SE  df z.ratio p.value
 Obs - Lev         0.0558 0.113 Inf   0.492  0.8751
 Obs - (Lev+5FU)   0.3776 0.121 Inf   3.111  0.0053
 Lev - (Lev+5FU)   0.3218 0.124 Inf   2.601  0.0252

Results are averaged over the levels of: sex, perfor, surg 
Results are given on the log (not the response) scale. 
P value adjustment: tukey method for comparing a family of 3 estimates  

La première partie correspond aux moyennes marginales estimées du log-risque pour chaque groupe de traitement, tandis que la seconde partie présente les comparaisons entre les groupes.

Ces résultats montrent que le traitement combiné Lev+5FU est associé à un risque réduit par rapport aux autres traitements (Obs et Lev), lorsque les résultats sont ajustés pour les autres variables du modèle.

Pour aller plus loin :

Si le sujet de la modélisation des données de survie avec R vous intéresse, je vous recommande l’ouvrage suivant

Le mot de la fin

J’espère que cet article permettra au plus grand nombre de découvrir et/ou de mieux comprendre les modèles de Cox et leur application en analyse de survie. N’hésitez pas à laisser un commentaire pour partager vos questions ou expériences, cela nous aide à progresser ensemble !

Poursuivez votre lecture

Vous souhaitez vous former en biostatistiques ?

Retrouvez les informations de ma formation de remise à niveau en biostatistiques avec R ici :

👉 https://delladata.fr/formation-remise-niveau-biostats-r/

Et les informations de ma formation de Remise à niveau en biostatistiques avec Jamovi et Jasp là
👉
https://delladata.fr/formation-remise-niveau-biostats-jamovi/

Vous souhaitez soutenir mon travail ?​

5 réponses

  1. Merci pour ce sujet, très facile à suivre.
    J’ai une question sur le plan méthodo: je suis dans le milieu de la cancérologie donc biensur les analyses de survie sont le gros de mon sujet.

    Pour choisir les variables à faire entrer dans une analyse de Cox multivariée comme ci-dessus, est ce correct de realiser les analyses de Cox univariées une par une (et choisir celle avec une p value < 0,02 ou autre) ? si oui y a t'il un moyen "rapide" de faire ces multiples analyses sous R (plutot que 1 à 1)?
    et alors quelle est la place du Log rank dans tout cela?

    Merci beaucoup par avance pour votre aide précieuse
    Camille

  2. # visualiser les résidus de martingale en fonction des covariables
    survminer::ggcoxfunctional(
    fit,
    data = colon2,
    # x = « agel »,
    # y = « status »,
    fun = « martingale »,
    funlabel = « Résidus de Martingale »
    )

Laisser un commentaire

Votre adresse e-mail ne sera pas publiée. Les champs obligatoires sont indiqués avec *

Fonctions statistiques R

Aide mémoire off'R ;)

Enregistrez vous pour recevoir gratuitement mes fiches « aide mémoire » (ou cheat sheets) qui vous permettront de réaliser facilement les principales analyses biostatistiques avec le logiciel R et pour être informés des mises à jour du site.