Analyses de survie

L’analyse de survie est centrale dans les essais cliniques, où de nombreuses mesures concernent des “temps avant événements”. De nombreux essais, particulièrement dans le domaine de la cancérologie, s’intéressent par exemple à :

• la survie globale (OS : Overall Survival),
• ou à la survie sans progression (PFS : Progression-Free Survival).

L’objectif de ce TD est d’aborder les méthodes de bases utilisées pour analyser ces “temps d’attente”, quels qu’ils soient.

On aura besoin des packages suivants :

# comme toujours, nous utiliserons dans ce tutoriel la syntaxe du "tidyverse"
library(tidyverse)
# nous utiliserons "cowplot" pour faire des figures multiples facilement
library(cowplot)
# ainsi que les deux packages suivants dédiés aux analyses de survie
library(survival)
library(survminer)
knitr::opts_chunk$set(echo = T, fig.width=10, warning=F, message=F)

Exemple du TD

On prend ici un exemple fictif d’essai clinique randomisé contrôlé dans le cancer du pancréas, en deux groupes, l’un recevant un traitement et l’autre recevant un placebo.

1. L’objectif principal consiste à montrer la supériorité du traitement sur la survie globale des patients.
2. L’objectif secondaire consiste à montrer la supériorité du traitement sur la survie sans progression.

L’inclusion des patients s’est faite sur une durée de deux ans. A l’issue de la période d’inclusion, l’étude a duré 3 ans. Certains patients peuvent être suivis 5 ans (premiers inclus) ou 3 (derniers inclus). Mais certains patients peuvent également être perdus de vus (peut-être ont-ils déménagé ailleurs, peut-être ne répondent-ils plus au téléphone…).

Pour chaque patient, un certain nombre de temps sont enregistrés dans la base de données :

d <- read.csv("data_survival.csv")
head(d)

##   id     arm age date_rando date_death  date_prog date_last_follow_up
## 1  1 control  38 2018-09-08 2019-12-15       <NA>          2019-12-15
## 2  2 control  53 2018-01-20 2020-03-06 2019-03-02          2020-03-06
## 3  3 control  58 2018-02-08       <NA> 2019-05-05          2023-01-15
## 4  4 control  54 2019-01-11 2019-09-22 2019-07-14          2019-09-22
## 5  5 control  44 2018-05-25 2020-06-18       <NA>          2020-06-18
## 6  6 control  45 2018-06-01 2019-07-23       <NA>          2019-07-23

Prenez 3 min pour explorer la façon dont ce jeu de données est construit. Dans quels cas a-t-on une date de décès inconnue ? Quelle est la date de début de l’essai ? De clôture de l’essai ?

Le package lubridate et ses fonctions ymd(), dmy() ou mdy() est bien utile pour interpréter les chaînes de caractères comme des dates, de façon à pouvoir les manipuler pour en tirer des durées entre deux dates.

d <- d %>%
    mutate(
        date_death = ymd(date_death),
        date_rando = ymd(date_rando),
        date_last_follow_up = ymd(date_last_follow_up),
        date_prog = ymd(date_prog)
    )

Premiers essais naïfs

Avant de s’intéresser à ce qu’on peut faire de plus fancy à partir de ces temps d’attente, réfléchissons à ce que nous serions déjà en mesure de faire avec les outils que nous connaissons déjà.

Proposez ci-dessous une analyse de la probabilité de décès à 6 mois entre les deux groupes.

d <- d %>% mutate(
    time_before_death = as.numeric(date_death - date_rando),
    is_alive_6m = (date_last_follow_up - date_rando) > 365/2,
    is_alive_24m = (date_last_follow_up - date_rando) > 365*2,
)

table(d$arm, d$is_alive_6m)

##            
##             FALSE TRUE
##   control      20  280
##   treatment    14  286

prop.test(table(d$arm, d$is_alive_6m))

## 
##  2-sample test for equality of proportions with continuity correction
## 
## data:  table(d$arm, d$is_alive_6m)
## X-squared = 0.77946, df = 1, p-value = 0.3773
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  -0.02029847  0.06029847
## sample estimates:
##     prop 1     prop 2 
## 0.06666667 0.04666667

