Dans le cadre de ce TD, nous allons utiliser les données veteran disponible dans le package survival

Importer le jeu de données veteran dans le package survival

#installer survival d'abord
#install.packages("survival")
library(survival)

data(package = "survival", veteran)
## Warning in data(package = "survival", veteran): data set 'veteran' not found

Quelle est la structure des données et mettre au bon format les variables qui ne sont pas

str(veteran)
## 'data.frame':    137 obs. of  8 variables:
##  $ trt     : num  1 1 1 1 1 1 1 1 1 1 ...
##  $ celltype: Factor w/ 4 levels "squamous","smallcell",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ time    : num  72 411 228 126 118 10 82 110 314 100 ...
##  $ status  : num  1 1 1 1 1 1 1 1 1 0 ...
##  $ karno   : num  60 70 60 60 70 20 40 80 50 70 ...
##  $ diagtime: num  7 5 3 9 11 5 10 29 18 6 ...
##  $ age     : num  69 64 38 63 65 49 69 68 43 70 ...
##  $ prior   : num  0 10 0 10 10 0 10 0 0 0 ...

Réponse

trt est codé en numéric alors qu’elle devrait être qualitative binaire

Recodons cette variable au bon format

veteran$traitement <- factor(
  veteran$trt,
  levels = c(1, 2),
  labels = c("standard", "test"))

Présentez quelques statistiques descriptives de ces données

Réponse

summary(veteran)
##       trt             celltype       time           status      
##  Min.   :1.000   squamous :35   Min.   :  1.0   Min.   :0.0000  
##  1st Qu.:1.000   smallcell:48   1st Qu.: 25.0   1st Qu.:1.0000  
##  Median :1.000   adeno    :27   Median : 80.0   Median :1.0000  
##  Mean   :1.496   large    :27   Mean   :121.6   Mean   :0.9343  
##  3rd Qu.:2.000                  3rd Qu.:144.0   3rd Qu.:1.0000  
##  Max.   :2.000                  Max.   :999.0   Max.   :1.0000  
##      karno          diagtime           age            prior          traitement
##  Min.   :10.00   Min.   : 1.000   Min.   :34.00   Min.   : 0.00   standard:69  
##  1st Qu.:40.00   1st Qu.: 3.000   1st Qu.:51.00   1st Qu.: 0.00   test    :68  
##  Median :60.00   Median : 5.000   Median :62.00   Median : 0.00                
##  Mean   :58.57   Mean   : 8.774   Mean   :58.31   Mean   : 2.92                
##  3rd Qu.:75.00   3rd Qu.:11.000   3rd Qu.:66.00   3rd Qu.:10.00                
##  Max.   :99.00   Max.   :87.000   Max.   :81.00   Max.   :10.00

Construisons l’histogramme et le box plot de l’âge et de la variable time

par(mfrow = c(2, 2))
hist(veteran$age, 
     main = "age", col = "blue")
hist(veteran$time, , 
     main = "Temps (en jours)",
     col = "red")
boxplot(veteran$age, 
        main = "age", col = "blue")
boxplot(veteran$time, main = "Temps",
        col = "red")

par(mfrow = c(1,1))
#effectifs 
table(veteran$celltype)
## 
##  squamous smallcell     adeno     large 
##        35        48        27        27
#proportion 
(prp <- prop.table(table(veteran$celltype)))
## 
##  squamous smallcell     adeno     large 
## 0.2554745 0.3503650 0.1970803 0.1970803
par(mfrow = c(1,3))
barplot(table(veteran$celltype), 
        main = "effectifs", 
        col = "pink")
pie(prp, main = "proportion")
barplot(prp, 
        main = "proportion",
        col = "purple")

par(mfrow = c(1,1))

Pour la variable status et la variable

table(veteran$status)
## 
##   0   1 
##   9 128
prop.table(table(veteran$status))
## 
##          0          1 
## 0.06569343 0.93430657
cat(round(prop.table(table(veteran$status))[2]*100, 2),"% des individus ont présenté l'évènement.")
## 93.43 % des individus ont présenté l'évènement.

