Étapes pour effectuer une analyse de survie dans R

 Étapes pour effectuer une analyse de survie dans R


[Cetarticleaétépubliépourlapremièrefoisle[Thisarticlewasfirstpublishedon[Cetarticleaétépubliépourlapremièrefoisle[Thisarticlewasfirstpublishedon R-posts.com, et a gentiment contribué aux R-blogueurs]. (Vous pouvez signaler un problème concernant le contenu de cette page ici)


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.0002020952

Ces 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_curve 

Author 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.





Source de l’article

A découvrir