Puis une analyse de la probabilité de décès à 2 ans entre les deux groupes.

table(d$arm, d$is_alive_24m)

##            
##             FALSE TRUE
##   control     118  182
##   treatment    96  204

prop.test(table(d$arm, d$is_alive_24m))

## 
##  2-sample test for equality of proportions with continuity correction
## 
## data:  table(d$arm, d$is_alive_24m)
## X-squared = 3.2032, df = 1, p-value = 0.07349
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  -0.006432081  0.153098747
## sample estimates:
##    prop 1    prop 2 
## 0.3933333 0.3200000

Puis une analyse des temps d’attente avant décès avec un test de Wilcoxon.

wilcox.test(time_before_death ~ arm, data=d)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  time_before_death by arm
## W = 18038, p-value = 0.3132
## alternative hypothesis: true location shift is not equal to 0

On a deux reproches principaux à faire à ces approches naïves :

• d’une part, si on analyse des proportions de survivants, il nous faut définir UN temps auquel on regarde la survie, ce qui est loin d’être évident à décider avant de faire l’étude. Répéter l’analyse pour un grand nombre de temps différents serait une très mauvaise pratique amenant à une inflation du risque de première espèce.
• d’autre part, on perd toute l’information à propos des durées de survie des patients qui ont été suivi un peu mais pas suffisamment pour qu’on ait observé l’événement. C’est cette information qu’on va chercher à capter grâce à la notion de censure qu’on introduit juste après.

Censure

Soit \(T_i :=\) le temps d’attente avant décès de l’individu \(i\).

Lorsqu’on traite de données de temps d’attente dans la vraie vie, on se retrouve souvent dans des cas où on ne peut pas observer exactement le temps de survenue d’un événément, mais où on possède malgré tout une information parcellaire à propos de ce temps. On dit que nos données sont censurées (à droite/à gauche).

• censure à gauche : si on sait que \(T_i < t_i^+\).
• censure à droite : si on sait que \(T_i > t_i^-\).
• censure à gauche et à droite : si on sait que \(T_i \in (t_i^-, t_i^+)\).

Les biais de définition de t0

Dans quel cas se trouve-t-on ici avec nos perdus de vue ?

La plupart du temps, dans les essais cliniques, nous aurons des données qui correspondent à des censures à droite, c’est à dire à des temps de follow-up sans événement, qui nous apportent tout de même une information énorme sur ce qu’il ne s’est pas passé jusqu’à ce temps.

Toutes les méthodes qui suivent vont donc parvenir à gérer cette censure à droite.

Survie et estimation de Kaplan-Meier

Le problème de l’estimation de la courbe de survie consiste à estimer la courbe :

\[ S(t) = \mathbb{P}(T > t) \]

à partir de données potentiellement censurées à droite.

Le principe de l’approche de Kaplan-Meier est le suivant :

• le temps est découpé en sous-intervalles pendant lesquels il ne se passe aucun décès.
• sur un intervalle \(i\) avec \(n_i\) individus susceptibles de décéder au début de l’intervalle et \(d_i\) décès à la fin de l’intervalle, la probabilité de survivre est estimée comme \((1 - \frac{d_i}{n_i})\).
• donc la survie jusqu’à ce temps correspond au produit des probabilités de survie de tous les intervalles précédents.

La formule d’estimation de la courbe de survie par la méthode de Kaplan-Meier s’écrit donc de façon plus concise comme :

\[ \hat{S}(t) = \prod_{i: t_i < t} \left( 1 - \frac{d_i}{n_i} \right) \]

La censure se prend en compte très naturellement ici en réduisant \(n_i\) à chaque fois qu’un individu est perdu de vue. Un “tick” est alors visuellement noté sur la courbe de survie pour signifier qu’un individu est perdu sans être décédé.

Avec R, on utilise la fonction survfit et ggsurvplot qui prennent en entrée :

• un temps avant survenue d’un événément (time)
• une valeur binaire 1 si un décès survient ou 0 si une censure survient à ce temps (event)
• une éventuelle covariable correspondant à des groupes pour lesquels on souhaite tracer des courbes de survie différentes.

