Fonction de survie globale

kpm=survfit(base~1,)
summary(kpm)

## 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,main ='Courbe de survie globale',xlab='time',ylab="SKM")

Fonction de survie groupe 1 (chimiothérapie maintenue)

kpm1=survfit(base~x,data=data,subset=x=="Maintained")
summary(kpm1)

## 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,,main =' Courbe de survie (maintenue)')

Fonction de survie groupe 2 (chimiothérapie non maintenue)

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)')

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.

la statistique de test est \(T_n=\frac{U^{2}}{V}\)

avec $U= {k=0}^{n}w_i (d{B,i}-e_{B,i}) $ et $V=_{k=0}^{n} $

*wi une pondération dependant du choix du test

• \(d_{A,i}\) , \(d_{B,i}\) le nombre d’évenements observés au temps \(t_i^{*}\) dans les groupes A et B. ( \(d_i=d_{A,i} +d_{B,i}\))

• \(e_{A,i}\),\(e_{B,i}\) le nombre d’évenements attendus au temps \(t_i^{*}\) dans les groupes A et B sous H0

• \(R_{A,i}\) , \(R_{B,i}\) le nombre d’individus a risque juste avant \(t_i^{*}\) dans les groupes A et B avec \(R_i=R_{A,i}+R_{B,i}\)

sous (H0), \(T_n \hookrightarrow \chi^2 (1)\)

et la region critique du test est telle que \(R_{\alpha} =\{T_n ≥ c\alpha\}\) où \(c\alpha\) est le quantile d’ordre \(1 − \alpha\) de la loi de \(\chi^2 (1)\).

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 n’ existe une différence statistiquement significative entre les probabilités de survie des deux groupes.Mais ce test peut etre ilusoire pour tirer des conclusions car nous avons un faible effectif de chaque groupe.

Estimateur de Nelson-Alen

Nelson-Alen estime le risque cumulé par une fonction en escalier:

\(\widehat{S}(t)= \sum_{t^{*}_j < t}\frac{d_j}{Y^{*}_j}\)

library(survival)
fit=survfit(base~1)
a=summary(fit)
p=cumsum(a$n.event/a$n.risk)
plot(a$time,p, type="s", main= "Estimateur de Nelson-Alen du risque cumulé") 

Estimateur de Breslow

on sait que la formule (1) peut se réecrire de la maniere suivante :

\(H(u)=-ln(S(u))\)

Breslow estime donc la fonction de risque cumulé par :

\(H(u)=-ln(\widehat{S}(u))\)

fit=summary(survfit(Surv(time,status)~1,type="fh"))
fit

## Call: survfit(formula = Surv(time, status) ~ 1, type = "fh")
## 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     5     23       2    0.915  0.0575       0.8088        1.000
##     8     21       2    0.830  0.0775       0.6910        0.997
##     9     19       1    0.787  0.0844       0.6381        0.971
##    12     18       1    0.745  0.0899       0.5878        0.944
##    13     17       1    0.702  0.0943       0.5397        0.914
##    18     14       1    0.654  0.0995       0.4852        0.881
##    23     13       2    0.557  0.1057       0.3840        0.808
##    27     11       1    0.509  0.1070       0.3367        0.768
##    30      9       1    0.455  0.1083       0.2855        0.725
##    31      8       1    0.402  0.1079       0.2372        0.680
##    33      7       1    0.348  0.1060       0.1917        0.632
##    34      6       1    0.295  0.1023       0.1493        0.582
##    43      5       1    0.241  0.0966       0.1101        0.529
##    45      4       1    0.188  0.0887       0.0745        0.474
##    48      2       1    0.114  0.0784       0.0296        0.439

plot(fit$time,-log(fit$surv),type="s",col="red",xlab='time',,main="Estimateur de Breslow du risque cumulé")

On remarque que les estimateurs de Nelson-Aalen et Breslow sont quasiment égaux.

Determination de nk,Rk,fk,Sk,pk

• Détermination de \(n_k\)

\(n_k\):nombre de cas exposés au risque après une durée correspondant à la borne inférieure de la classe k.

\(N_{k}=N_{k-1}-nouvelemploi_{k-1}-Autresraisons_{k-1}\)

nk=rep(0,length(l)-1)
nk[1]=202
for (i in 2:length(nk)){
  nk[i]=nk[i-1]-Nouvelemploi[i-1]-Autresraisons[i-1]
}

• Détermination de \(R_{k}\)

\(R_{k}\):nombre de cas exposés utilisé pour l’estimation du risque. \(R_{k}=n_{k}-\frac{Autres raisons_{k}}{2}\)

Rk=rep(0,length(l)-1)
for(i in 1:length(Rk)){
  Rk[i]=nk[i]-Autresraisons[i]/2
}

• Détermination de \(S_{k}\)

\(S_{k}\):estimation de la probabilité de rester au chômage pour une durée égale ou supérieure à celle de la classe k. \(S_{k}= \prod_{a_{i} \geq a_{j}}(1-\frac{nouvelemploi_{i-1}}{R_{i-1}})\)

