Analyse de survie avec R | R-blogueurs
Vous souhaitez partager votre contenu sur les R-blogueurs? cliquez ici si vous avez un blog, ou ici si vous ne le faites pas.
Avec des racines remontant au moins à 1662 lorsque John Graunt, un marchand londonien, a publié un vaste ensemble d’inférences basées sur les enregistrements de mortalité, l’analyse de survie est l’un des sous-domaines les plus anciens de la statistique [1]. Les méthodes de base des tables de survie, y compris les techniques de traitement des données censurées, étaient connues avant 1700 [2]. Au début du XVIIIe siècle, les anciens maîtres, de Moivre travaillant sur les rentes et Daniel Bernoulli étudiant les risques concurrents pour ses travaux sur l’inoculation de la variole, ont développé les fondations de la modélisation du time-to-event. [2]. Aujourd’hui, les modèles d’analyse de survie sont importants dans l’ingénierie, l’assurance, le marketing et la médecine et bien d’autres domaines d’application. Il n’est donc pas surprenant que la vue des tâches R sur l’analyse de survie, une liste organisée, organisée et annotée de packages et de fonctions R pertinents, soit formidable.
Regarder la vue des tâches sur un petit écran, c’est un peu comme se tenir trop près d’un mur de briques – gauche-droite, haut-bas, des briques tout autour. C’est un édifice fantastique qui donne une idée des contributions importantes que les développeurs de R ont apportées à la fois à la théorie et à la pratique de l’analyse de survie. Aussi bien organisé soit-il, cependant, j’imagine que même les experts en analyse de survie ont besoin de temps pour trouver leur chemin dans cette vue des tâches. (Je m’en voudrais de ne pas mentionner que nous devons tous beaucoup de gratitude à Arthur Allignol et Aurielien Latouche, les responsables de la vue de tâches.) Les nouveaux arrivants, les personnes novices dans R ou novices dans l’analyse de survie ou les deux, doivent trouver cela écrasant. C’est donc en pensant aux nouveaux arrivants que j’offre la trajectoire mince suivante à travers la vue des tâches qui ne repose que sur quelques packages: survie, KMsurv, Oisurv et ranger
Le paquet de survie, qui a commencé sa vie comme un Paquet S à la fin des années 90, est la pierre angulaire de l’ensemble de l’édifice R Survival Analysis. Non seulement le package lui-même est riche en fonctionnalités, mais l’objet créé par le Surv() , qui contient le temps de panne et les informations de censure, est la structure de base des données d’analyse de survie dans R.
KMsurv contient des ensembles de données intéressants tirés du texte classique de John Klein et Melvin Moeschberger, Techniques d’analyse de survie pour les données censurées et tronquées.
Ma principale raison de choisir le package OIsurv était d’attirer l’attention sur le guide très utile Analyse de survie dans R, produit par les gens de OpenIntro.
ranger pourrait être la surprise dans ma très courte liste de packages de survie. le ranger() La fonction est bien connue pour être une implémentation rapide de l’algorithme Random Forests pour la construction d’ensembles d’arbres de classification et de régression. Mais c’était une nouvelle pour moi que ranger() construit également des modèles de survie. Benchmarks indique que ranger() convient à la création de modèles de temps à événement avec les grands ensembles de données de grande dimension importants pour les applications de marketing Internet. Puisque ranger() utilise la norme Surv() objets de survie, c’est un outil idéal pour se familiariser avec l’analyse de survie à l’ère de l’apprentissage automatique.
*** Charger et transformer les données Ce premier bloc de code charge les packages requis avec la trame de données de greffe de moelle osseuse du package KMsurv. J’ai choisi ceci car il y a un certain nombre de covariables et aucune valeur manquante. La seule préparation des données consiste à transformer les variables appropriées en facteurs.
library(survival)
library(dplyr)
library(OIsurv) # Aumatically loads KMsurv
library(ranger)
library(ggplot2)
#------------
data(bmt)
# sapply(bmt,class)
# Prepare new data frame for modeling
bmt2 <- select(bmt,-c(t2,d2,d3))
bmt2 <- mutate(bmt2,
group = as.factor(group),
d1 = as.factor(d1),
da = as.factor(da),
dc = as.factor(dc),
dp = as.factor(dp),
z3 = as.factor(z3),
z4 = as.factor(z4),
z5 = as.factor(z5),
z6 = as.factor(z6),
z8 = as.factor(z8),
z9 = as.factor(z9),
z10 = as.factor(z10)
)
head(bmt2)
## group t1 d1 ta da tc dc tp dp z1 z2 z3 z4 z5 z6 z7 z8 z9 z10 ## 1 1 2081 0 67 1 121 1 13 1 26 33 1 0 1 1 98 0 1 0 ## 2 1 1602 0 1602 0 139 1 18 1 21 37 1 1 0 0 1720 0 1 0 ## 3 1 1496 0 1496 0 307 1 12 1 26 35 1 1 1 0 127 0 1 0 ## 4 1 1462 0 70 1 95 1 13 1 17 21 0 1 0 0 168 0 1 0 ## 5 1 1433 0 1433 0 236 1 12 1 32 36 1 1 1 1 93 0 1 0 ## 6 1 1377 0 1377 0 123 1 12 1 22 31 1 1 1 1 2187 0 1 0
Analyse de Kaplan Meier
La première chose à faire est d'utiliser Surv() pour construire l'objet de survie standard. La variable t1 enregistre l'heure de la mort ou l'heure censurée; d1 indique que le patient est décédé (d1 = 1) ou que le patient a survécu jusqu'à la fin de l'étude (d1 = 0). Notez qu'un «+» après l'heure dans l'impression de y_bmt indique la censure. La formule y_bmt ~ 1 instruit le survfit() pour ajuster un modèle avec interception uniquement, et produit le Kaplan-Meier estimation.
le confBands() fonction des estimations du package OIsurv Hall-Wellner bandes de confiance. Il s'agit d'estimations simultanées à grand échantillon et représentées en rouge. Ils contrastent avec les bandes de confiance ponctuelles rendues sous forme de lignes pointillées noires.
# Kaplan Meier Survival Curve y_bmt <- Surv(bmt$t1, bmt$d1) y_bmt
## [1] 2081+ 1602+ 1496+ 1462+ 1433+ 1377+ 1330+ 996+ 226+ 1199+ 1111+ ## [12] 530+ 1182+ 1167+ 418 417 276 156 781 172 487 716 ## [23] 194 371 526 122 1279 110 243 86 466 262 162 ## [34] 262 1 107 269 350 2569+ 2506+ 2409+ 2218+ 1857+ 1829+ ## [45] 1562+ 1470+ 1363+ 1030+ 860+ 1258+ 2246+ 1870+ 1799+ 1709+ 1674+ ## [56] 1568+ 1527+ 1324+ 957+ 932+ 847+ 848+ 1850+ 1843+ 1535+ 1447+ ## [67] 1384+ 414 2204 1063 481 105 641 390 288 522 79 ## [78] 1156 583 48 431 1074 393 10 53 80 35 1499+ ## [89] 704 653 222 1356+ 2640+ 2430+ 2252+ 2140+ 2133+ 1238+ 1631+ ## [100] 2024+ 1345+ 1136+ 845+ 491 162 1298 121 2 62 265 ## [111] 547 341 318 195 469 93 515 183 105 128 164 ## [122] 129 122 80 677 73 168 74 16 248 732 105 ## [133] 392 63 97 153 363
fit1_bmt <- survfit(y_bmt ~ 1) summary(fit1_bmt)
## Call: survfit(formula = y_bmt ~ 1) ## ## time n.risk n.event survival std.err lower 95% CI upper 95% CI ## 1 137 1 0.993 0.00727 0.979 1.000 ## 2 136 1 0.985 0.01025 0.966 1.000 ## 10 135 1 0.978 0.01250 0.954 1.000 ## 16 134 1 0.971 0.01438 0.943 0.999 ## 35 133 1 0.964 0.01602 0.933 0.995 ## 48 132 1 0.956 0.01748 0.923 0.991 ## 53 131 1 0.949 0.01881 0.913 0.987 ## 62 130 1 0.942 0.02003 0.903 0.982 ## 63 129 1 0.934 0.02117 0.894 0.977 ## 73 128 1 0.927 0.02222 0.884 0.972 ## 74 127 1 0.920 0.02322 0.875 0.966 ## 79 126 1 0.912 0.02415 0.866 0.961 ## 80 125 2 0.898 0.02588 0.848 0.950 ## 86 123 1 0.891 0.02668 0.840 0.944 ## 93 122 1 0.883 0.02744 0.831 0.939 ## 97 121 1 0.876 0.02817 0.822 0.933 ## 105 120 3 0.854 0.03017 0.797 0.915 ## 107 117 1 0.847 0.03078 0.788 0.909 ## 110 116 1 0.839 0.03137 0.780 0.903 ## 121 115 1 0.832 0.03193 0.772 0.897 ## 122 114 2 0.818 0.03300 0.755 0.885 ## 128 112 1 0.810 0.03350 0.747 0.879 ## 129 111 1 0.803 0.03399 0.739 0.872 ## 153 110 1 0.796 0.03445 0.731 0.866 ## 156 109 1 0.788 0.03490 0.723 0.860 ## 162 108 2 0.774 0.03575 0.707 0.847 ## 164 106 1 0.766 0.03615 0.699 0.841 ## 168 105 1 0.759 0.03653 0.691 0.834 ## 172 104 1 0.752 0.03690 0.683 0.828 ## 183 103 1 0.745 0.03726 0.675 0.821 ## 194 102 1 0.737 0.03760 0.667 0.815 ## 195 101 1 0.730 0.03793 0.659 0.808 ## 222 100 1 0.723 0.03825 0.651 0.802 ## 243 98 1 0.715 0.03856 0.644 0.795 ## 248 97 1 0.708 0.03886 0.636 0.788 ## 262 96 2 0.693 0.03943 0.620 0.775 ## 265 94 1 0.686 0.03969 0.612 0.768 ## 269 93 1 0.678 0.03995 0.604 0.761 ## 276 92 1 0.671 0.04019 0.597 0.755 ## 288 91 1 0.664 0.04042 0.589 0.748 ## 318 90 1 0.656 0.04063 0.581 0.741 ## 341 89 1 0.649 0.04084 0.574 0.734 ## 350 88 1 0.642 0.04104 0.566 0.727 ## 363 87 1 0.634 0.04122 0.558 0.720 ## 371 86 1 0.627 0.04140 0.551 0.713 ## 390 85 1 0.619 0.04156 0.543 0.706 ## 392 84 1 0.612 0.04172 0.535 0.699 ## 393 83 1 0.605 0.04186 0.528 0.693 ## 414 82 1 0.597 0.04199 0.520 0.686 ## 417 81 1 0.590 0.04212 0.513 0.679 ## 418 80 1 0.583 0.04223 0.505 0.671 ## 431 79 1 0.575 0.04234 0.498 0.664 ## 466 78 1 0.568 0.04243 0.490 0.657 ## 469 77 1 0.560 0.04252 0.483 0.650 ## 481 76 1 0.553 0.04259 0.476 0.643 ## 487 75 1 0.546 0.04266 0.468 0.636 ## 491 74 1 0.538 0.04271 0.461 0.629 ## 515 73 1 0.531 0.04276 0.453 0.622 ## 522 72 1 0.524 0.04280 0.446 0.615 ## 526 71 1 0.516 0.04282 0.439 0.607 ## 547 69 1 0.509 0.04285 0.431 0.600 ## 583 68 1 0.501 0.04287 0.424 0.593 ## 641 67 1 0.494 0.04288 0.416 0.585 ## 653 66 1 0.486 0.04288 0.409 0.578 ## 677 65 1 0.479 0.04286 0.402 0.571 ## 704 64 1 0.471 0.04284 0.394 0.563 ## 716 63 1 0.464 0.04281 0.387 0.556 ## 732 62 1 0.456 0.04277 0.380 0.548 ## 781 61 1 0.449 0.04272 0.372 0.541 ## 1063 52 1 0.440 0.04276 0.364 0.533 ## 1074 51 1 0.432 0.04278 0.355 0.524 ## 1156 48 1 0.423 0.04282 0.346 0.515 ## 1279 42 1 0.413 0.04297 0.336 0.506 ## 1298 41 1 0.402 0.04308 0.326 0.496 ## 2204 9 1 0.358 0.05696 0.262 0.489
cb <- confBands(y_bmt, type = "hall")
plot(fit1_bmt,
main = 'Kaplan Meyer Plot with confidence bands')
lines(cb, col = "red",lty = 3)
legend(1000, 0.99, legend = c('K-M survival estimate',
'pointwise intervals', 'Hall-Werner conf bands'), lty = 1:3)
Modèle de risques proportionnels de Cox
Ensuite, je place un Risques proportionnels Cox modèle, qui utilise plusieurs des variables contenues dans l'ensemble de données. Je ne prétends pas avoir le meilleur modèle ici, ni même un très bon. Mon intention est simplement de montrer le package de survie coxph() et montrez comment les covariables entrent dans la formule. Notez, cependant, que ce modèle atteint un R2 valeur de 0,8, et que trois variables - ta, tc, et dc1 - apparaissent comme significatifs.
# Fit Cox Model
form <- formula(y_bmt ~ group + ta + tc + tp + dc + dp +
z1 + z2 + z3 + z4 + z5 + z6 + z7 + z8 + z10)
cox_bmt <- coxph(form,data = bmt2)
summary(cox_bmt)
## Call: ## coxph(formula = form, data = bmt2) ## ## n= 137, number of events= 81 ## ## coef exp(coef) se(coef) z Pr(>|z|) ## group2 0.2540084 1.2891826 0.4342118 0.585 0.5586 ## group3 0.5280287 1.6955865 0.4197516 1.258 0.2084 ## ta -0.0017962 0.9982054 0.0003736 -4.808 1.52e-06 *** ## tc -0.0060005 0.9940175 0.0010512 -5.708 1.14e-08 *** ## tp -0.0005709 0.9994293 0.0012622 -0.452 0.6510 ## dc1 -3.5832062 0.0277865 0.4938320 -7.256 3.99e-13 *** ## dp1 -1.2266468 0.2932744 0.4832251 -2.538 0.0111 * ## z1 0.0109077 1.0109674 0.0230293 0.474 0.6358 ## z2 -0.0155304 0.9845895 0.0200202 -0.776 0.4379 ## z31 0.0162146 1.0163468 0.2480277 0.065 0.9479 ## z41 0.0334753 1.0340419 0.2722755 0.123 0.9021 ## z51 0.0792485 1.0824733 0.2856339 0.277 0.7814 ## z61 -0.0702048 0.9322029 0.2557463 -0.275 0.7837 ## z7 0.0000440 1.0000440 0.0004235 0.104 0.9173 ## z81 0.4618598 1.5870228 0.3386689 1.364 0.1726 ## z101 0.5145619 1.6729055 0.3330479 1.545 0.1223 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## exp(coef) exp(-coef) lower .95 upper .95 ## group2 1.28918 0.7757 0.55044 3.01937 ## group3 1.69559 0.5898 0.74478 3.86023 ## ta 0.99821 1.0018 0.99747 0.99894 ## tc 0.99402 1.0060 0.99197 0.99607 ## tp 0.99943 1.0006 0.99696 1.00190 ## dc1 0.02779 35.9887 0.01056 0.07314 ## dp1 0.29327 3.4098 0.11375 0.75613 ## z1 1.01097 0.9892 0.96635 1.05764 ## z2 0.98459 1.0157 0.94670 1.02399 ## z31 1.01635 0.9839 0.62506 1.65258 ## z41 1.03404 0.9671 0.60642 1.76319 ## z51 1.08247 0.9238 0.61842 1.89474 ## z61 0.93220 1.0727 0.56470 1.53887 ## z7 1.00004 1.0000 0.99921 1.00087 ## z81 1.58702 0.6301 0.81716 3.08218 ## z101 1.67291 0.5978 0.87092 3.21338 ## ## Concordance= 0.926 (se = 0.034 ) ## Rsquare= 0.806 (max possible= 0.995 ) ## Likelihood ratio test= 224.6 on 16 df, p=0 ## Wald test = 107.5 on 16 df, p=1.332e-15 ## Score (logrank) test = 211.7 on 16 df, p=0
cox_fit_bmt <- survfit(cox_bmt) # plot(cox_fit_bmt)
Modèle de forêts aléatoires
Pour donner une idée de la portée des capacités de R à travailler avec les données Time to Event, j'utilise le ranger() pour adapter un modèle d'ensemble de forêts aléatoires aux données. Notez que ranger() crée un modèle pour chaque observation de l'ensemble de données. Le bloc de code suivant construit le modèle en utilisant les mêmes variables que le modèle de Cox ci-dessus et chacune des 137 courbes de survie calculées pour l'ensemble de données bmt, ainsi qu'une courbe de valeurs moyennes.
# ranger model
# ranger model
r_fit_bmt <- ranger(form,
data = bmt2,
importance = "permutation",
seed = 1234)
# Average the survival models
death_times <- r_fit_bmt$unique.death.times
surv_prob <- data.frame(r_fit_bmt$survival)
avg_prob <- sapply(surv_prob,mean)
# Plot the survival models for each patient
plot(r_fit_bmt$unique.death.times,r_fit_bmt$survival[1,], type = "l",
ylim = c(0,1),
col = "red",
xlab = "death times",
ylab = "survival",
main = "Patient Survival Curves")
for(n in c(2:137)){
lines(r_fit_bmt$unique.death.times, r_fit_bmt$survival[n,], type = "l", col = "red")
}
lines(death_times, avg_prob, lwd = 2)
legend(100, 0.2, legend = c('Averages - black'))
Ici, nous montrons le classement de l'importance des variables calculé par la méthode de permutation, qui est ranger()Par défaut pour les données de survie. Notez que ta, tc, et dc sont les mêmes trois principales variables marquées dans le modèle de Cox.
Une mesure de l'erreur de prédiction calculée à partir de C-index de Harrell. Cet indice est défini comme «… la proportion de toutes les paires de patients utilisables dans lesquelles les prévisions et les résultats sont concordants» (cf. [8] p 370), où les prédictions pour les couples sont concordantes si les temps de survie prédits sont plus grands pour les patients qui ont vécu plus longtemps. Notez que l’indice c de Harrell peut être considéré comme une généralisation de la recherche du sont sous une courbe ROC. (Pour les résultats binaires, le c-index de Harrell se réduit à Statistique de Wilcoxon-Mann-Whitney ce qui, à son tour, équivaut au calcul de l'aire sous la courbe ROC.)
vi <- data.frame(sort(round(r_fit_bmt$variable.importance, 4), decreasing = TRUE)) names(vi) <- "importance" head(vi)
## importance ## ta 0.1259 ## tc 0.0688 ## dc 0.0190 ## tp 0.0117 ## dp 0.0092 ## z2 0.0046
cat("Prediction Error = 1 - Harrell's c-index = ", r_fit_bmt$prediction.error)
## Prediction Error = 1 - Harrell's c-index = 0.09304771
Enfin, nous traçons les courbes de survie calculées pour les trois modèles sur le même graphique. À noter que la courbe «ad hoc» des courbes de survie moyennes calculées par le modèle Ranger suit assez bien la courbe de Kaplan-Meier.
# Set up for ggplot
km <- rep("KM", length(fit1_bmt$time))
km_df <- data.frame(fit1_bmt$time,fit1_bmt$surv,km)
names(km_df) <- c("Time","Surv","Model")
cox <- rep("Cox",length(cox_fit_bmt$time))
cox_df <- data.frame(cox_fit_bmt$time,cox_fit_bmt$surv,cox)
names(cox_df) <- c("Time","Surv","Model")
rf <- rep("RF",length(r_fit_bmt$unique.death.times))
rf_df <- data.frame(r_fit_bmt$unique.death.times,avg_prob,rf)
names(rf_df) <- c("Time","Surv","Model")
plot_df <- rbind(km_df,cox_df,rf_df)
p <- ggplot(plot_df, aes(x = Time, y = Surv, color = Model))
p + geom_line() + ggtitle("Comparison of Survival Curves")
Pour une très belle présentation du type de modélisation d'analyse prédictive de survie qui peut être réalisée avec ranger, n'oubliez pas de jeter un œil à Manuel Amunategui Publier et vidéo.
Cette excursion en quatre packages ne fait allusion qu'aux outils d'analyse de survie disponibles dans R, mais elle illustre une partie de la richesse de la plate-forme R qui est en développement et en amélioration continus depuis près de vingt ans. L'utilisation du Surv() La fonction montre comment le code open source permet à des générations de développeurs de s'appuyer sur le travail de leurs prédécesseurs. le ranger packages fournit un exemple pratique de la manière dont R peut incorporer du code C ++ rapide et s'adapter au monde des applications d'apprentissage automatique, et l'utilisation accidentelle d'options telles que les bandes de confiance Hall-Wellner et le c-index de Harrell donne une idée de la profondeur statistique qui sous-tend presque tout R.
Les références
Pour plus de commodité, j'ai rassemblé les références utilisées tout au long de l'article ici.
[1] Piratage, Ian. (2006) L'émergence de la probabilité: une étude philosophique des premières idées sur l'induction des probabilités et l'inférence statistique. Cambridge University Press, 2e éd. p11[2] Andersen, P.K., Keiding, N. (1998) Analyse de survie (vue d'ensemble) Encyclopédie de la biostatistique 6. Wiley, p 4452-4461
[3] Kaplan, E.L. Et Meier, P. (1958). Estimation non paramétrique des observations incomplètes, J American Stats Assn. 53, 457–481, 562–563.
[4] Cox, D.R. (1972). Modèles de régression et tables de mortalité (avec discussion), Journal of the Royal Statistical Society (B) 34, 187–220.
[5] Diez, David. Analyse de survie dans R. OpenIntro [6] Klein, John P et Moeschberger, Melvin L. (1997) Techniques d'analyse de survie pour les données censurées et tronquées, Springer.
[7] Wright, Marvin et Ziegler, Andreas. (2017) * [ranger: A Fast Implementation of Random Forests for High Dimensional Data in C++ and R] (https://www.jstatsoft.org/article/view/v077i01)*, JSS Vol 77, Numéro 1.
[8] Harrell, Frank, Lee, Kerry et Mark, Daniel. (1996) Modèles pronostiques multivariables: problèmes liés à l'élaboration de modèles, à l'évaluation des hypothèses et de l'adéquation, ainsi qu'à la mesure et à la réduction des erreurs. Statistiques en médecine, Vol 15, 361-387
[9] Amunategui, Manuel. Ensembles de survie: classification Survival Plus pour des prédictions temporelles améliorées dans R
en relation