Préparez le jeu de données pour avoir time et event proprement définis (regardez l’aide de Surv) et utilisez le morceau de code qui suit pour tracer les courbes de survie de Kaplan-Meier.

Sentez vous libres de jouer avec les différentes options listées. Les amateurs de ggplot2 pourront accéder au graphe ggplot2 modifiable avec l’objet final$plot.

d <- d %>% mutate(
    time = as.numeric(date_last_follow_up - date_rando),
    event = ifelse(!is.na(date_death), 1, 0)
)

km_os <- survfit(Surv(time, event) ~ arm, data=d)

graph_os <- ggsurvplot(
    km_os,
    pval = F, 
    conf.int = FALSE,
    risk.table = F, # Add risk table
    risk.table.col = "strata", # Change risk table color by groups
    #linetype = "strata", # Change line type by groups
    surv.median.line = "hv", # Specify median survival
    ggtheme = theme_bw() # Change ggplot2 theme
)

graph_os

Notez que ce graphique repose sur plusieurs hypothèses qui seront également implicites dans la suite :

1. On suppose que les individus qui sont censurés ont les mêmes chances de survie que les autres, i.e. que la censure n’est pas informative d’un futur décès,
2. On suppose que les survies sont identiques pour les sujets recrutés tôt et tard dans l’étude, de façon à pouvoir ré-aligner tout le monde sur un t0 qui correspond à la randomisation.

Il nous reste un dernier item à discuter en lien avec ces analyses de survie, la notion de médiane de survie qui est une quantité pratique à se représenter que tout le monde aimerait partager. Oui mais voilà, même une notion aussi simple s’avère être propice à faire n’importe quoi avec nos problèmes de temps d’attente !

Testez ce que ça donne de prendre naïvement la médiane des temps de décès.

Puis testez ce que ça donne de prendre naïvement la médiane des temps de follow-up.

d %>% group_by(arm) %>% 
    summarise(
        median_t_death = median(time_before_death, na.rm=T),
        median_t_follow = median(time)
    )

## # A tibble: 2 × 3
##   arm       median_t_death median_t_follow
##   <chr>              <dbl>           <dbl>
## 1 control              714            871 
## 2 treatment            776           1132.

Réfléchissez à la signification de ces deux quantités. En réalité, ce qu’on cherche quand on s’intéresse à la médiane de survie, c’est la valeur \(t_m\) telle que \(S(t_m) = 1/2\).

On l’obtient avec le code qui suit :

surv_median(km_os)

##          strata median lower upper
## 1   arm=control    902   836  1072
## 2 arm=treatment   1213  1145  1319

Hazard et hazard ratio

Le “hazard” \(h(t)\), parfois aussi appelé “hazard rate” est le taux de survenue de l’événement au temps \(t\) sachant qu’on a survécu jusqu’au temps \(t\).

\[ h(t) := lim_{dt \rightarrow 0} ~ \frac{\mathbb{P}(T \in (t, t+dt) | T > t)}{dt} \]

Le hazard est positif en tout temps \(t\) mais n’est pas nécessairement monotone.

Si on parle tant de cette quantité en analyse de survie, c’est que les principaux tests statistiques proposés dans le domaine pour comparer les survies de deux groupes 1 et 2 reposent sur l’hypothèse suivante :

\[ \exists K \in \mathbb{R}^+, \forall t, ~ \frac{h_1(t)}{h_2(t)} = K \]

C’est l’hypothèse des hazard proportionnels ou du hazard ratio constant.

Notez qu’on ne fait pas d’hypothèse forte sur la forme de la fonction du hazard, mais qu’on a besoin de faire l’hypothèse que cette forme est similaire dans les deux groupes d’intérêt.

Pour se fixer un peu plus les idées, voici deux illustrations de relation entre hazards et survie dans le cadre de deux lois théoriques différentes.

Hazards constants et loi exponentielle

