Dans le cadre de ce TD, nous allons utiliser les données veteran
disponible 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
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"))
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.
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
plot(fit, col = "blue")
grid()
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.
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
)