library(survival)
library(survminer)
library(ggplot2)
library(ggfortify)
library(plotly)Análisis de Supervivencia para pacientes en diálisis peritoneal: Estimadores y Comparaciones
Introducción
Se presenta la primera parte de un Análisis de Supervivencia Clásico para los datos del artículo de Borges (2005), Este análisis consiste en la estimación de la función de supervivencia a través del estimador de Kaplan & Meier (1958), o algunas variantes (Fleming & Harrington, 1991), presentando diversas alternativas gráficas, y la comparación de curvas de supervivencia a través del Test de los Logaritmos de los Rangos (Log Rank Test) (Gehan, 1965) y el Test de Peto & Peto (1972).
El análisis se llevó a cabo usando funciones del paquete survival (Thernau, 2023), en la versión 4.3.1 del lenguaje R (R Core Team, 2023), ejecutadas en la versión 2024.4.2.764 del Ambiente de Desarrollo Integrado (IDE), RStudio (RStudio Team, 2022), y los gráficos fueron complementados con funciones del paquete survminer (Kassambara, Kosinski & Biecek, 2021). o con funciones de ggplot2 (Wickham, 2016), ggfortify (Tang, Horikoshi & Li, 2016; Horikoshi & Tang, 2016) y plotly (Sievert, 2020).
Aspectos preliminares
Carga de los paquetes
Lo primero que hay que hacer es cargar los paquetes utilizados para el análisis:
Lectura de los datos
Los datos a analizar, se encuetran disponibles en el archivo dpa.txt, que luego de ser descargado de enlace , debe ser copiado en la carpeta de trabajo, para proceder a la lectura de los mismos:
dpa<-read.table("dpa.txt",header=TRUE)Análisis de supervivencia
Estimadores
Existen varios estimadores para la función de supervivencia, siendo el más conocido y utilizado el de Kaplan & Meier (1958), pero para el cual existen otras opciones como el de Fleming y Harrington, o el de Fleming y Harrington modificado (Fleming & Harrington, 1991).
Estimador de Kaplan y Meier
Ajuste del modelo
# Ajuste:
km1<-survfit(Surv(meses,censor2)~1, data = dpa)
# Resumen del modelo
print(km1)Call: survfit(formula = Surv(meses, censor2) ~ 1, data = dpa)
n events median 0.95LCL 0.95UCL
[1,] 246 64 61 55 NA# Función de supervivencia (tabla):
summary(km1)Call: survfit(formula = Surv(meses, censor2) ~ 1, data = dpa)
time n.risk n.event survival std.err lower 95% CI upper 95% CI
0 246 2 0.992 0.00573 0.9807 1.000
1 240 1 0.988 0.00704 0.9740 1.000
3 228 4 0.970 0.01102 0.9490 0.992
4 221 1 0.966 0.01182 0.9431 0.989
5 215 1 0.962 0.01259 0.9372 0.987
6 209 1 0.957 0.01334 0.9311 0.983
7 202 1 0.952 0.01409 0.9250 0.980
8 197 1 0.947 0.01483 0.9187 0.977
9 193 1 0.942 0.01554 0.9125 0.973
10 188 1 0.937 0.01625 0.9061 0.970
11 180 3 0.922 0.01831 0.8866 0.958
12 171 1 0.916 0.01898 0.8800 0.954
13 161 1 0.911 0.01970 0.8729 0.950
14 151 4 0.887 0.02257 0.8435 0.932
15 144 1 0.880 0.02324 0.8361 0.927
17 135 3 0.861 0.02532 0.8127 0.912
19 124 1 0.854 0.02605 0.8044 0.907
20 119 2 0.840 0.02752 0.7873 0.895
21 115 3 0.818 0.02956 0.7617 0.878
22 110 3 0.795 0.03143 0.7361 0.859
23 104 1 0.788 0.03205 0.7274 0.853
25 94 1 0.779 0.03278 0.7177 0.846
26 90 1 0.771 0.03354 0.7077 0.839
28 81 1 0.761 0.03445 0.6966 0.832
30 78 2 0.742 0.03623 0.6739 0.816
31 75 4 0.702 0.03933 0.6291 0.784
33 63 1 0.691 0.04025 0.6164 0.775
34 59 1 0.679 0.04124 0.6031 0.765
37 53 2 0.654 0.04348 0.5737 0.745
39 50 1 0.641 0.04453 0.5589 0.734
42 43 1 0.626 0.04592 0.5418 0.722
47 38 1 0.609 0.04757 0.5227 0.710
52 33 1 0.591 0.04958 0.5011 0.696
55 31 1 0.572 0.05152 0.4791 0.682
59 26 2 0.528 0.05616 0.4283 0.650
60 22 1 0.504 0.05850 0.4012 0.632
61 21 1 0.480 0.06044 0.3748 0.614
65 18 1 0.453 0.06268 0.3455 0.594
77 11 1 0.412 0.06920 0.2963 0.573
96 7 1 0.353 0.08054 0.2258 0.552
110 4 2 0.177 0.09701 0.0601 0.518# Función de supervivencia (gráfico):
plot(km1,xlab="Meses",ylab="Supervivencia", main="Gráfico No. 1. Estimador de Kaplan y Meier")
Variantes del estimador de Kaplan y Meier: Estimadores de Fleming y Harrington
# Estimador de Fleming y Harrington:
fh <- survfit(Surv(meses,censor2)~1, type ="fleming-harrington", data = dpa)
# Estimador de Fleming y Harrington:
fh.2 <- survfit(Surv(meses,censor2)~1, type ="fh2", data =dpa)
# Gráfico comparativo de los estimadores:
plot(km1,xlab="Meses", conf.int= FALSE,ylab="Supervivencia", main="Gráfico No. 2. Estimadores de la función de supervivencia")
lines(fh, lty=2, conf.int= FALSE)
lines(fh.2, lty=3, conf.int= FALSE)
legend(75,1,legend=c("KM","FH", "FH2"),lty=c(1,2,3))
Gráficos con el paquete survminer
Gráfico del estimador de Kaplan y Meier:
ggsurvplot(km1, data = dpa)
Gráfico de la historia de ocurrencia de los eventos:
ggsurvplot(km1, data = dpa, fun = "event")
Gráfico de la función de riesgo acumulada:
ggsurvplot(km1, data = dpa, fun = "cumhaz")
Gráficos con ggplot
Gráfico con ggfortify:
autoplot(km1)
Gráfico con ggfortify y plotly (gráfico interactivo):
km1.plot <- autoplot(km1)
ggplotly(km1.plot)Warning in geom2trace.default(dots[[1L]][[1L]], dots[[2L]][[1L]], dots[[3L]][[1L]]): geom_GeomConfint() has yet to be implemented in plotly.
If you'd like to see this geom implemented,
Please open an issue with your example code at
https://github.com/ropensci/plotly/issuesComparación de funciones de supervivencia
Estimación del modelo:
# Ajuste del modelo:
Km2 <- survfit(Surv(meses,censor2)~diabetes, data = dpa)
# Gráfico del modelo:
plot(Km2,xlab="Meses",ylab="Supervivencia", main="Gráfico No. 2. Estimador de Kaplan y Meier \n para pacientes diab?ticos y no diab?ticos",lty=c(1,2),mark.time=FALSE)
legend(75,0.9,legend=c("No diabéticos","Diabéticos"),lty=c(1,2))
# Medidas resumen del modelo:
print(Km2)Call: survfit(formula = Surv(meses, censor2) ~ diabetes, data = dpa)
n events median 0.95LCL 0.95UCL
diabetes=0 212 49 77 59 NA
diabetes=1 34 15 30 23 NA# Tabla de los detalles del modelo:
summary(Km2)Call: survfit(formula = Surv(meses, censor2) ~ diabetes, data = dpa)
diabetes=0
time n.risk n.event survival std.err lower 95% CI upper 95% CI
0 212 2 0.991 0.00664 0.978 1.000
3 198 4 0.971 0.01185 0.948 0.994
4 191 1 0.965 0.01283 0.941 0.991
5 185 1 0.960 0.01378 0.934 0.988
6 179 1 0.955 0.01471 0.926 0.984
7 173 1 0.949 0.01563 0.919 0.981
8 168 1 0.944 0.01653 0.912 0.977
9 165 1 0.938 0.01739 0.905 0.973
10 160 1 0.932 0.01824 0.897 0.969
11 154 2 0.920 0.01991 0.882 0.960
12 146 1 0.914 0.02075 0.874 0.955
14 129 4 0.885 0.02447 0.839 0.935
15 122 1 0.878 0.02532 0.830 0.929
17 114 1 0.870 0.02624 0.820 0.923
19 105 1 0.862 0.02727 0.810 0.917
20 100 1 0.854 0.02833 0.800 0.911
21 98 1 0.845 0.02935 0.789 0.904
22 95 3 0.818 0.03221 0.757 0.884
25 81 1 0.808 0.03336 0.745 0.876
26 78 1 0.798 0.03450 0.733 0.868
31 67 3 0.762 0.03863 0.690 0.842
34 54 1 0.748 0.04041 0.673 0.831
37 49 2 0.717 0.04415 0.636 0.809
39 46 1 0.702 0.04586 0.617 0.798
42 39 1 0.684 0.04809 0.596 0.785
47 34 1 0.664 0.05070 0.571 0.771
52 30 1 0.642 0.05362 0.545 0.756
55 28 1 0.619 0.05639 0.517 0.740
59 24 1 0.593 0.05964 0.487 0.722
60 21 1 0.565 0.06313 0.453 0.703
61 20 1 0.536 0.06598 0.421 0.683
65 17 1 0.505 0.06924 0.386 0.660
77 11 1 0.459 0.07666 0.331 0.637
96 7 1 0.393 0.08945 0.252 0.614
110 4 2 0.197 0.10803 0.067 0.577
diabetes=1
time n.risk n.event survival std.err lower 95% CI upper 95% CI
1 34 1 0.971 0.0290 0.9154 1.000
11 26 1 0.933 0.0460 0.8473 1.000
13 23 1 0.893 0.0593 0.7838 1.000
17 21 2 0.808 0.0784 0.6678 0.977
20 19 1 0.765 0.0850 0.6154 0.951
21 17 2 0.675 0.0959 0.5110 0.892
23 14 1 0.627 0.1005 0.4579 0.858
28 11 1 0.570 0.1063 0.3954 0.821
30 10 2 0.456 0.1115 0.2824 0.736
31 8 1 0.399 0.1112 0.2311 0.689
33 7 1 0.342 0.1089 0.1832 0.638
59 2 1 0.171 0.1326 0.0374 0.782Gráficos con survminer:
Gráfico del estimador de Kaplan y Meier según diábetes:
ggsurvplot(Km2, data = dpa)
Gráfico de la historia de ocurrencia de los eventos:
ggsurvplot(Km2, data = dpa, fun = "event")
Gráfico de la función de riesgo acumulada:
ggsurvplot(Km2, data = dpa, fun = "cumhaz")
Un gráfico aún mas bonito:
ggsurvplot(Km2, data = dpa,
conf.int = TRUE,
pval = TRUE,
fun = "pct",
risk.table = TRUE,
size = 1,
linetype = "strata",
palette = c("#E7B800",
"#2E9FDF"),
title = "Estimador de Kaplan-Meier para pacientes en diálisis peritoneal según diabetes",
legend = "bottom",
legend.title = "Estimador de KM según diabetes",
legend.labs = c("No diabético",
"Diabético"))
Gráficos con ggplot:
Gráfico con ggfortify:
autoplot(Km2)
Gráfico con ggfortify y plotly (gráfico interactivo):
km2.plot <- autoplot(Km2)
ggplotly(km2.plot)Warning in geom2trace.default(dots[[1L]][[2L]], dots[[2L]][[1L]], dots[[3L]][[1L]]): geom_GeomConfint() has yet to be implemented in plotly.
If you'd like to see this geom implemented,
Please open an issue with your example code at
https://github.com/ropensci/plotly/issues
Warning in geom2trace.default(dots[[1L]][[2L]], dots[[2L]][[1L]], dots[[3L]][[1L]]): geom_GeomConfint() has yet to be implemented in plotly.
If you'd like to see this geom implemented,
Please open an issue with your example code at
https://github.com/ropensci/plotly/issuesComparación de funciones de supervivencia:
Test de los Logaritmos de los Rangos (Log-Rank Test):
survdiff(Surv(meses,censor2)~diabetes, data = dpa)Call:
survdiff(formula = Surv(meses, censor2) ~ diabetes, data = dpa)
N Observed Expected (O-E)^2/E (O-E)^2/V
diabetes=0 212 49 56.46 0.986 8.6
diabetes=1 34 15 7.54 7.386 8.6
Chisq= 8.6 on 1 degrees of freedom, p= 0.003 Variantes:
# Log-Rank Test (opción por defecto:
survdiff(Surv(meses,censor2)~diabetes, rho = 0, data = dpa)Call:
survdiff(formula = Surv(meses, censor2) ~ diabetes, data = dpa,
rho = 0)
N Observed Expected (O-E)^2/E (O-E)^2/V
diabetes=0 212 49 56.46 0.986 8.6
diabetes=1 34 15 7.54 7.386 8.6
Chisq= 8.6 on 1 degrees of freedom, p= 0.003 # Test de Peto y Peto
survdiff(Surv(meses,censor2)~diabetes, rho = 1, data = dpa)Call:
survdiff(formula = Surv(meses, censor2) ~ diabetes, data = dpa,
rho = 1)
N Observed Expected (O-E)^2/E (O-E)^2/V
diabetes=0 212 38.5 44.34 0.774 7.28
diabetes=1 34 12.2 6.39 5.375 7.28
Chisq= 7.3 on 1 degrees of freedom, p= 0.007 # Otra variante:
survdiff(Surv(meses,censor2)~diabetes, rho = 0.5, data = dpa)Call:
survdiff(formula = Surv(meses, censor2) ~ diabetes, data = dpa,
rho = 0.5)
N Observed Expected (O-E)^2/E (O-E)^2/V
diabetes=0 212 43.1 49.67 0.879 7.98
diabetes=1 34 13.5 6.92 6.308 7.98
Chisq= 8 on 1 degrees of freedom, p= 0.005 Referencias:
Borges, R. E. (2005). Análisis de supervivencia de pacientes con diálisis peritoneal. Revista Colombiana de Estadística, 28(2), 243-259.
Fleming, T. R., & Harrington, D. P. (1991). Counting Processes and Survival Analysis. N.Y: John Wiley & Sons.
Gehan, E. A. (1965). A generalized Wilcoxon test for comparing arbitrarily singly-censored samples, Biometrika, 52: 203–223.
Horikoshi, M., & Tang, Y. (2016). ggfortify: Data Visualization Tools for Statistical Analysis Results. URL: https://CRAN.R-project.org/package=ggfortify.
Kaplan, E. L., & Meier, P. (1958). Nonparametric estimation from incomplete observations. Journal of the American statistical association, 53(282), 457-481.
Kassambara A, Kosinski M, & Biecek P (2021). survminer: Drawing Survival Curves using ‘ggplot2’_. R package version 0.4.9, URL: https://CRAN.R-project.org/package=survminer.
Peto, R., & Peto, J. (1972). Asymptotically efficient rank invariant test procedures. Journal of the Royal Statistical Society: Series A (General), 135(2): 185-198.
R Core Team (2023). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. URL: https://www.R-project.org/
RStudio Team (2024). RStudio: Integrated Development Environment for R. RStudio, PBC, Boston, MA URL: http://www.rstudio.com/.
Sievert. C. (2020). Interactive Web-Based Data Visualization with R, plotly, and shiny. Boca Ratón, FLA: Chapman and Hall/CRC.
Tang, Y., Horikoshi, M., & Li, W. (2016). ggfortify: Unified Interface to Visualize Statistical Result of Popular R Packages. The R Journal 8(2): 478-489.
Therneau T (2023). A Package for Survival Analysis in R. R package version 3.5-5, URL: https://CRAN.R-project.org/package=survival. }
Wickham, H. (2016). ggplot2: Elegant Graphics for Data Analysis. N.Y.: Springer-Verlag.