Pour voir ce que ça donne avec une distribution que vous avez des chances de connaître, la distribution exponentielle. Comme vous l’avez sans doute déjà entendu dire, la loi exponentielle est “sans mémoire”, c’est à dire que le taux de survenue d’un événément, sachant que l’événement n’est pas encore survenu (le hazard) est constant ! Voici donc une première illustration à se mettre en tête :

t <- seq(0,10,0.1)
rates <- c(0.5,1,2)
distrib_exp <- data.frame(
    t = rep(t, each=length(rates)),
    density = dexp(rep(t,each=length(rates)), rate=rep(rates, length(t))),
    survival = pexp(rep(t,each=length(rates)), rate=rep(rates, length(t)), lower.tail=F),
    lambda = as.factor(rep(rates, length(t))),
    hazard = rep(rates, length(t))
)

p1 <- distrib_exp %>% ggplot(aes(x=t, y=density, color=lambda)) +
    xlab("t") +
    ylab("Densité : P(t < T < t+dt)/dt") +
    geom_line() +
    theme_bw() +
    theme(legend.position = "none") + 
    scale_color_manual(values=c("#548687", "#ff5555", "#8fbc94"))

p2 <- distrib_exp %>% ggplot(aes(x=t, y=survival, color=lambda)) +
    xlab("t") +
    ylab("Survie : P(T>t)") +
    geom_line() +
    theme_bw() +
    theme(legend.position = "none") + 
    scale_color_manual(values=c("#548687", "#ff5555", "#8fbc94"))

p3 <- distrib_exp %>% ggplot(aes(x=t, y=hazard, color=lambda)) + 
    xlab("t") +
    ylab("Hazard : P(t < T < t+dt | T > t)/dt") +
    geom_line() +
    theme_bw()  + 
    scale_color_manual(values=c("#548687", "#ff5555", "#8fbc94"))

plot_grid(p1, p2, p3, ncol = 3, nrow = 1)

Hazards d’une loi de Weibull

Histoire de ne pas retenir que ce qui se passe dans le cas archi-simple de la loi exponentielle, voici une seconde loi très employée dans le domaine de l’analyse de survie : la loi de Weibull.

Elle a une propriété rigolote, qui est uniquement propre à cette famille de distribution : on peut l’utiliser à la fois dans le cadre d’une hypothèse de hazards proportionnels et dans le cadre d’une hypothèse d’accelerated failure time, c’est à dire que la structure de la distribution est telle qu’une augmentation de hazard d’un certain facteur accélère la survie d’un autre facteur.

Remarquez aussi que la loi exponentielle appartient à la famille des loi de Weibull, et que la même propriété était donc vraie sur le graphe ci-dessus. Simplement à présent on peut spécifier des hazards non-constants dans le cadre d’une loi de Weibull.

t <- seq(0,5,0.01)
shape <- 1.6
scales <- c(1, 2)
distrib_weibull <- data.frame(
    t = rep(t, each=length(scales)),
    density = dweibull(rep(t,each=length(scales)), shape=shape, scale=rep(scales, length(t))),
    survival = pweibull(rep(t,each=length(scales)), shape=shape, scale=rep(scales, length(t)), lower.tail=F),
    scale = as.factor(rep(scales, length(t))) 
)
distrib_weibull <- distrib_weibull %>% mutate(hazard = density / survival)

p1 <- distrib_weibull %>% ggplot(aes(x=t, y=density, color=scale)) +
    xlab("t") +
    ylab("Densité : P(t < T < t+dt)/dt") +
    geom_line() +
    theme_bw() +
    theme(legend.position = "none") + 
    scale_color_manual(values=c("#548687", "#ff5555", "#8fbc94"))

p2 <- distrib_weibull %>% ggplot(aes(x=t, y=survival, color=scale)) +
    xlab("t") +
    ylab("Survie : P(T > t)") +
    geom_line() +
    theme_bw() +
    theme(legend.position = "none",
          plot.title = element_text(hjust=0.5)) +
    ggtitle("S_1(t) = S_0(0.5 x t)") + 
    scale_color_manual(values=c("#548687", "#ff5555", "#8fbc94"))

