Chargement des packages utilisés

library(knitr) #pour le rapport Rmarkdown
library(survival)

1. Apercu du jeu de données

Nous recuperons le jeu de données aml du package Survival ensuite nous l’attachons à l’aide de la fonction attach . il est imporant de noter que La fonction attach permet de ne pas à avoir à préciser à chaque fois le jeu de donnée sur lequel on travaille pour acceder a une de ses colonnes

data=aml
attach(data)
kable(aml)

Il est important de faire un peu de statistque descriptive sur les données

kable(summary(data))

n=nrow(data)
nb_cens=sum(status==0) #nombre d'observations censurées
nb_noncens=n-nb_cens#nombre d'observations non censurées

Nb_main=sum(x=='Maintained') #maintained
Nb_nonmain=sum(x=='Nonmaintained')# nonmaintained

L’etude porte dur 23 patients. Nous remarquons que le maximum des temps de survie est de 161 et le minimum est de 5, nous avons aussi 5observations censurées et 18 observations non censurées.De plus il ya 11 patient dont le traitement par chimiothérapie est maintenu , pour les 12 , le traitement n’est pas maintenu.

2. Création d’un objet de type survie

Nous allons crées un objet de type survie: ce objet nous permettra d’utiliser les fonctions du package survival. Pour cela , on utilise la fonction Surv.

base=Surv(time,status)
base #données censurées marquées par un +

##  [1]   9   13   13+  18   23   28+  31   34   45+  48  161+   5    5    8    8 
## [16]  12   16+  23   27   30   33   43   45

Les données censurés sont marqués par un +. on a comme remarqué dans la question précédente, 5 observations censurées.

3. Estimateur de Kapla-Meier

Nous allons calculé l’estimateur de kapla-Meir de la survie globale c’est a dire utilisant le jeu de données en entie, puis calculer l’estimateur pour les 2 groupes formés à partir de la variable x (Maintained et Nonmaintained).

rappelons la formule de l’estimateur de kapla-Meier: \(\widehat{S}(t)=\prod_j^{ t^{*}_j < t }\left(1-\frac{d_j}{Y^{*}_j}\right)\)

avec \(t^{*}_j\) les temps de survie de l’évenement étudié, d_j le nombre de mort en \(t^{*}_j\) et \(Y^{*}_j\) le nombre d’individus à risque.

Pour calculer cet estimateur, on utilise la fonction survit du package survival R.

kpm=survfit(base~1,)
summary(kpm,main =' Courbe de survie globale')

## Call: survfit(formula = base ~ 1)
## 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     5     23       2   0.9130  0.0588       0.8049        1.000
##     8     21       2   0.8261  0.0790       0.6848        0.996
##     9     19       1   0.7826  0.0860       0.6310        0.971
##    12     18       1   0.7391  0.0916       0.5798        0.942
##    13     17       1   0.6957  0.0959       0.5309        0.912
##    18     14       1   0.6460  0.1011       0.4753        0.878
##    23     13       2   0.5466  0.1073       0.3721        0.803
##    27     11       1   0.4969  0.1084       0.3240        0.762
##    30      9       1   0.4417  0.1095       0.2717        0.718
##    31      8       1   0.3865  0.1089       0.2225        0.671
##    33      7       1   0.3313  0.1064       0.1765        0.622
##    34      6       1   0.2761  0.1020       0.1338        0.569
##    43      5       1   0.2208  0.0954       0.0947        0.515
##    45      4       1   0.1656  0.0860       0.0598        0.458
##    48      2       1   0.0828  0.0727       0.0148        0.462

plot(kpm)

kpm1=survfit(base~x,data=data,subset=x=="Maintained")
summary(kpm1,main =' Courbe de survie (maintenue)')

## Call: survfit(formula = base ~ x, data = data, subset = x == "Maintained")
## 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     9     11       1    0.909  0.0867       0.7541        1.000
##    13     10       1    0.818  0.1163       0.6192        1.000
##    18      8       1    0.716  0.1397       0.4884        1.000
##    23      7       1    0.614  0.1526       0.3769        0.999
##    31      5       1    0.491  0.1642       0.2549        0.946
##    34      4       1    0.368  0.1627       0.1549        0.875
##    48      2       1    0.184  0.1535       0.0359        0.944

plot(kpm1)

kpm2=survfit(base~x,data=data,subset=x=="Nonmaintained")
summary(kpm2)

