library(survival)
library(survminer)
library(ggplot2)
library(ggfortify)
library(plotly)Ejemplo de Análisis de Supervivencia usando los datos de AML
Introducción
Se presenta un análisis de supervivencia clásico para los datos de Leucemia Mielógena Aguda (AML) que son presentados en el libro de Miller (1982), y que fueron trabajados originalmente por Embury, Elias, Heller et al. (1977), que provienen de un ensayo clínico controlado para evaluar la eficacia del mantenimiento de la quimioterapia para esta enfermedad, luego de alcanzar la remisión de la misma a través de un tratamiento de quimioterapia, en el cual los pacientes que ingresaron al estudio fueron asignados aleatoriamente a uno de dos grupos: (i) Un grupo que recibió un quimioterapia de mantenimiento y, (ii) Otro grupo que no la recibió. El pobjetivo inicial del ensayo era ver si el mantenimiento en quimioterapia prolongaba el tiempo hasta la recaída por la enfermedad.
Los datos están disponibles en el paquete survival (Thernau, 2023), que contiene además las principales funciones de éste y la mayoría de los análisis de supervivencias en el lenguaje R (R Core Team, 2023), cuya versión utilizada para el mismo fue la 4.3.1 ejecutadas en la versión 2024.4.2.764 del Ambiente de Desarrollo Integrado (IDE), RStudio (RStudio Team, 2022), las cuales fueron complementadas 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).
Los datos fueron analizados mediante métodos clásicos de Análisis de Supervivencia, incluyendo:
- La estimación de la función de supervivencia a través del estimador de Kaplan & Meier (1958), o algunas de sus variantes (Fleming & Harrington, 1991).
- La comparación de las funciones de supervivencia a través del Test de los Logaritmos de los Rangos (Log-Rank Test) (Gehan, 1965) o sus variantes (Prentice, 1978), entre ellos en test de Peto & Peto (1972).
Análisis
Carga de los paquetes:
Lo primero que hay que hacer es cargar los paquetes utilizados para el análisis:
Visualización de los datos:
Una descripción de la base de datos puede obtenerse a través de:
help(aml)starting httpd help server ... donecon lo cual se despliega la descripción en la ventana de ayuda.
Si se desea visualizar los datos basta con escribir el nombre:
aml time status x
1 9 1 Maintained
2 13 1 Maintained
3 13 0 Maintained
4 18 1 Maintained
5 23 1 Maintained
6 28 0 Maintained
7 31 1 Maintained
8 34 1 Maintained
9 45 0 Maintained
10 48 1 Maintained
11 161 0 Maintained
12 5 1 Nonmaintained
13 5 1 Nonmaintained
14 8 1 Nonmaintained
15 8 1 Nonmaintained
16 12 1 Nonmaintained
17 16 0 Nonmaintained
18 23 1 Nonmaintained
19 27 1 Nonmaintained
20 30 1 Nonmaintained
21 33 1 Nonmaintained
22 43 1 Nonmaintained
23 45 1 NonmaintainedY si se quiere visualizar en el formato de un análisis de supervivencia, se debe usar la función Surv, donde para el caso de datos censurados por la derecha, basta con declarar el tiempo y el estatus de la censura:
Surv(aml$time,aml$status) [1] 9 13 13+ 18 23 28+ 31 34 45+ 48 161+ 5 5 8 8
[16] 12 16+ 23 27 30 33 43 45 Y si se quiere visualizar la información de acuerdo a los grupos de tratamiento de mantenimiento, se deben hacer un filtrado para cada grupo y usar posteriormente la función Surv:
Surv(aml[aml$x=='Maintained',]$time,aml[aml$x=='Maintained',]$status) [1] 9 13 13+ 18 23 28+ 31 34 45+ 48 161+Surv(aml[aml$x=='Nonmaintained',]$time,aml[aml$x=='Nonmaintained',]$status) [1] 5 5 8 8 12 16+ 23 27 30 33 43 45 Estimación de la función de supervivencia:
En este caso debemos estimar la función de supervivencia para cada grupo de la terapia de mantenimiento.
La opción por defecto es a través del estimador de Kaplan y Meier, para lo cual conviene guardar el ajuste en un objeto que llamaremos km:
km <- survfit(Surv(time,status)~x, data = aml)A partir del cual podemos obtener las medidas resumen, mediante:
kmCall: survfit(formula = Surv(time, status) ~ x, data = aml)
n events median 0.95LCL 0.95UCL
x=Maintained 11 7 31 18 NA
x=Nonmaintained 12 11 23 8 NAY podemos obtener además el detalle de la función de supervivencia en cada tiempo, mediante:
summary(km)Call: survfit(formula = Surv(time, status) ~ x, data = aml)
x=Maintained
time n.risk n.event survival std.err lower 95% CI upper 95% CI
9 11 1 0.909 0.0867 0.7541 1.000
13 10 1 0.818 0.1163 0.6192 1.000
18 8 1 0.716 0.1397 0.4884 1.000
23 7 1 0.614 0.1526 0.3769 0.999
31 5 1 0.491 0.1642 0.2549 0.946
34 4 1 0.368 0.1627 0.1549 0.875
48 2 1 0.184 0.1535 0.0359 0.944
x=Nonmaintained
time n.risk n.event survival std.err lower 95% CI upper 95% CI
5 12 2 0.8333 0.1076 0.6470 1.000
8 10 2 0.6667 0.1361 0.4468 0.995
12 8 1 0.5833 0.1423 0.3616 0.941
23 6 1 0.4861 0.1481 0.2675 0.883
27 5 1 0.3889 0.1470 0.1854 0.816
30 4 1 0.2917 0.1387 0.1148 0.741
33 3 1 0.1944 0.1219 0.0569 0.664
43 2 1 0.0972 0.0919 0.0153 0.620
45 1 1 0.0000 NaN NA NAQue podemos visualizar mediante su gráfico, que puede ser obtenido mediante:
plot(km,xlab="Semanas", ylab="Supervivencia", main="Estimador de Kaplan y Meier de aml para pacientes \n mantenidos y no mantenidos en quimioterapia",lty=c(1,2),mark.time=FALSE)
legend(75,0.9,legend=c("No mantenidos","Mantenidos"),lty=c(1,2))
Pero tambien podemos graficar otras funciones, como por ejemplo, la ocurrencia de los eventos:
plot(km, fun = "event", xlab="Semanas", ylab="Supervivencia", main="Ocurrencia de los eventos de aml para pacientes \n mantenidos y no mantenidos en quimioterapia",lty=c(1,2),mark.time=FALSE)
legend(75,0.6,legend=c("No mantenidos","Mantenidos"),lty=c(1,2))
La función de riesgo acumulada:
plot(km,fun = "cumhaz", xlab="Semanas", ylab="Supervivencia", main="Estimador de la función de riesgo acumulada de aml para \n pacientes mantenidos y no mantenidos en quimioterapia",lty=c(1,2),mark.time=FALSE)
legend(75,0.9,legend=c("No mantenidos","Mantenidos"),lty=c(1,2)) 
Gráficos usando el paquete survminer:
Los tres gráficos anteriores pueden ser obtenidos a través de la función ggsurvplot, obteniéndose gráficos más bonitos, para cada caso, mediante:
# Gráfico del estimador de Kaplan y Meier según tratamiento:
ggsurvplot(km, data = aml)
# Gráfico de los eventos acumulados según tratamiento usando el estmador de Kaplan y Meier:
ggsurvplot(km, data = aml, fun = "event")
# Gráfico de la función de riesgo acumulado estimada según tratamiento:
ggsurvplot(km, data = aml, fun = "cumhaz")
Y el gráfico para el estimador de Kaplan y Meier, puede aún mejorarse, mediante:
ggsurvplot(km, data = aml,
conf.int = TRUE,
pval = TRUE,
fun = "pct",
risk.table = TRUE,
size = 1,
linetype = "strata",
palette = c("#E7B800",
"#2E9FDF"),
legend = "bottom",
legend.title = "Estimador de KM según tratamiento adicional",
legend.labs = c("No mantenido",
"Mantenido"))
Gráficos usando ggplot2:
Gráfico usando ggfortify:
autoplot(km)
Gráficos usando ggfortify y ggplotly (gráfico interactivo):
km.plot <- autoplot(km)
ggplotly(km.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/issuesVariantes del estimador de Kaplan y Meier:
Existen algunas variantes al estimador de Kaplan y Meier, una de ellas es usando el estimador de Fleming y Harrington, obtenido a través del estimador de Nelson y Aalen para la función de riesgo acumulado, siendo los códigos y salidas correspondientes:
fh <- survfit(Surv(time,status)~x, data = aml, , type ="fleming-harrington")A partir del cual podemos obtener las medidas resumen, mediante:
fhCall: survfit(formula = Surv(time, status) ~ x, data = aml, type = "fleming-harrington")
n events median 0.95LCL 0.95UCL
x=Maintained 11 7 34 23 NA
x=Nonmaintained 12 11 27 8 NAY podemos obtener además el detalle de la función de supervivencia en cada tiempo, mediante:
summary(fh)Call: survfit(formula = Surv(time, status) ~ x, data = aml, type = "fleming-harrington")
x=Maintained
time n.risk n.event survival std.err lower 95% CI upper 95% CI
9 11 1 0.913 0.083 0.764 1.000
13 10 1 0.826 0.112 0.634 1.000
18 8 1 0.729 0.134 0.508 1.000
23 7 1 0.632 0.147 0.400 0.998
31 5 1 0.517 0.159 0.283 0.945
34 4 1 0.403 0.160 0.185 0.876
48 2 1 0.244 0.156 0.070 0.853
x=Nonmaintained
time n.risk n.event survival std.err lower 95% CI upper 95% CI
5 12 2 0.8465 0.0998 0.67190 1.000
8 10 2 0.6930 0.1276 0.48313 0.994
12 8 1 0.6116 0.1361 0.39543 0.946
23 6 1 0.5177 0.1439 0.30022 0.893
27 5 1 0.4239 0.1452 0.21663 0.829
30 4 1 0.3301 0.1400 0.14379 0.758
33 3 1 0.2365 0.1276 0.08219 0.681
43 2 1 0.1435 0.1055 0.03394 0.606
45 1 1 0.0528 0.0655 0.00463 0.601Que podemos visualizar mediante su gráfico, que puede ser obtenido mediante:
plot(fh,xlab="Semanas", ylab="Supervivencia", main="Estimador de Fleming y Harrington de aml para pacientes \n mantenidos y no mantenidos en quimioterapia",lty=c(1,2),mark.time=FALSE)
legend(75,0.9,legend=c("No mantenidos","Mantenidos"),lty=c(1,2))
Pero tambien podemos graficar otras funciones, como por ejemplo, la ocurrencia de los eventos:
plot(fh, fun = "event", xlab="Semanas", ylab="Supervivencia", main="Ocurrencia de los eventos de aml para pacientes \n mantenidos y no mantenidos en quimioterapia \n a través del estimador de Fleming y Harrington",lty=c(1,2),mark.time=FALSE)
legend(75,0.6,legend=c("No mantenidos","Mantenidos"),lty=c(1,2))
La función de riesgo acumulada:
plot(fh,fun = "cumhaz", xlab="Semanas", ylab="Supervivencia", main="Estimador de Nelson y Aalen de la función de riesgo acumulada de aml \n para pacientes mantenidos y no mantenidos en quimioterapia",lty=c(1,2),mark.time=FALSE)
legend(75,0.9,legend=c("No mantenidos","Mantenidos"),lty=c(1,2))
Gráficos usando el paquete survminer:
Los tres gráficos anteriores pueden ser obtenidos a través de la función ggsurvplot, obteniéndose gráficos más bonitos, para cada caso, mediante:
# Gráfico del estimador de Fleming y Harrington según tratamiento:
ggsurvplot(fh, data = aml) 
# Gráfico de los eventos acumulados según tratamiento usando el estmador de Fleming y Harrington:
ggsurvplot(fh, data = aml, fun = "event") 
# Gráfico de la función de riesgo acumulado estimada a través del estimador de Nelson y Aalen según tratamiento:
ggsurvplot(fh, data = aml, fun = "cumhaz")
Y el gráfico para el estimador de Fleming y Harrington, puede aún mejorarse, mediante:
ggsurvplot(fh, data = aml,
conf.int = TRUE,
pval = TRUE,
fun = "pct",
risk.table = TRUE,
size = 1,
linetype = "strata",
palette = c("#E7B800",
"#2E9FDF"),
legend = "bottom",
legend.title = "Estimador de FH según tratamiento adicional",
legend.labs = c("No mantenido",
"Mantenido"))
Gráficos usando ggplot2:
Gráfico usando ggfortify:
autoplot(fh)
Gráficos usando ggfortify y ggplotly (gráfico interactivo):
fh.plot <- autoplot(fh)
ggplotly(fh.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/issuesY podemos además usar una variante del estimador de Fleming y Harrington, mediante:
fh2 <- survfit(Surv(time,status)~x, data = aml, , type ="fh2")A partir del cual podemos obtener las medidas resumen, mediante:
fh2Call: survfit(formula = Surv(time, status) ~ x, data = aml, type = "fh2")
n events median 0.95LCL 0.95UCL
x=Maintained 11 7 34 23 NA
x=Nonmaintained 12 11 27 8 NAY podemos obtener además el detalle de la función de supervivencia en cada tiempo, mediante:
summary(fh2)Call: survfit(formula = Surv(time, status) ~ x, data = aml, type = "fh2")
x=Maintained
time n.risk n.event survival std.err lower 95% CI upper 95% CI
9 11 1 0.913 0.083 0.764 1.000
13 10 1 0.826 0.112 0.634 1.000
18 8 1 0.729 0.134 0.508 1.000
23 7 1 0.632 0.147 0.400 0.998
31 5 1 0.517 0.159 0.283 0.945
34 4 1 0.403 0.160 0.185 0.876
48 2 1 0.244 0.156 0.070 0.853
x=Nonmaintained
time n.risk n.event survival std.err lower 95% CI upper 95% CI
5 12 2 0.8401 0.1036 0.65971 1.000
8 10 2 0.6802 0.1318 0.46525 0.994
12 8 1 0.6003 0.1384 0.38200 0.943
23 6 1 0.5081 0.1446 0.29093 0.887
27 5 1 0.4160 0.1447 0.21041 0.823
30 4 1 0.3240 0.1388 0.13994 0.750
33 3 1 0.2322 0.1260 0.08013 0.673
43 2 1 0.1408 0.1039 0.03315 0.598
45 1 1 0.0518 0.0644 0.00453 0.592Que podemos visualizar mediante su gráfico, que puede ser obtenido mediante:
plot(fh2,xlab="Semanas", ylab="Supervivencia", main="Estimador de Fleming y Harrington modificado de aml para \n pacientes mantenidos y no mantenidos en quimioterapia",lty=c(1,2),mark.time=FALSE)
legend(75,0.9,legend=c("No mantenidos","Mantenidos"),lty=c(1,2))
Pero tambien podemos graficar otras funciones, como por ejemplo, la ocurrencia de los eventos:
plot(fh2, fun = "event", xlab="Semanas", ylab="Supervivencia", main="Ocurrencia de los eventos de aml para pacientes \n mantenidos y no mantenidos en quimioterapia \n a través del estimador de Fleming y Harrington modificado",lty=c(1,2),mark.time=FALSE)
legend(75,0.6,legend=c("No mantenidos","Mantenidos"),lty=c(1,2))
La función de riesgo acumulada:
plot(fh2,fun = "cumhaz", xlab="Semanas", ylab="Supervivencia", main="Estimador de Nelson y Aalen modificado de la función de riesgo acumulada \n de aml para pacientes mantenidos y no mantenidos en quimioterapia",lty=c(1,2),mark.time=FALSE)
legend(75,0.9,legend=c("No mantenidos","Mantenidos"),lty=c(1,2))
Gráficos usando el paquete survminer:
Los tres gráficos anteriores pueden ser obtenidos a través de la función ggsurvplot, obteniéndose gráficos más bonitos, para cada caso, mediante:
# Gráfico del estimador de Fleming y Harrington modificado según tratamiento:
ggsurvplot(fh2, data = aml) 
# Gráfico de los eventos acumulados según tratamiento usando el estmador de Fleming y Harrington modificado:
ggsurvplot(fh2, data = aml, fun = "event") 
# Gráfico de la función de riesgo acumulado estimada a través del estimador de Nelson y Aalen modificado según tratamiento:
ggsurvplot(fh2, data = aml, fun = "cumhaz")
Y el gráfico para el estimador de Flemin y Harrington modificado, puede aún mejorarse, mediante:
ggsurvplot(fh2, data = aml,
conf.int = TRUE,
pval = TRUE,
fun = "pct",
risk.table = TRUE,
size = 1,
linetype = "strata",
palette = c("#E7B800",
"#2E9FDF"),
legend = "bottom",
legend.title = "Estimador de FH modificado según tratamiento adicional",
legend.labs = c("No mantenido",
"Mantenido"))
Gráficos usando ggplot2:
Gráfico usando ggfortify:
autoplot(fh2)
Gráficos usando ggfortify y ggplotly (gráfico interactivo):
fh2.plot <- autoplot(fh2)
ggplotly(fh2.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 las funciones de supervivencia para ambos grupos
Para comparar las funciones de supervivencia del grupo de mantenidos en quimioterapia y el grupo de no mantenidos usamos la función survdiff, y el código correspondiente es:
survdiff(Surv(time,status)~x, data = aml)Call:
survdiff(formula = Surv(time, status) ~ x, data = aml)
N Observed Expected (O-E)^2/E (O-E)^2/V
x=Maintained 11 7 10.69 1.27 3.4
x=Nonmaintained 12 11 7.31 1.86 3.4
Chisq= 3.4 on 1 degrees of freedom, p= 0.07 Y la salida obtenida corresponde a la del Test de los Logaritmo de los Rangos (Log-Rank Test), a partir de la cual, podemos afirmar con un nivel de significación del 5% que no existe evidencia significativa de que las funciones de supervivencia sean distintas para ambos grupos (p-valor = 0,07).
Si se incluye y modifica el argumento rho (originalmente rho = 0), se pueden obtener una gran variedad de tests, por ejemplo, con rho = 1, se obtiene la salida del Test de Peto y Peto:
survdiff(Surv(time,status)~x, rho = 1, data = aml)Call:
survdiff(formula = Surv(time, status) ~ x, data = aml, rho = 1)
N Observed Expected (O-E)^2/E (O-E)^2/V
x=Maintained 11 3.85 6.14 0.859 2.78
x=Nonmaintained 12 7.18 4.88 1.081 2.78
Chisq= 2.8 on 1 degrees of freedom, p= 0.1 A partir de la cual, podemos afirmar con un nivel de significación del 5% que no existe evidencia significativa de que las funciones de supervivencia sean distintas para ambos grupos en los tiempos iniciales (p-valor = 0,1).
Referencias
Embury, S. H., Elias, L., Heller, P. H., Hood, C. E., Greenberg, P. L., & Schrier, S. L. (1977). Remission maintenance therapy in acute myelogenous leukemia. Western Journal of Medicine, 126(4): 267-272.
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.
Miller, R. G. (1981). Survival Anaysis. N.Y.: John Wiley & Sons.
Peto, R., & Peto, J. (1972). Asymptotically efficient rank invariant test procedures. Journal of the Royal Statistical Society: Series A (General), 135(2): 185-198.
Prentice, R. L. (1978). Linear rank tests with right censored data, Biometrika, 65 (1): 167–179.
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.