p3 <- distrib_weibull %>% ggplot(aes(x=t, y=hazard, color=scale)) + 
    xlab("t") +
    ylab("Hazard : P(t < T < t+dt | T > t)/dt") +
    geom_line() +
    theme_bw() + 
    labs(color = "sigma") +
    theme(plot.title = element_text(hjust=0.5)) +
    ggtitle("h_1(t) = 0.33 h_0(t)") + 
    scale_color_manual(values=c("#548687", "#ff5555", "#8fbc94"))

plot_grid(p1, p2, p3, ncol = 3, nrow = 1)

Test du logrank

Le test du logrank permet de répondre à la question suivante :

• H0 : HR = 1
• H1 : HR \(\neq\) 1.

La principale condition d’application pour le réaliser consiste à supposer que nos hazards sont proportionnels.

La syntaxe R pour le réaliser suit.

survdiff(Surv(time, event) ~ arm, data=d)

## Call:
## survdiff(formula = Surv(time, event) ~ arm, data = d)
## 
##                 N Observed Expected (O-E)^2/E (O-E)^2/V
## arm=control   300      213      182      5.14      9.63
## arm=treatment 300      180      211      4.45      9.63
## 
##  Chisq= 9.6  on 1 degrees of freedom, p= 0.002

Que concluez-vous de ce test ?

Modèle de Cox

Si \((X_1^{(i)}, X_2^{(i)}, ..., X_p^{(i)})\) sont \(p\) variables explicatives mesurées pour l’individu \(i\), le modèle de Cox consiste à relier le hazard de l’individu \(i\) avec ces covariables selon la forme suivante :

\[ h_i(t) = h_0(t) e^{\sum_{j=1}^p \beta_j X_j^{(i)}} \]

Comme pour d’autres modèles de régression que vous connaissez déjà, vous pouvez imaginer inclure dans ce modèle tout type de variable :

• avec une variable numérique, l’exponentielle du paramètre correspond au facteur multiplicatif du hazard ratio associé à une augmentation de +1 unité de la variable.
• avec une variable catégorielle, on a un paramètre par modalité (moins la modalité de référence), et l’exponentielle du paramètre correspond au facteur multiplicatif du hazard ratio associé au fait d’avoir cette modalité.

Le code suivant vous donne la syntaxe permettant de fitter le modèle de Cox sur nos données, pour réaliser un test similaire au test du logrank présenté précédemment. Notez que ces deux tests reposent sur des théorèmes différents, mais qu’ils sont “asymptotiquement équivalents”, c’est à dire qu’ils devraient vous donner des p-values similaires pour un nombre d’observations assez grand. En pratique, pour éviter toute multiplication de tests, on décide AVANT de regarder les données quel test on réalisera, on l’écrit dans le plan d’analyse statistique, puis on l’applique.

coxph(Surv(time, event) ~ arm, data=d)

## Call:
## coxph(formula = Surv(time, event) ~ arm, data = d)
## 
##                 coef exp(coef) se(coef)      z     p
## armtreatment -0.3134    0.7309   0.1014 -3.091 0.002
## 
## Likelihood ratio test=9.6  on 1 df, p=0.001949
## n= 600, number of events= 393

Obtient-on un résultat cohérent avec le test du logrank dans notre cas ?

Modifiez le code précédent pour rajouter des covariables au modèle. Observez-vous un effet de l’âge des patients sur le hazard ?

coxph(Surv(time, event) ~ arm + age, data=d)

## Call:
## coxph(formula = Surv(time, event) ~ arm + age, data = d)
## 
##                   coef exp(coef)  se(coef)      z       p
## armtreatment -0.315136  0.729690  0.101444 -3.107 0.00189
## age           0.004706  1.004718  0.007534  0.625 0.53215
## 
## Likelihood ratio test=9.99  on 2 df, p=0.006781
## n= 600, number of events= 393

Evénements intercurrents et PFS