## Call: survfit(formula = base ~ x, data = data, subset = x == "Nonmaintained")
## 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     5     12       2   0.8333  0.1076       0.6470        1.000
##     8     10       2   0.6667  0.1361       0.4468        0.995
##    12      8       1   0.5833  0.1423       0.3616        0.941
##    23      6       1   0.4861  0.1481       0.2675        0.883
##    27      5       1   0.3889  0.1470       0.1854        0.816
##    30      4       1   0.2917  0.1387       0.1148        0.741
##    33      3       1   0.1944  0.1219       0.0569        0.664
##    43      2       1   0.0972  0.0919       0.0153        0.620
##    45      1       1   0.0000     NaN           NA           NA

plot(kpm2,main =' Courbe de survie (non maintenue)')

4.Courbes de survie des deux groupes sur un meme grpahique

kpm3=survfit(Surv(time,status)~ x,data=data)
plot(kpm3,col=c('red','black'),mark=4,cex=2, main =' Courbes de survie des deux groupes')

legend("topright", legend=c("Maintained", "Nonmaintained"),
       col=c("red", "blue"), lty=c(1,1))

5. Moyenne, médiane, quartiles de survie

Nous allons calculé la moyenne, médiane, quartiles de survie dans les deux groupes et dans la pouplation globale

# moyenne, mediane et quantile de survie ( 2 groupe)
model<- survfit(Surv(time,status) ~ x)
a=quantile(model)
#moyenne groupe 1
moy_main=mean(mean(data[which(x=='Maintained'),1]))
#moyenne groupe 2
moy_nonmain=mean(mean(data[which(x=='Nonmaintained'),1]))

#moyenne, mediane et quantile de survie ( global)
model1<- survfit(Surv(time,status) ~ 1, data = data)
#moyenne de survi
moy_glob=mean(model1$time)
#quantile
b=quantile(model1)

#tableau recapitulatif
tab=data.frame(moyenne=c(moy_main,moy_nonmain,moy_glob),mediane=c(a$quantile[,2],b$quantile[2]),fisrt_quant=c(a$quantile[,1],b$quantile[1]),third_quant=c(a$quantile[,3],b$quantile[3]))
rownames(tab)=c('maintenue','non-maintenue','globale')

kable(tab)

Comparaison des parametres

6. Comparaison des deux courbes de survie

Si l’on compare les deux courbes de survie, constate que les patients du groupe dont le traitement de chimiothérapie a été maintenu présentent une probabilité de survie plus grande que ceux du groupe dont le traitement de chimiothérapie n’a pas été maintenu car la courbe des patients “Maintained” se situe toujours en-dessous de celle des patients “Nonmaintained”. ce graphique recèle quelques hypothèses potentiellement difficiles à prouver. par exemple on est pas sur que la date de survenu de l’evenement soit effectivment la bonne . En effetSi l’événement étudié etait par exemple «la mort», il n’y aurait pas de doute mais dans notre cas l’événement considéré est «l’apparition d’une récidive» dans une étude sur le cancer (leucémie mylione), si on suppose la recidive a été observé lors d’un examen de routine, alors on neconnaît pas le moment précis de la récidive, mais seulement l’intervalle de temps pendant lequel elle a eu lieu. nous nous trouvons donc face à un biais, notamment une surestimation du délai effectif de survenue de la récidive et par conséquent de la probabilité de survie sans récidive.

Dans la question suivante on utilisera un test statisque approprié ( test de log- Rank) pour comparer les deux courbes de survie.

7. Mise en place un test de comparaison des deux courbes de survie.

Test de log-Rank

Dans les fonctions de survie, le test convenant le mieux est le test du log-rank, qui fonctionne comme un test d’adéquation et compare le nombre de décès observés et le nombre de décès attendus dans chaque tranche temporelle, sous l’hypothèse nulle qu’il n’y ait pas de différence entre les deux courbes. De plus c’est le plus performant des tests quand les deux courbes de survies ne se croise pas ce qui est le cas dans cette étude.

( Ajouter la statistique de test , sa loi sous ho ,…)

Pour implementer le test de log-rank sur R , on utilise la fonction _survdiff du package survival.

survdiff(Surv(time,status)~x)

## Call:
## survdiff(formula = Surv(time, status) ~ x)
## 
##                  N Observed Expected (O-E)^2/E (O-E)^2/V
## x=Maintained    11        7    10.69      1.27       3.4
## x=Nonmaintained 12       11     7.31      1.86       3.4
## 
##  Chisq= 3.4  on 1 degrees of freedom, p= 0.07

La valeur p du test de log-rank est égale à 0,07 qui est superieur à 0,05, ce qui signifie que l’on a prouvé qu’il existe une différence statistiquement significative entre les probabilités de survie des deux groupes.

8. Deux estimateurs du risque cumulé pour la population globale.