Analyse de survie univariée

Estimateur de Kaplan-Meier

On utilisera la fonction survfit(). Lire l’aide de cette fonction ?survfit

formula <- Surv(time, status)~1
fit <- survfit(formula, data = veteran)

Les sorties de la fonction survfit peuvent être obtenues avec les fonctions print ou summary.

print(fit)
## Call: survfit(formula = formula, data = veteran)
## 
##        n events median 0.95LCL 0.95UCL
## [1,] 137    128     80      52     105
cat("La survie médiane est de 80 jours ou de",80/365.25, "ans")
## La survie médiane est de 80 jours ou de 0.2190281 ans

n = 137 individus events = 128 personnes ont présentés l’évènement d’intérêt

Courbe de Kaplan-Meier

Avec la fonction plot basique de R

plot(fit, col = "blue")
grid()

Avec les fonctions du package survminer (ggsurvplot par exemple)

Installons d’abord survminer.

#en fonction de l'ordinateur
#install.packages("survminer")
library("survminer")
## Loading required package: ggplot2
## Loading required package: ggpubr
## 
## Attaching package: 'survminer'
## The following object is masked from 'package:survival':
## 
##     myeloma

Courbe de survie basique

fit1 <- survfit(Surv(time, status) ~ 1, data = veteran)
ggsurvplot(fit = fit1, data = veteran)

Sur la même courbe nous rajoutons des éléments plus informatifs

ggsurvplot(fit1, data = veteran,
 surv.median.line = "hv", # Add medians survival
 conf.int = TRUE,
 # Add risk table
 risk.table = TRUE,
 tables.height = 0.2,
 tables.theme = theme_cleantable(),
 ggtheme = theme_bw() # Change ggplot2 theme
)

Lorsque l’évènement d’intérêt est le dècès, la courbe de survie est appelé courbe de survie globale.

Lorsque l’évènement d’intérêt est le fait de ne pas avoir une récidive (rechute), la courbe de survie est appelé courbe de survie sans récidive ou disease free survival.

Analyse de survie bivariée

Estimateur de Kaplan-Meier

Es ce qu’il y a une difference de survie entre les deux groupes de traitement. On utilisera un test de Log-rank avec la fonction survdiff de survival.

** Réponse**

Il s’agir de comparer des distributions de survies.

H0 : il n’y a pas de differences de survies entre les deux groupes.

H1 : il y a une differences de survies entre les deux groupes.

Prédiction

fit2 <- survfit(Surv(time, status) ~ traitement,
                data = veteran)

print(fit2)
## Call: survfit(formula = Surv(time, status) ~ traitement, data = veteran)
## 
##                      n events median 0.95LCL 0.95UCL
## traitement=standard 69     64  103.0      59     132
## traitement=test     68     64   52.5      44      95

Test de logrank comparant des distributions de survie dans les deux groupes de traitement :

(fit3 <- survdiff(Surv(time, status) ~ traitement,
                  data = veteran))
## Call:
## survdiff(formula = Surv(time, status) ~ traitement, data = veteran)
## 
##                      N Observed Expected (O-E)^2/E (O-E)^2/V
## traitement=standard 69       64     64.5   0.00388   0.00823
## traitement=test     68       64     63.5   0.00394   0.00823
## 
##  Chisq= 0  on 1 degrees of freedom, p= 0.9

Conclusion

  • \(P<0.05\)

  • Test non significatif

  • Non rejet de H0 au risque \(\beta = 20\%\)

  • Il n’ya pas de differences de survie entre les deux traitements

ggsurvplot(fit = fit2,
           data = veteran,surv.median.line = "hv", # Add medians survival

 # Change legends: title & labels
 legend.title = "Traitement",
 legend.labs = c("standard", "test"),
 # Add p-value and tervals
 pval = TRUE,
           risk.table = TRUE,
 tables.height = 0.2,
 tables.theme = theme_cleantable(),
 ggtheme = theme_bw() # Change ggplot2 theme
)

Pour aller loin sur cette analyse

https://rpubs.com/xvalda/survival