Étapes pour effectuer une analyse de survie dans R

Envie de partager votre contenu sur les R-bloggers ? cliquez ici si vous avez un blog, ou ici si vous ne le faites pas.
Une autre façon d’analyser ?
Alors qu’il existe tant d’outils et de techniques de modélisation de prédiction, pourquoi avons-nous un autre domaine connu sous le nom d’analyse de survie ? En tant que l’une des branches les plus populaires de la statistique, l’analyse de survie est un moyen de prédiction à différents moments. C’est-à-dire que, alors que d’autres modèles de prédiction prédisent si un événement se produira, l’analyse de survie prédit si l’événement se produira à un moment spécifié. Ainsi, il nécessite une composante temporelle pour la prédiction et, en conséquence, prédit le moment où un événement se produira. Cela aide à comprendre la durée prévue des événements et fournit des informations beaucoup plus utiles. On peut penser aux domaines naturels d’application de l’analyse de survie qui incluent les sciences biologiques où l’on peut prédire le temps nécessaire pour que les bactéries ou autres organismes cellulaires se multiplient à une taille particulière ou le temps prévu de désintégration des atomes. Certaines applications intéressantes incluent la prédiction de l’heure à laquelle une machine tombera en panne et une maintenance sera nécessaire
À quel point cela devient-il difficile..
Il n’est pas facile d’appliquer d’emblée les concepts de l’analyse de survie. Il faut d’abord comprendre les façons dont il peut être utilisé. Cela inclut les courbes de Kaplan-Meier, créant la fonction de survie grâce à des outils tels que les arbres de survie ou les forêts de survie et le test du log-rank. Examinons chacun d’eux un par un dans R. Nous utiliserons le package de survie dans R comme exemple de départ. Le package de survie a la fonction surv() qui est le centre de l’analyse de survie.
# install.packages("survival") # Loading the package library("survival")
Le package contient un exemple de jeu de données à des fins de démonstration. L’ensemble de données est pbc qui contient une étude de 10 ans de 424 patients atteints de cirrhose biliaire primitive (CBP) lorsqu’ils sont traités à la clinique Mayo. Un point à noter ici à partir de la description de l’ensemble de données est que sur 424 patients, 312 ont participé à l’essai du médicament D-pénicillamine et les 112 autres ont consenti à ce que leurs mesures de base soient enregistrées et suivies pour la survie mais n’ont pas participé à l’essai. 6 de ces 112 cas ont été perdus. Nous sommes particulièrement intéressés par les fonctionnalités « temps » et « état » dans l’ensemble de données. Le temps représente le nombre de jours après l’enregistrement et le statut final (qui peut être censuré, transplantation hépatique ou mort). Puisqu’il s’agit de survie, nous considérerons le statut comme mort ou non-mort (greffe ou censuré). De plus amples détails sur l’ensemble de données peuvent être lus à partir de la commande :
#Dataset description ?pbc
Nous commençons par une application directe de la fonction Surv() et la passons à la fonction survfit(). La fonction Surv() prendra les paramètres de temps et d’état et en créera un objet de survie. La fonction survfit() prend un objet de survie (celui que Surv() produit) et crée les courbes de survie.
#Fitting the survival model survival_func=survfit(Surv(pbc$time,pbc$status == 2)~1) survival_func Call: survfit(formula = Surv(pbc$time, pbc$status == 2) ~ 1) n events median 0.95LCL 0.95UCL 418 161 3395 3090 3853
La fonction nous donne le nombre de valeurs, le nombre de positifs en statut, le temps médian et les valeurs de l’intervalle de confiance à 95%. Le modèle peut également être tracé.
#Plot the survival model plot(survival_func)
Comme prévu, le graphique nous montre les probabilités décroissantes de survie au fil du temps. Les lignes en pointillés sont les intervalles de confiance supérieur et inférieur. Dans la fonction survfit() ici, nous avons passé la formule comme ~ 1 qui indique que nous demandons à la fonction d’ajuster le modèle uniquement sur la base de l’objet de survie et donc d’avoir une interception. La sortie ainsi que les intervalles de confiance sont en fait des estimations de Kaplan-Meier. Cette estimation est importante dans l’analyse de survie de la recherche médicale. Les estimations Kaplan – Meier sont basées sur le nombre de patients (chaque patient comme une ligne de données) par rapport au nombre total qui survivent pendant un certain temps après le traitement. (ce qui est l’événement). On peut représenter la fonction de Kaplan – Meier par la formule :
Ŝ Here, di the number of events and ni is the total number of people at risk at time ti
Que faire du graphique ?
Contrairement à d’autres techniques d’apprentissage automatique où l’on utilise des échantillons de test et fait des prédictions sur eux, la courbe d’analyse de survie est une courbe qui s’explique d’elle-même. D’après la courbe, nous voyons que la possibilité de survivre environ 1000 jours après le traitement est d’environ 0,8 ou 80%. Nous pouvons également définir la probabilité de survie pour différents nombres de jours après le traitement. Dans le même temps, nous avons également les plages d’intervalles de confiance qui montrent la marge d’erreur attendue. Par exemple, dans le cas d’un exemple de survie à 1000 jours, l’intervalle de confiance supérieur atteint environ 0,85 ou 85 % et descend à environ 0,75 ou 75 %. Affichez la plage de données, qui est de 10 ans ou environ 3500 jours, les calculs de probabilité sont très erratiques et vagues et ne doivent pas être repris. Par exemple, si l’on veut connaître la probabilité de survivre 4 500 jours après le traitement, alors que le graphique Kaplan – Meier ci-dessus montre une plage comprise entre 0,25 et 0,55 qui est elle-même une grande valeur pour tenir compte du manque de données, les données ne sont toujours pas suffisamment et de meilleures données devraient être utilisées pour faire une telle estimation.
Modèles alternatifs : modèle de risque proportionnel de Cox
Le package de survie contient également une fonction de risque proportionnel cox coxph() et utilise d’autres fonctionnalités dans les données pour créer un meilleur modèle de survie. Bien que les données aient des valeurs manquantes non traitées, je saute le traitement des données et j’ajuste directement le modèle. Dans la pratique, cependant, il faut étudier les données et chercher des moyens de traiter les données de manière appropriée afin que les meilleurs modèles possibles soient ajustés. Comme l’intention de cet article est de familiariser les lecteurs avec la fonction plutôt que de la traiter, l’application de la fonction est l’étape de raccourci que je prends.
# Fit Cox Model Cox_model = coxph(Surv(pbc$time,pbc$status==2) ~.,data=pbc) summary(Cox_model) Call: coxph(formula = Surv(pbc$time, pbc$status == 2) ~ ., data = pbc) n= 276, number of events= 111 (142 observations deleted due to missingness) coef exp(coef) se(coef) z Pr(>|z|) id -2.729e-03 9.973e-01 1.462e-03 -1.866 0.06203 . trt -1.116e-01 8.944e-01 2.156e-01 -0.518 0.60476 age 3.191e-02 1.032e+00 1.200e-02 2.659 0.00784 ** sexf -3.822e-01 6.824e-01 3.074e-01 -1.243 0.21378 ascites 6.321e-02 1.065e+00 3.874e-01 0.163 0.87038 hepato 6.257e-02 1.065e+00 2.521e-01 0.248 0.80397 spiders 7.594e-02 1.079e+00 2.448e-01 0.310 0.75635 edema 8.860e-01 2.425e+00 4.078e-01 2.173 0.02980 * bili 8.038e-02 1.084e+00 2.539e-02 3.166 0.00155 ** chol 5.151e-04 1.001e+00 4.409e-04 1.168 0.24272 albumin -8.511e-01 4.270e-01 3.114e-01 -2.733 0.00627 ** copper 2.612e-03 1.003e+00 1.148e-03 2.275 0.02290 * alk.phos -2.623e-05 1.000e+00 4.206e-05 -0.624 0.53288 ast 4.239e-03 1.004e+00 1.941e-03 2.184 0.02894 * trig -1.228e-03 9.988e-01 1.334e-03 -0.920 0.35741 platelet 7.272e-04 1.001e+00 1.177e-03 0.618 0.53660 protime 1.895e-01 1.209e+00 1.128e-01 1.680 0.09289 . stage 4.468e-01 1.563e+00 1.784e-01 2.504 0.01226 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 exp(coef) exp(-coef) lower .95 upper .95 id 0.9973 1.0027 0.9944 1.000 trt 0.8944 1.1181 0.5862 1.365 age 1.0324 0.9686 1.0084 1.057 sexf 0.6824 1.4655 0.3736 1.246 ascites 1.0653 0.9387 0.4985 2.276 hepato 1.0646 0.9393 0.6495 1.745 spiders 1.0789 0.9269 0.6678 1.743 edema 2.4253 0.4123 1.0907 5.393 bili 1.0837 0.9228 1.0311 1.139 chol 1.0005 0.9995 0.9997 1.001 albumin 0.4270 2.3422 0.2319 0.786 copper 1.0026 0.9974 1.0004 1.005 alk.phos 1.0000 1.0000 0.9999 1.000 ast 1.0042 0.9958 1.0004 1.008 trig 0.9988 1.0012 0.9962 1.001 platelet 1.0007 0.9993 0.9984 1.003 protime 1.2086 0.8274 0.9690 1.508 stage 1.5634 0.6397 1.1020 2.218 Concordance= 0.849 (se = 0.031 ) Rsquare= 0.462 (max possible= 0.981 ) Likelihood ratio test= 171.3 on 18 df, p=0 Wald test = 172.5 on 18 df, p=0 Score (logrank) test = 286.1 on 18 df, p=0
La sortie du modèle de Cox est similaire à la façon dont une sortie de régression linéaire apparaît. Le R2 n’est que de 46%, ce qui n’est pas élevé et nous n’avons aucune caractéristique très significative. Les caractéristiques les plus importantes semblent être l’âge, la bilirubine (bili) et l’albumine. Voyons à quoi ressemble l’intrigue.
#Create a survival curve from the cox model Cox_curve
With more data, we get a different plot and this one is more volatile. Compared to the Kaplan – Meier curve, the cox-plot curve is higher for the initial values and lower for the higher values. The major reason for this difference is the inclusion of variables in cox-model. The plots are made by similar functions and can be interpreted the same way as the Kaplan – Meier curve.Going traditional : Using survival forests
Random forests can also be used for survival analysis and the ranger package in R provides the functionality. However, the ranger function cannot handle the missing values so I will use a smaller data with all rows having NA values dropped. This will reduce my data to only 276 observations.#Using the Ranger package for survival analysis Install.packages("ranger") library(ranger) #Drop rows with NA values pbc_nadrop=pbc[complete.cases(pbc), ] #Fitting the random forest ranger_model
Let’s look at the variable importance plot which the random forest model calculates.#Get the variable importance data.frame(sort(ranger_model$variable.importance,decreasing = TRUE)) sort.ranger_model.variable.importance..decreasing...TRUE. bili 0.0762338981 copper 0.0202733989 albumin 0.0165070226 age 0.0130134413 edema 0.0122113704 ascites 0.0115315711 chol 0.0092889960 protime 0.0060215073 id 0.0055867915 ast 0.0049932803 stage 0.0030225398 hepato 0.0029290675 trig 0.0028869184 platelet 0.0012958105 sex 0.0010639806 spiders 0.0005210531 alk.phos 0.0003291581 trt -0.0002020952Ces nombres peuvent être différents pour différentes courses. Dans mon exemple, nous voyons que la bilirubine est la caractéristique la plus importante.
Leçons apprises : Conclusion
Bien que les données d'entrée pour l'estimation Kaplan - Meier du package Survival, le modèle Cox et le modèle ranger soient toutes différentes, nous comparerons les méthodologies en les traçant sur le même graphique à l'aide de ggplot.
#Comparing models library(ggplot2) #Kaplan-Meier curve dataframe #Add a row of model name km
We see here that the Cox model is the most volatile with the most data and features. It is higher for lower values and drops down sharply when the time increases. The survival forest is of the lowest range and resembles Kaplan-Meier curve. The difference might be because of Survival forest having less rows. The essence of the plots is that there can be different approaches to the same concept of survival analysis and one may choose the technique based on one’s comfort and situation. A better data with processed data points and treated missing values might fetch us a better R2 and more stable curves. At the same time, they will help better in finding time to event cases such as knowing the time when a promotion’s effect dies down, knowing when tumors will develop and become significant and lots of other applications with a significant chunk of them being from medical science. Survival, as the name suggests, relates to surviving objects and is thus related to event occurrence in a completely different way than machine learning. It is important to know this technique to know more and more ways data can help us in solving problems, with time involved in this particular case. Hope this article serves the purpose of giving a glimpse of survival analysis and the feature rich packages available in R.Here is the complete code for the article:
# install.packages("survival") # Loading the package library("survival") #Dataset description ?pbc #Fitting the survival model survival_func=survfit(Surv(pbc$time,pbc$status == 2)~1) survival_func #Plot the survival model plot(survival_func) # Fit Cox Model Cox_model = coxph(Surv(pbc$time,pbc$status==2) ~.,data=pbc) summary(Cox_model) #Create a survival curve from the cox model Cox_curveAuthor Bio:
This article was contributed by Perceptive Analytics. Madhur Modi, Chaitanya Sagar, Vishnu Reddy and Saneesh Veetil contributed to this article.Perceptive Analytics provides data analytics, data visualization, business intelligence and reporting services to e-commerce, retail, healthcare and pharmaceutical industries. Our client roster includes Fortune 500 and NYSE listed companies in the USA and India.
Related