
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.
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ù :
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 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.
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.
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.
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.
Pour que le modèle de Cox soit valide et que ses résultats soient interprétables, deux conditions principale doivent être respectées :
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 :
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() 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) 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.
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)) 
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.
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 
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.
Les résultats du modèles de Cox peuvent être représentés comme ceci :
library(GGally)
ggcoef_model(fit, exponentiate = TRUE)

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)) 
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.
Si le sujet de la modélisation des données de survie avec R vous intéresse, je vous recommande l’ouvrage suivant
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 !
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/

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.
5 réponses
Merci beaucoup pour ce résumé très riche sur le modele de Cox.
Bonne fête de fin d’année.
Super comme d’habitude. C’est expliqué dans un langage accessible.
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
Bonjour Camille,
Faire des analyses univariées une par une est une approche classique. Par contre c’est généralement une pvalue à 0.2 qu’on utilise, pas 0.02.
Le package finalfit permet de faire toutes les analyses univariées : https://finalfit.org/articles/survival.html
Vous pouvez aussi lire cet article : https://delladata.fr/calculer-et-rapporter-les-odds-ratio-sans-se-fatiguer/
Bonne continuation
# 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 »
)