lk=rep(0,length(l)-1)
lk[1]=1
for (i in 2:length(lk)){
  lk[i]=1-Nouvelemploi[i-1]/Rk[i-1]
}
Sk=rep(0,length(lk))
Sk=lk[1]
for(i in 2:length(lk)){
  for(j in 2:i){
    Sk[i]=Sk[j-1]*lk[i]
  }
}

• Détermination de \(p_{k}\)

\(p_{k}\):estimation de la probabilité de rester au chômage au moins jusqu’au début et au plus tard jusqu’à la fin de la classe k. \(p_{k}= \frac{S_{k}*\frac{nouvelemploi_{k}}{R_{k}}}{amplitude_{k}}\)

pk=rep(0,length(l))
amplitude=c(2,2,2,2,7)
for (i in 1:length(pk)){
  pk[i]=Sk[i]*(Nouvelemploi[i]/Rk[i])/amplitude[i]
}
pk=pk[-5]

• Détermination de \(h_{k}\)

\(h_{k}\) :estimation du risque instantané par intervalle. \(h_{k}=\frac{\frac{nouvelemploi_{k}}{R_{k}}}{amplitude_{k}(1-\frac{\frac{nouvelemploi_{k}}{R_{k}}}{2})}\)

hk=rep(0,length(pk))
for(i in 1:length(hk)){
  hk[i]=(pk[i]/(1-(Nouvelemploi[i]/Rk[i])/2))/Sk[i]
}

Tableau complèté

donnee=cbind(chomage,nk,Rk,hk,Sk,pk)
kable(donnee)

Chargement des données

#importation
load(file="C:/Users/User/Desktop/Survival Time/Survival Time/chomage_tp.rda")

Construction des variables neccesaires pour utiliser lifetab

tableau=Dataset
Tempsenmois=cut(tableau$temps,l)
tableau1=table(Tempsenmois)
tableau1=as.data.frame(tableau1)
nouveauemploi=rep(0,nrow(tableau1)-1)
for(i in 1:(length(l)-1)){
  nouveauemploi[i]=nrow(tableau[(tableau$temps<l[i+1] & tableau$temps>=l[i]) & tableau$statut=="evenement",])
}

autresraisons=tableau1$Freq-nouveauemploi
tableau1=cbind(tableau1,nouveauemploi,autresraisons)
kable(tableau1)

Construction de la table de survie avec lifetab

n=nrow(tableau)
surv=lifetab(l,n,autresraisons,nouveauemploi)
kable(surv)

Representation de la fonction de survie**

x=rep(l,each=2)
x=x[-c(1,12)]
y=rep(surv$surv,each=2)
plot(x,y,type="l",col="red",xlab="intervalle de temps",ylab="fonction survie")
title(main = "Courbe de survie globale", cex=.6)

Comparaison des résultats avec ceux obtenus précédemment et commentaire.

les résultats sont presque identique , il existe une légère différence au niveau l’intervalle [8,15]. cette différence est due au faite que les informations sur l’intervalle [8,15] du jeu du données et du tableau ne concorde pas.

Table de survie pour les hommes

tableau=Dataset
tableau=tableau[tableau$genre=="H",]
Tempsenmois=cut(tableau$temps,l)
tableau1=table(Tempsenmois)
tableau1=as.data.frame(tableau1)
nouveauemploi=rep(0,nrow(tableau1)-1)
for(i in 1:(length(l)-1)){
  nouveauemploi[i]=nrow(tableau[(tableau$temps<l[i+1] & tableau$temps>=l[i]) & tableau$statut=="evenement",])
}

autresraisons=tableau1$Freq-nouveauemploi
tableau1=cbind(tableau1,nouveauemploi,autresraisons)
kable(tableau1)

n=nrow(tableau)
survH=lifetab(l,n,autresraisons,nouveauemploi)
kable(survH)

Table de survie pour les femmes

tableau=Dataset
tableau=tableau[tableau$genre=="F",]
Tempsenmois=cut(tableau$temps,l)
tableau1=table(Tempsenmois)
tableau1=as.data.frame(tableau1)
nouveauemploi=rep(0,nrow(tableau1)-1)
for(i in 1:(length(l)-1)){
  nouveauemploi[i]=nrow(tableau[(tableau$temps<l[i+1] & tableau$temps>=l[i]) & tableau$statut=="evenement",])
}

autresraisons=tableau1$Freq-nouveauemploi
tableau1=cbind(tableau1,nouveauemploi,autresraisons)
kable(tableau1)

n=nrow(tableau)
survF=lifetab(l,n,autresraisons,nouveauemploi)
kable(survF)

construction des deux courbes de survie

x=rep(l,each=2)
x=x[-c(1,12)]
yH=rep(survH$surv,each=2)
plot(x,yH,type="l",col="red",xlab="intervalle de temps",ylab="fonction survie",ylim=seq(0,1,1))

x=rep(l,each=2)
x=x[-c(1,12)]
yF=rep(survF$surv,each=2)
lines(x,yF,type="l",col="blue",xlab="intervalle de temps",ylab="fonction survie")
title(main = "Courbe de survie (Homme vs femme)", cex=.6)


legend("topright", legend=c("Homme ", "Femmme"),
       col=c("red", "blue"), lty=c(1,1))

la courbe de survie des hommes est au dessus de la courbe de survie des femmes ce qui veut dire que les hommes retrouvent plus rapidement du travail que les femmes.