Análisis de supervivencia
Lina Angélica Buitrago Reyes - lina.buitrago2@unisabana.edu.co
1 Introducción
El análisis de supervivencia se diseño con el objetivo de analizar, en estudios longitudinales, los tiempos hasta que ocurra un evento de interés, no necesariamente la muerte.
Un ejemplo podría ser un estudio de cohorte en el cual se evaluen dos intervenciones de salud relacionadas con estilo de vida y el evento de interés sea “alcanzar el peso ideal”; así estaríamos analizando el tiempo (en días o meses o años) desde el inicio de la intervención hasta el momento en que el paciente alcanza su peso ideal. ¿En este caso cuál tratamiento resultaría más efectivo?
2 Objetivos
Formalizar los conocimientos sobre análisis de supervivencia, definiendo los conceptos básicos como función de supervivencia, hazard, etc.
Comprender los métodos paramétricos y no paramétricos para la estimación de la función de supervivencia.
Entender,aplicar y validar apropiadamente el modelo semiparamétrico de Cox.
3 Competencias
En este módulo el estudiante conocerá los diferentes métodos para estimar funciones de supervivencia y su aplicación en la práctica.
4 Conceptos preliminares
4.1 Datos censurados
En muchas ocasiones, no es posible observar la ocurrencia del evento en todos los participantes durante el periodo del estudio, cuando esto ocurre decimos que el dato no observado es un datos censurado. Esto puede ocurrir, en general, por las siguientes razones:
Debido a la longitud del periodo de seguimiento: El evento no ocurre antes de que el estudio finalice.
Por las pérdidas en el seguimiento durante el estudio.
El participante sale del estudio por otras razones (fallece y el evento de interés no es la muerte, reacciones adversas, etc).
4.1.1 Ejemplo: Tomado de Kleinbaum, D. G., & Klein, M. (2010)
Supongamos que se quiere evaluar el tiempo hasta la recaída (en semanas) en paciente en remisión de leucemia. La siguiente gráfica representa lo ocurrido con los 5 primeros pacientes. (X representa la ocurrencia del evento)
4.1.2 Tipos de censura
Censura a la derecha: El verdadero tiempo de supervivencia es mayor o igual al tiempo observado.
Censura a la izquierda: El verdadero tiempo de supervivencia es menor o igual al tiempo observado.
Censura de intervalo: El verdadero tiempo de supervivencia está dentro de un intervalo conocido.
Notación: Definiremos la v.a. \(d\) como la indicadora de la ocurrencia del evento, así:
\[d=1,\text{ si ocurre el evento}\] \[d=0,\text{ si el dato es censurado}\]
4.2 Función de supervivencia
Dadas \(T\) la variable que representa el tiempo hasta el evento, \(t\) un posible valor de dicha variable, \(f_T(t)\) la función de densidad de \(T\) y \(F_T(t)\) su correspondiente función de distribución, la función de supervivencia en el tiempo \(t\) se define como la probabilidad de que el tiempo al evento sea mayor a \(t\):
\[S(t)=P(T>t)=1−F_T(t),\text{ }0\leq t<\infty\]
4.2.1 Propiedades
- \(S(0)=1-F_T(0)=1\)
- Es una función no creciente: una curva describiría el comportamiento poblacional (izq) y una función escalonada su comportamiento muestral (der).
- \(S(\infty)=1-F_T(\infty)=0\)
4.3 Función de peligro o tasa de falla (\(h(t)\))
La función de peligro cuantifica el riesgo instantáneo, en momento determinado, de que ocurra el evento por unidad de tiempo, dado que el individuo a sobrevivido hasta un tiempo determinado:
\[h(t)=\lim_{Δt→0}\frac{P(t<T<t+Δt|T≥t)}{Δt}=\frac{f_T(t)}{S(t)}\]
4.3.1 Propiedades
- \(h(t)\geq 0\)
- Como resultado se obtiene la tasa media de fallas instantánea por unidad de tiempo, es decir que su magnitud depende de la unidad de tiempo en la que se esté trabajando, entre más pequeña sea esa unidad de tiempo, más grande será la tasa de falla.
4.4 Relación entre \(S(t)\) y \(h(t)\)
\[S(t)=exp\left[-\int_{u=0}^{u=t} h(u)du\right]\] \[h(t)=-\frac{S'(t)}{S(t)}\] Ejemplo: Si la tasa de falla instantánea es constante, \(h(t)=\lambda\):
\[S(t)=e^{-\lambda}\] # Estimador de Kaplan-Meier
4.5 Ejemplo
## Loading required package: ggplot2## Loading required package: ggpubr## Loading required package: magrittr#Tiempos
t<-c(6,6,6,7,10,13,16,22,23,6,9,10,11,17,19,20,25,32,32,34,35,1,1,2,2,3,4,4,5,5,8,8,8,8,11,11,12,12,15,17,22,23)
#Indicadora del evento. 1:ocurre el evento, 0:censura
d<-c(rep(1,9),rep(0,12),rep(1,21))
grupo<-c(rep(1,21),rep(2,21))
datos<-data.frame(t,d,grupo)
#Declaración de los datos como tipo Surv
surv<-Surv(t, d)
#Estimación de Kaplan-Meier para todos los datos
surv.km <- survfit(surv ~ 1, type = "kaplan-meier")
summary(surv.km)## Call: survfit(formula = surv ~ 1, type = "kaplan-meier")
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 42 2 0.952 0.0329 0.8901 1.000
## 2 40 2 0.905 0.0453 0.8202 0.998
## 3 38 1 0.881 0.0500 0.7883 0.985
## 4 37 2 0.833 0.0575 0.7279 0.954
## 5 35 2 0.786 0.0633 0.6709 0.920
## 6 33 3 0.714 0.0697 0.5899 0.865
## 7 29 1 0.690 0.0715 0.5628 0.845
## 8 28 4 0.591 0.0764 0.4588 0.762
## 10 23 1 0.565 0.0773 0.4325 0.739
## 11 21 2 0.512 0.0788 0.3783 0.692
## 12 18 2 0.455 0.0796 0.3227 0.641
## 13 16 1 0.426 0.0795 0.2958 0.615
## 15 15 1 0.398 0.0791 0.2694 0.588
## 16 14 1 0.369 0.0784 0.2437 0.560
## 17 13 1 0.341 0.0774 0.2186 0.532
## 22 9 2 0.265 0.0765 0.1507 0.467
## 23 7 2 0.189 0.0710 0.0909 0.395#Gráfica
ggsurvplot(fit = surv.km, data=datos,conf.int = T, title = "Estimación de la curva de supervivencia",
xlab = "t", ylab = "S(t)", legend.title = "Estimación",
legend.labs = "Kaplan-Meier")
#Para encontrar las estadísticas descriptivas del tiempo de supervivencia
print(surv.km,print.rmean = T)## Call: survfit(formula = surv ~ 1, type = "kaplan-meier")
##
## n events *rmean *se(rmean) median 0.95LCL
## 42.00 30.00 15.34 1.86 12.00 8.00
## 0.95UCL
## 22.00
## * restricted mean with upper limit = 35#Estimación de Kaplan-Meier para cada grupo
survgr.km <- survfit(surv ~ grupo, type = "kaplan-meier")
summary(survgr.km)## Call: survfit(formula = surv ~ grupo, type = "kaplan-meier")
##
## grupo=1
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 6 21 3 0.857 0.0764 0.720 1.000
## 7 17 1 0.807 0.0869 0.653 0.996
## 10 15 1 0.753 0.0963 0.586 0.968
## 13 12 1 0.690 0.1068 0.510 0.935
## 16 11 1 0.627 0.1141 0.439 0.896
## 22 7 1 0.538 0.1282 0.337 0.858
## 23 6 1 0.448 0.1346 0.249 0.807
##
## grupo=2
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 21 2 0.9048 0.0641 0.78754 1.000
## 2 19 2 0.8095 0.0857 0.65785 0.996
## 3 17 1 0.7619 0.0929 0.59988 0.968
## 4 16 2 0.6667 0.1029 0.49268 0.902
## 5 14 2 0.5714 0.1080 0.39455 0.828
## 8 12 4 0.3810 0.1060 0.22085 0.657
## 11 8 2 0.2857 0.0986 0.14529 0.562
## 12 6 2 0.1905 0.0857 0.07887 0.460
## 15 4 1 0.1429 0.0764 0.05011 0.407
## 17 3 1 0.0952 0.0641 0.02549 0.356
## 22 2 1 0.0476 0.0465 0.00703 0.322
## 23 1 1 0.0000 NaN NA NA#Gráfica
ggsurvplot(fit = survgr.km, data=datos,conf.int = T, title = "Estimación de la curva de supervivencia",
xlab = "t", ylab = "S(t)", legend.title ="Estimación",
legend.labs = c("Kaplan-Meier grupo1","Kaplan-Meier grupo2"))
# Prueba Log-rank
La prueba log-rank permite comparar las funciones de supervivencia de diferentes grupos: \[H_0:S_1(t)=S_2(t)=...=S_k(t)\] \[H_1:S_j(t_0)≠S_j'(t0) \text{ p.a. }j\neq j'\text{ y tiempo }t_0\]
## Call:
## survdiff(formula = surv ~ grupo, data = datos)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## grupo=1 21 9 19.3 5.46 16.8
## grupo=2 21 21 10.7 9.77 16.8
##
## Chisq= 16.8 on 1 degrees of freedom, p= 4e-055 Modelos paramétricos
Es posible trabajar con un modelo paramétrico cuando la distribución del tiempo hasta el evento se ajusta a una distribución conocida, \(T\sim f_T(t)\). Las más comunes son la exponencial, la log normal, la Gamma y la Weibull.
5.1 Exponencial
La función de densidad está dada por: \[f_T(t)=\lambda e^{-\lambda t} I_{(0,\infty)}(t)\] La función de supervivencia y la función de peligro (constante):
\[S(t)=P(T>t)=e^{-\lambda t},\qquad h(t)=\lambda\] Además: \[{\displaystyle E(T)={\frac {1}{\lambda }},\qquad V(T)={\frac {1}{\lambda ^{2}}}, \qquad Med(T)=\ln(2)\lambda}\] El parámetro de la distribución puede estimarse vía máxima verosimilitud.
5.1.1 Ejemplo
- Supongamos que \(T\sim exp(0.1)\) (anual), así \(S(5)=P(T>5)=1-F_t(5)\):
## [1] 0.6065307Es decir que bajo ese modelo probabilístico, la probabilidad de sobrevivir 5 años es del 60.6%.
- Si queremos estimar el parámetro correspondiente de la distribución con base en un conjunto de datos:
#Seguiremos trabajando con los datos del ejemplo anterior
#El primer paso sería determinar si los datos se ajustarían a una distribución exponencial
ks.test(datos$t,pexp)## Warning in ks.test(datos$t, pexp): ties should not be present for the
## Kolmogorov-Smirnov test##
## One-sample Kolmogorov-Smirnov test
##
## data: datos$t
## D = 0.86264, p-value < 2.2e-16
## alternative hypothesis: two-sided#en este caso trabajando con un 5% de significancia, los datos no tienen distribución exponencial, sin embargo seguiremos con el ajuste, solo a modo de ejemplo
library(flexsurv)
exp<-flexsurvreg(surv ~ 1, data = datos, dist = "exp")
exp## Call:
## flexsurvreg(formula = surv ~ 1, data = datos, dist = "exp")
##
## Estimates:
## est L95% U95% se
## rate 0.0555 0.0388 0.0793 0.0101
##
## N = 42, Events: 30, Censored: 12
## Total time at risk: 541
## Log-likelihood = -116.7667, df = 1
## AIC = 235.5333##
## time est lcl ucl
## 1 1 0.9460566 0.92542954 0.9626930
## 2 2 0.8950231 0.85641984 0.9267778
## 3 3 0.8467425 0.79255622 0.8922025
## 4 4 0.8010664 0.73345494 0.8589172
## 5 5 0.7578542 0.67876087 0.8268735
## 6 6 0.7169729 0.62814536 0.7960254
## 7 7 0.6782970 0.58130427 0.7663281
## 8 8 0.6417073 0.53795615 0.7377387
## 9 9 0.6070915 0.49784051 0.7102159
## 10 10 0.5743429 0.46071632 0.6837199
## 11 11 0.5433609 0.42636049 0.6582123
## 12 12 0.5140502 0.39456660 0.6336564
## 13 13 0.4863206 0.36514359 0.6100166
## 14 15 0.4352682 0.31271621 0.5653499
## 15 16 0.4117883 0.28939682 0.5442584
## 16 17 0.3895751 0.26781637 0.5239537
## 17 19 0.3486787 0.22936326 0.4855887
## 18 20 0.3298698 0.21225953 0.4674728
## 19 22 0.2952411 0.18178328 0.4332435
## 20 23 0.2793148 0.16822762 0.4170805
## 21 25 0.2499932 0.14407347 0.3865409
## 22 32 0.1695696 0.08375053 0.2962172
## 23 34 0.1517687 0.07172561 0.2745275
## 24 35 0.1435818 0.06637700 0.2642857#Gráfica de la función de supervivencia
plot(exp, type = "survival", ci = F,
conf.int = T, col = 1, col.obs = 2, xlab = "t", ylab = "S(t)")
#Para comparar los ajustes con los diferentes modelos probabilísticos:
#Distribuciones a comparar:
distr<-c("exp","weibull","llogis")
#función que calcula el ajuste con alguna distribución "f", para los datos de supervivencia "data"
comp<-sapply(distr,function(x) flexsurvreg(surv~1, dist = x), simplify = F)
comp## $exp
## Call:
## flexsurvreg(formula = surv ~ 1, dist = x)
##
## Estimates:
## est L95% U95% se
## rate 0.0555 0.0388 0.0793 0.0101
##
## N = 42, Events: 30, Censored: 12
## Total time at risk: 541
## Log-likelihood = -116.7667, df = 1
## AIC = 235.5333
##
##
## $weibull
## Call:
## flexsurvreg(formula = surv ~ 1, dist = x)
##
## Estimates:
## est L95% U95% se
## shape 1.141 0.850 1.532 0.171
## scale 17.904 13.082 24.502 2.866
##
## N = 42, Events: 30, Censored: 12
## Total time at risk: 541
## Log-likelihood = -116.4054, df = 2
## AIC = 236.8108
##
##
## $llogis
## Call:
## flexsurvreg(formula = surv ~ 1, dist = x)
##
## Estimates:
## est L95% U95% se
## shape 1.57 1.16 2.12 0.24
## scale 11.83 8.38 16.69 2.08
##
## N = 42, Events: 30, Censored: 12
## Total time at risk: 541
## Log-likelihood = -115.3511, df = 2
## AIC = 234.7022plot(comp[[1]], ci = F, conf.int = F, lty = 2,
xlab = "t", ylab = "S(t)")
plot(comp[[2]], ci = F, conf.int = F, lty = 3, add=T,col=3)
plot(comp[[3]], ci = F, conf.int = F, lty = 4, add=T,col=4)
legend("topright", c("KM", distr), lty = 1:(length(distr) + 1), col = 1:(length(distr) +
1))
6 Bibliografía
Kleinbaum, D. G., & Klein, M. (2010). Survival analysis (Vol. 3). New York: Springer.
Colosimo, E. A., & Giolo, S. (2006). Análise de Sobrevivência Aplicada. 1ª edição. São Paulo: Editora Edgard Blücher.