# Exemples du chapitre 3. Courbes de survie ################################ # Temps discret: Breast cancer # ################################ # Méthode actuarielle library(KMsurv) options(digits=3) # vecteur représentant les années tis <- c(0,1,2,3,4,5,6,7,8,9,10,11,12) # variable représentant le nombre initial de sujets à risque ninit <- 1207 # vecteur représentant le nombre de sujets censurés chaque année nlost <- c(129,183,147,166,153,121,91,59,39,25,19,3) # vecteur représentant le nombre de sujet subissant l'événement chaque année nevent <- c(2,15,14,20,8,5,7,0,1,0,0,0) # Table de survie survie_bc1 <- lifetab(tis,ninit,nlost,nevent) survie_bc1 # Rem: Les probabilités de survie sont calculées à partir des probabilités # de chaque année et non du hasard instantané. # Plot de la fonction de survie zy <- kronecker(survie_bc1$surv,c(1,1)) zx <- c(tis[1],kronecker(tis[2:12],c(1,1)),tis[13]) win.graph(width=8, height=6) plot(zx,zy,cex=1.3,cex.axis=1.3,cex.lab=1.3,col="blue2",lwd=2,type="l",ylim=c(0.85,1), main="Fonction de survie", xlab="Temps en années", ylab="Survie cumulée") ########################### # Même chose, mais par mois # vecteur représentant les mois tis12 <- c(0,12,24,36,48,60,72,84,96,108,120,132,144) # Table de survie survie_bc2 <- lifetab(tis12,ninit,nlost,nevent) survie_bc2 ########################################################### # Plot de la fonction de survie avec calcul du quantile 90% zy <- kronecker(survie_bc1$surv,c(1,1)) zx <- c(tis[1],kronecker(tis[2:12],c(1,1)),tis[13]) win.graph(width=8, height=6) plot(zx,zy,cex=1.3,cex.axis=1.3,cex.lab=1.3,col="blue2",lwd=2,type="l",ylim=c(0.85,1), main="Fonction de survie", xlab="Temps en années", ylab="Survie cumulée") lines(c(-1,7,7),c(0.9,0.9,0.8),lwd=2,col="red2",lty=2) ############################ # Temps discret: Yamaguchi # ############################ # vecteur représentant les intervalles en mois tis <- c(0,12,24,36,48) # variable représentant le nombre initial de sujets à risque ninit <- 265 # vecteur représentant le nombre de sujets censurés chaque année nlost <- c(0,0,3,155) # vecteur représentant le nombre de sujet subissant l'événement chaque année nevent <- c(59,23,18,7) # Table de survie globale survie_yam1 <- lifetab(tis,ninit,nlost,nevent) survie_yam1 # Plot de la fonction de survie zy <- kronecker(survie_yam1$surv,c(1,1)) zx <- c(tis[1],kronecker(tis[2:4],c(1,1)),tis[5]) win.graph(width=8, height=6) plot(zx,zy,cex=1.3,cex.axis=1.3,cex.lab=1.3,col="blue2",lwd=2,type="l", main="Fonction de survie", xlab="Temps en mois", ylab="Survie cumulée") ######################################### # Tables de survie séparées hommes/femmes # Hommes # vecteur représentant les intervalles en mois tis <- c(0,12,24,36,48) # variable représentant le nombre initial de sujets à risque ninit <- 114 # vecteur représentant le nombre de sujets censurés chaque année nlost <- c(0,0,1,71) # vecteur représentant le nombre de sujet subissant l'événement chaque année nevent <- c(22,10,8,2) # Table de survie survie_yam2H <- lifetab(tis,ninit,nlost,nevent) survie_yam2H # Femmes # vecteur représentant les intervalles en mois tis <- c(0,12,24,36,48) # variable représentant le nombre initial de sujets à risque ninit <- 151 # vecteur représentant le nombre de sujets censurés chaque année nlost <- c(0,0,2,84) # vecteur représentant le nombre de sujet subissant l'événement chaque année nevent <- c(37,13,10,5) # Table de survie survie_yam2F <- lifetab(tis,ninit,nlost,nevent) survie_yam2F # Plot des deux fonctions de survie zyH <- kronecker(survie_yam2H$surv,c(1,1)) zyF <- kronecker(survie_yam2F$surv,c(1,1)) zx <- c(tis[1],kronecker(tis[2:4],c(1,1)),tis[5]) win.graph(width=8, height=6) plot(zx,zyH,cex=1.3,cex.axis=1.3,cex.lab=1.3,col="blue2",lwd=2,type="l",ylim=c(0.5,1), main="Fonctions de survie", xlab="Temps en mois", ylab="Survie cumulée") lines(zx,zyF,col="red2",lwd=2) legend(35,0.95, c("Hommes", "Femmes"), col=c("blue2","red2"),lwd=2, lty = 1) ################################################## # Intervalle de confiance pour une courbe de survie # Yamaguchi # vecteur représentant les intervalles en mois tis <- c(0,12,24,36,48) # Avec les erreurs standards données par R yam_ci_inf_1 <- survie_yam1$surv - 1.96*survie_yam1$se.surv yam_ci_sup_1 <- survie_yam1$surv + 1.96*survie_yam1$se.surv yam_ci_inf_1 yam_ci_sup_1 # Avec la formule donnée par Altman Altman_se <- sqrt((1-survie_yam1$surv)/survie_yam1$nsubs) yam_ci_inf_2 <- survie_yam1$surv - 1.96*Altman_se yam_ci_sup_2 <- survie_yam1$surv + 1.96*Altman_se yam_ci_inf_2 yam_ci_sup_2 # Plot de la courbe de survie avec l'intervalle de confiance (Altman) # Plot de la fonction de survie zy <- kronecker(survie_yam1$surv,c(1,1)) zx <- c(tis[1],kronecker(tis[2:4],c(1,1)),tis[5]) win.graph(width=8, height=6) plot(zx,zy,cex=1.3,cex.axis=1.3,cex.lab=1.3,col="blue2",lwd=2,type="l",ylim=c(0.5,1), main="Fonction de survie", xlab="Temps en mois", ylab="Survie cumulée") # CI 95% zy_inf <- kronecker(yam_ci_inf_2,c(1,1)) zy_sup <- kronecker(yam_ci_sup_2,c(1,1)) lines(zx,zy_inf,col="red2",lwd=2,lty=2) lines(zx,zy_sup,col="red2",lwd=2,lty=2) legend(35,0.95, c("Survie", "95% CI (Altman)"), col=c("blue2","red2"),lwd=2, lty = c(1,2)) ####################### # Temps discret: TREE # ####################### load("TREE_4var.rda") # Courbe de survie globale avec CI # vecteur représentant les intervalles en années tis <- c(0:7) # variable représentant le nombre initial de sujets à risque ninit <- nrow(Dataset) # vecteur représentant le nombre de sujets censurés chaque année nlost <- rep(0,7) for (i in 1:7) { nlost[i] <- length(which(Dataset$evenement==0 & Dataset$duree==i)) } # vecteur représentant le nombre de sujet subissant l'événement chaque année nevent <- rep(0,7) for (i in 1:7) { nevent[i] <- length(which(Dataset$evenement==1 & Dataset$duree==i)) } # Table de survie survie_tree1 <- lifetab(tis,ninit,nlost,nevent) survie_tree1 # CI # Avec les erreurs standards données par R tree_ci_inf_1 <- survie_tree1$surv - 1.96*survie_tree1$se.surv tree_ci_sup_1 <- survie_tree1$surv + 1.96*survie_tree1$se.surv tree_ci_inf_1 tree_ci_sup_1 # Plot de la courbe de survie avec l'intervalle de confiance # Plot de la fonction de survie zy <- kronecker(survie_tree1$surv,c(1,1)) zx <- c(tis[1],kronecker(tis[2:7],c(1,1)),tis[8]) win.graph(width=8, height=6) plot(zx,zy,cex=1.3,cex.axis=1.3,cex.lab=1.3,col="blue2",lwd=2,type="l",ylim=c(0.8,1), main="Fonction de survie", xlab="Temps en années", ylab="Survie cumulée") # CI 95% zy_inf <- kronecker(tree_ci_inf_1,c(1,1)) zy_sup <- kronecker(tree_ci_sup_1,c(1,1)) lines(zx,zy_inf,col="red2",lwd=2,lty=2) lines(zx,zy_sup,col="red2",lwd=2,lty=2) legend(35,0.95, c("Survie", "95% CI"), col=c("blue2","red2"),lwd=2, lty = c(1,2)) # Courbes de survie par sexe avec CI # Hommes # vecteur représentant les intervalles en années tis <- c(0:7) # variable représentant le nombre initial de sujets à risque ninit <- length(which(Dataset$genre=="male")) # vecteur représentant le nombre de sujets censurés chaque année nlost <- rep(0,7) for (i in 1:7) { nlost[i] <- length(which(Dataset$evenement==0 & Dataset$duree==i & Dataset$genre=="male")) } # vecteur représentant le nombre de sujet subissant l'événement chaque année nevent <- rep(0,7) for (i in 1:7) { nevent[i] <- length(which(Dataset$evenement==1 & Dataset$duree==i & Dataset$genre=="male")) } # Table de survie survie_tree2H <- lifetab(tis,ninit,nlost,nevent) survie_tree2H # CI # Avec les erreurs standards données par R tree_ci_inf_2H <- survie_tree2H$surv - 1.96*survie_tree2H$se.surv tree_ci_sup_2H <- survie_tree2H$surv + 1.96*survie_tree2H$se.surv tree_ci_inf_2H tree_ci_sup_2H # Femmes # vecteur représentant les intervalles en années tis <- c(0:7) # variable représentant le nombre initial de sujets à risque ninit <- length(which(Dataset$genre=="female")) # vecteur représentant le nombre de sujets censurés chaque année nlost <- rep(0,7) for (i in 1:7) { nlost[i] <- length(which(Dataset$evenement==0 & Dataset$duree==i & Dataset$genre=="female")) } # vecteur représentant le nombre de sujet subissant l'événement chaque année nevent <- rep(0,7) for (i in 1:7) { nevent[i] <- length(which(Dataset$evenement==1 & Dataset$duree==i & Dataset$genre=="female")) } # Table de survie survie_tree2F <- lifetab(tis,ninit,nlost,nevent) survie_tree2F # CI # Avec les erreurs standards données par R tree_ci_inf_2F <- survie_tree2F$surv - 1.96*survie_tree2F$se.surv tree_ci_sup_2F <- survie_tree2F$surv + 1.96*survie_tree2F$se.surv tree_ci_inf_2F tree_ci_sup_2F # Plot des deux courbes de survie avec les intervalles de confiance # Plot des fonctions de survie # Hommes zy <- kronecker(survie_tree2H$surv,c(1,1)) zx <- c(tis[1],kronecker(tis[2:7],c(1,1)),tis[8]) win.graph(width=8, height=6) plot(zx,zy,cex=1.3,cex.axis=1.3,cex.lab=1.3,col="blue2",lwd=2,type="l",ylim=c(0.7,1), main="Fonctions de survie", xlab="Temps en années", ylab="Survie cumulée") # Femmes zy <- kronecker(survie_tree2F$surv,c(1,1)) lines(zx,zy,col="red2",lwd=2) # CI 95% # Hommes zy_inf <- kronecker(tree_ci_inf_2H,c(1,1)) zy_sup <- kronecker(tree_ci_sup_2H,c(1,1)) lines(zx,zy_inf,col="blue2",lwd=2,lty=2) lines(zx,zy_sup,col="blue2",lwd=2,lty=2) # Femmes zy_inf <- kronecker(tree_ci_inf_2F,c(1,1)) zy_sup <- kronecker(tree_ci_sup_2F,c(1,1)) lines(zx,zy_inf,col="red2",lwd=2,lty=2) lines(zx,zy_sup,col="red2",lwd=2,lty=2) legend(0,0.75, c("Survie hommes (95% CI)", "Survie femmes (95% CI)"), col=c("blue2","red2"),lwd=2) ############################## # Temps discret: 2 immeubles # ############################## load("2immeubles.rda") # Courbes de survie par immeuble # Immeuble 1 # vecteur représentant les intervalles en années tis <- c(0:10) # variable représentant le nombre initial de sujets à risque ninit <- Dataset$n[1] # vecteur représentant le nombre de sujets censurés chaque année nlost <- Dataset$w[1:10] # vecteur représentant le nombre de sujet subissant l'événement chaque année nevent <- Dataset$d[1:10] # Table de survie survie_immeuble1 <- lifetab(tis,ninit,nlost,nevent) survie_immeuble1 # Immeuble 2 # vecteur représentant les intervalles en années tis <- c(0:9) # variable représentant le nombre initial de sujets à risque ninit <- Dataset$n[11] # vecteur représentant le nombre de sujets censurés chaque année nlost <- Dataset$w[11:19] # vecteur représentant le nombre de sujet subissant l'événement chaque année nevent <- Dataset$d[11:19] # Table de survie survie_immeuble2 <- lifetab(tis,ninit,nlost,nevent) survie_immeuble2 # Plot des deux courbes de survie # Immeuble 1 tis <- c(0:10) zy <- kronecker(survie_immeuble1$surv,c(1,1)) zx <- c(tis[1],kronecker(tis[2:10],c(1,1)),tis[11]) win.graph(width=8, height=6) plot(zx,zy,cex=1.3,cex.axis=1.3,cex.lab=1.3,col="blue2",lwd=2,type="l",ylim=c(0,1), main="Fonctions de survie", xlab="Temps en années", ylab="Survie cumulée") # Immeuble 2 tis <- c(0:9) zy <- kronecker(survie_immeuble2$surv,c(1,1)) zx <- c(tis[1],kronecker(tis[2:9],c(1,1)),tis[10]) lines(zx,zy,col="red2",lwd=2) legend(0,0.2, c("Immeuble 1", "Immeuble 2"), col=c("blue2","red2"),lwd=2) ############################### # Temps continu: Kaplan-Meier # ############################### library(survival) ################################## # Exemple miniature (10 personnes) load("Exemple_miniature.rda") # Courbe de survie survie <- Surv(Dataset$date,Dataset$statut) survie_KMmin <- survfit(survie~1) survie_KMmin summary(survie_KMmin) # Calcul du temps moyen de survie nr <- length(survie_KMmin$time) mean_stime <- survie_KMmin$time[1] + (survie_KMmin$time[2:nr]-survie_KMmin$time[1:nr-1])%*%survie_KMmin$surv[1:nr-1] mean_stime # Plot de la fonction de survie win.graph(width=8, height=6) plot(survie_KMmin,cex=1.3,cex.axis=1.3,cex.lab=1.3,col="blue2",lwd=2, main="Fonction de survie", xlab="Durée en années avant la survenance de l'événement", ylab="Survie cumulée") legend(1,0.3, c("Fonction de survie", "Censure", "IC 95%"), col="blue2", lwd=c(2,1,2), lty = c(1,-1,2), pch = c(-1,3,-1)) ################### # Exemple Yamaguchi load("Yamaguchi.rda") # Courbe de survie globale survie <- Surv(Dataset$dur,Dataset$evt) survie_KMyam1 <- survfit(survie~1) survie_KMyam1 summary(survie_KMyam1) # Calcul du temps moyen de survie nr <- length(survie_KMyam1$time) mean_stime <- survie_KMyam1$time[1] + (survie_KMyam1$time[2:nr]-survie_KMyam1$time[1:nr-1])%*%survie_KMyam1$surv[1:nr-1] mean_stime # Plot de la fonction de survie win.graph(width=8, height=6) plot(survie_KMyam1,cex=1.3,cex.axis=1.3,cex.lab=1.3,col="blue2",lwd=2, main="Fonction de survie", xlab="Durée en mois avant la survenance de l'événement", ylab="Survie cumulée") legend(5,0.3, c("Fonction de survie", "Censure", "IC 95%"), col="blue2", lwd=c(2,1,2), lty = c(1,-1,2), pch = c(-1,3,-1)) ########################################### # 2 fonctions de survie en fonction du sexe survie_KMyam2 <- survfit(survie~strata(Dataset$sex)) survie_KMyam2 summary(survie_KMyam2) # Plot des 2 fonctions de survie win.graph(width=8, height=6) plot(survie_KMyam2,cex=1.3,cex.axis=1.3,cex.lab=1.3,col=c("blue2","red2"),lwd=2, lty=c(1,2), main="Fonctions de survie", xlab="Durée en mois avant la survenance de l'événement", ylab="Survie cumulée") legend(5,0.3, c("Hommes", "Femmes"),col=c("blue2","red2"),lwd=2, lty = c(1,2)) # Test # Log-rank survdiff(survie~Dataset$sex) # Peto & Peto modification of the Gehan-Wilcoxon test survdiff(survie~Dataset$sex,rho=1) ######################################################################### # Fonctions de survie en fonction de la note moyenne à l'école secondaire survie_KMyam3 <- survfit(survie~strata(Dataset$grd)) survie_KMyam3 summary(survie_KMyam3) # Plot des 5 fonctions de survie win.graph(width=8, height=6) plot(survie_KMyam3,cex=1.3,cex.axis=1.3,cex.lab=1.3,col=c("black","blue2","red2","green2","cyan2"),lwd=2, lty=c(1:5), main="Fonctions de survie", xlab="Durée en mois avant la survenance de l'événement", ylab="Survie cumulée") legend(5,0.3, c("A","A-B","B","B-C","C"),col=c("black","blue2","red2","green2","cyan2"),lwd=2, lty = c(1:5)) # Test # Log-rank survdiff(survie~Dataset$grd) ############################### # Temps continu: Nelson-Aalen # ############################### ################################## # Exemple miniature (10 personnes) load("Exemple_miniature.rda") # Courbe de survie survie <- Surv(Dataset$date,Dataset$statut) survie_NAmin <- survfit(survie~1,type="fleming-harrington") survie_NAmin summary(survie_NAmin) # Calcul du temps moyen de survie nr <- length(survie_NAmin$time) mean_stime <- survie_NAmin$time[1] + (survie_NAmin$time[2:nr]-survie_NAmin$time[1:nr-1])%*%survie_NAmin$surv[1:nr-1] mean_stime # Plot de la fonction de survie win.graph(width=8, height=6) plot(survie_NAmin,cex=1.3,cex.axis=1.3,cex.lab=1.3,col="blue2",lwd=2, main="Fonction de survie", xlab="Durée en années avant la survenance de l'événement", ylab="Survie cumulée") legend(1,0.3, c("Fonction de survie", "Censure", "IC 95%"), col="blue2", lwd=c(2,1,2), lty = c(1,-1,2), pch = c(-1,3,-1))