Vous souhaitez à présent vous attaquer à la deuxième question d’intérêt de l’étude : existe-t-il un impact du traitement sur l’évolution des tumeurs ? On souhaite donc regarder \(T_i' :=\) le temps avant progression de l’individu \(i\).

MAIS : si on fait cela, on va avoir deux problèmes :

• d’une part certains individus sont perdus de vue avant de progresser, ce qu’on est prêt à prendre en compte en censurant à droite.
• d’autre part certains individus risquent de décéder AVANT de progresser.

Lorsqu’un événement (ici la mort) empêche d’observer un autre événement (la progression), on parle d’événements intercurrents. Il existe plusieurs façons de prendre en compte ce problème, mais on ne s’intéressera ici qu’à la plus simple : lorsque nos deux événements sont tous deux négatifs (ou positifs) pour le patients, le plus simple est de définir notre temps d’intérêt comme étant le premier temps auquel l’un des deux événements survient.

Ici, on définit classiquement en cancérologie la PFS comme étant la survie sans progression (ni décès). Et le temps d’intérêt associé est donc le décès ou la progression : le premier des deux qui survient.

Estimez la PFS à l’aide d’un graphe de Kaplan-Meier.

Comparez la PFS dans les deux bras à l’aide d’un test de logrank ou d’un modèle de Cox.

Que concluez-vous ?

d <- d %>% mutate(
    time_prog_or_death = as.numeric(ifelse(!is.na(date_prog), date_prog - date_rando, date_last_follow_up - date_rando)),
    event_prog_or_death = ifelse(!is.na(date_prog) | !is.na(date_death), 1, 0)
)

km_pfs <- survfit(Surv(time_prog_or_death, event_prog_or_death) ~ arm, data=d)

graph_pfs <- ggsurvplot(
    km_pfs,
    pval = F, 
    conf.int = FALSE,
    risk.table = F, # Add risk table
    risk.table.col = "strata", # Change risk table color by groups
    #linetype = "strata", # Change line type by groups
    surv.median.line = "hv", # Specify median survival
    ggtheme = theme_bw() # Change ggplot2 theme
)

cox_pfs <- coxph(Surv(time_prog_or_death, event_prog_or_death) ~ arm, data=d)

Conclusion/Discussion

Dans les essais cliniques, la norme pour s’intéresser à un temps d’attente avant événement sera souvent le duo représentation graphique de Kaplan-Meier et modèle de Cox associé. L’utilisation de ces deux concepts stats est simple, mais il faudra toujours bien prendre garde :

• à bien spécifier les événements considérés comme censure,
• à bien réfléchir aux éventuels événements inter-currents et à la définition de notre temps d’intérêt.

Dans les études observationnelles, il est également possible de faire de l’analyse de survie. Les exemples ne manquent d’ailleurs pas dans la litérature.

MAIS ! Il vous faudra garder à l’esprit deux problèmes majeurs dans ces cas-là. Contrairement à l’essai clinique, la définition du \(t_0\) (randomisation dans l’essai clinique) est difficile à réaliser de façon “juste et équitable” pour les deux groupes. Les deux principaux problèmes rencontrés sont les suivants :

1. Si on définit les individus d’un groupe de traitement comme étant tous ceux qui ont déjà commencé le traitement avant une certaine date, alors on s’expose à un biais de déplétion des susceptibles. Typiquement, imaginez des fumeurs et des non-fumeurs. Les individus susceptibles de faire un cancer du poumon dans le groupe fumeur ont peut-être déjà tous fait leur cancer, ce qui risque de diminuer le nombre de cancers du poumon observés dans le groupe fumeur.

2. Si maintenant on souhaite éviter ce problème et qu’on définit quelqu’un dans le groupe traitement comme quelqu’un qui n’avait pas commencé le traitement à \(t_0\), mais le commence plus tard. Alors on obtient, cette fois, un biais de temps immortel, qui correspond au fait que les patients de ce groupe traitement ne peuvent pas décéder (par construction) avant la date de traitement qui les fait appartenir à leur groupe.

Un nombre é-nor-me d’études observationnelles sont concernées par l’un, l’autre, ou les deux problèmes à la fois, et plusieurs exemples historiques de désaccords entre études observationnelles et essais cliniques reposent sur ces deux problème de \(t_0\). Gardez bien ça en tête pour critiquer les futurs articles que vous lirez !