SurvivalCode
Vásquez Guerra Carlos Fernando
10/2020
Datos
Hay que considerar algunos puntos antes de iniciar el análisis de supervivencia
- En general, los datos que se proporcionarán para este tipo de análisis tendrán una variable que indicará la censura en sus observaciones. Aunque aveces se tendrá que crear esta en la limpieza de los datos.
Por ejemplo, en el inicio del libro se tienen las siguientes tablas respecto a los tiempos de infección (en meses) de pacientes en diálisis renal con diferentes procedimientos de cateterización. Los datos fueron tomados de (Klein and Moeschberger 2006)
| Colocación de catéter de forma quirúrgica | |||
|---|---|---|---|
| Tiempos de Infección: 1.5,3.5,4.5,4.5,5.5,8.5,8.5,9.5,10.5,11.5,15.5,16.5,18.5,23.5,26.5 | |||
| Observaciones Censuradas:2.5,2.5,3.5,3.5,3.5,4.5,5.5,6.5,6.5,7.5,7.5,7.5,7.5,8.5,9.5 |
| Colocación de catéter de forma percutánea | |||
|---|---|---|---|
| Tiempos de Infección: 0.5,0.5,0.5,0.5,0.5,0.5,2.5,2.5,3.5,6.5,15.5 | |||
| Observaciones Censuradas:0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,1.5,1.5,1.5,1.5,2.5 |
Dado esto, los datos en tal formato no nos ayudan para este tipo de análisis. Una manera útil para nosotros puede ser la siguiente
#SPC: Surgically_Placed_Catheter
SPC <- c(1.5,3.5,4.5,4.5,5.5,8.5,8.5,9.5,10.5,11.5,15.5,16.5,18.5,23.5,26.5,2.5,2.5,3.5,3.5,3.5,4.5,5.5,6.5,6.5,7.5,7.5,7.5,7.5,8.5,9.5)
# 1:censured, 0: not censured
SPC_censured <- c(rep(0, 15), rep(1, 15))
#PPC: Percutaneous Placed Catheter
PPC <- c(0.5,0.5,0.5,0.5,0.5,0.5,2.5,2.5,3.5,6.5,15.5, 0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,1.5,1.5,1.5,1.5,2.5)
# 1:censured, 0: not censured
PPC_censured <- c(rep(0, 11), rep(1, 15))Sólo para dar presentación vamos a cambiar los nombres de los vectores para indicar si la observación esta censurada o no y también los daremos en un orden aleatorio, ya que es común que las observaciones no tengan un orden ya que esto no es necesario.
set.seed(4)
names(SPC) <- SPC_censured
names(PPC) <- PPC_censured
sample(SPC) 1 0 1 0 0 0 1 1 1 1 0 0 0 1 1 1 0 0 1 0 1 0
6.5 15.5 3.5 4.5 8.5 16.5 7.5 5.5 4.5 9.5 8.5 3.5 1.5 8.5 7.5 2.5 18.5 5.5 6.5 23.5 3.5 26.5
0 1 0 1 0 1 1 0
9.5 3.5 11.5 2.5 4.5 7.5 7.5 10.5 sample(PPC) 0 1 1 1 1 1 1 0 0 0 0 1 0 1 0 1 0 0 1 1 1 1
0.5 0.5 1.5 1.5 0.5 0.5 0.5 0.5 3.5 6.5 0.5 0.5 2.5 0.5 0.5 1.5 0.5 15.5 0.5 2.5 1.5 0.5
1 0 0 1
0.5 0.5 2.5 0.5 Otra manera de tratar estos datos puede ser con un data frame con en la siguiente salida de código. La elección del formato dependerá del uso posterior, ya que algunas funciones destinadas a este análisis pueden requerir un indicador de censura en un vector.
tibble("Tiempos de infección" = SPC, "Censurados" = SPC_censured)tibble("Tiempos de infección" = PPC, "Censurados" = PPC_censured)- Siempre es ideal comenzar con un análisis descriptivo de los datos para considerar algunos puntos importantes que pueden tener repercusiones graves en un análisis posterior. Esto puede abarcar tanto gráficas como resúmenes estadísticos
Análisis de Supervivencia
A menos que el alumno se haya adelantado (👏), hasta el momento no sabemos como estimar la función de supervivencia, por lo que vamos a asumir que de inicio tenemos las siguientes funciones de supervivencia para los datos anteriores:
| Tiempo | Supervivencia |
|---|---|
| 1.5 | 1.0000000 |
| 2.5 | 0.9310345 |
| 3.5 | 0.8275862 |
| 4.5 | 0.7916042 |
| 5.5 | 0.7520240 |
| 6.5 | 0.6684658 |
| 7.5 | 0.5013493 |
| 8.5 | 0.4595702 |
| 9.5 | 0.4085069 |
Esta para el caso de los pacientes en diálisis renal donde el catéter se colocó de forma quirúrgica. Lo que en términos matemáticos sería lo siguiente:
\[ S(t) = \left\{ \begin{array}{ll} 1 & -\infty<t<2.5\\ 0.9310345 & 2.5\leq t<3.5\\ 0.8275862 & 3.5\leq t < 4.5\\ 0.7916042 & 4.5\leq t< 5.5\\ 0.7520240 & 5.5 \leq t< 6.5\\ 0.6684658 & 6.5 \leq t <7.5\\ 0.5013493 & 7.5 \leq t <8.5\\ 0.4595702 & 8.5 \leq t< 9.5\\ 0.4085069 & 9.5 \leq t \leq26.5\\ 0 & t> 26.5 \end{array} \right. \]
Para el caso de los pacientes en diálisis renal donde el catéter se colocó de forma percutánea se tienen los siguientes datos
Ahora, vamos a graficar las anteriores funciones de supervivencia
#Véase que vamos a agregar un punto al inicio y otro al final
SPCS <- tibble(Tiempo = c(0, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 26.5),
Surv = c(1,0.931, 0.828, 0.792, 0.752, 0.668, 0.501, 0.46, 0.409, 0))
PPCS <- tibble(Tiempo = c(0, 0.5, 1.5, 2.5, 3.5, 6.5, 15.5),
Surv = c(1, 0.6153846, 0.3692308, 0.3076923, 0.3076923, 0.3076923, 0))
plot_CS <- SPCS %>% ggplot(aes(x = Tiempo, y = Surv, color = "Quirúrgica")) +
geom_step() +
geom_step(data = PPCS, aes(color = "Percutánea")) +
labs(y = "S(t)") +
scale_color_manual(values = c("#848482", "#E23220")) +
general_theme
ggplotly(plot_CS) %>%
layout(legend = list(orientation = 'h', y = -0.2, title=list(text='Tipo de colocación del catéter: ')))¿Qué podemos concluir hasta ahora?
- ¿Existe alguna mejora con algún tipo de colocación del catéter? Claro, esto aún no lo podemos afirmar estadísticamente pero es un primer paso.
- Sólo con ver las gráficas, en qué tipo de colocación del catéter se obtiene un tiempo mediano de vida más grande?
- ¿En cuál procedimiento hay una mayor probabilidad de no obtener la infección en al menos 5 meses?
Aprovechando que tenemos la función de supervivencia, calculemos la función de riesgo. Recordemos que \(h(u_k) = 1-\frac{S(u_{k})}{S(u_{k-1})}\), por lo que se tienen los siguiente resultados para el caso cuando el catéter se colocó de forma quirúrgica
SPCS_risk <- SPCS %>% # Tomamos los datos
mutate(Surv_aux = lead(Surv)) %>% # Recorremos la supervivencia para facilitar los cálculos
mutate(quotien = Surv_aux/Surv) %>% # Hacemos el cociente
mutate(Risk_rate = 1-quotien) %>% #Obtenemos la tasa de fallo
head(-1) # El último valor no nos interesa ¿Por qué?
SPCS_risk %>% rename("$t$" = Tiempo, "$S(u_k)$" = Surv, "$S(u_{k-1})$" = Surv_aux,
"$S(u_k)/S(u_{k-1})$" = quotien, "$h(t)$" = Risk_rate) %>% #Para el diseño
style_wK() #Más diseño| \(t\) | \(S(u_k)\) | \(S(u_{k-1})\) | \(S(u_k)/S(u_{k-1})\) | \(h(t)\) |
|---|---|---|---|---|
| 0.0 | 1.000 | 0.931 | 0.9310000 | 0.0690000 |
| 2.5 | 0.931 | 0.828 | 0.8893663 | 0.1106337 |
| 3.5 | 0.828 | 0.792 | 0.9565217 | 0.0434783 |
| 4.5 | 0.792 | 0.752 | 0.9494949 | 0.0505051 |
| 5.5 | 0.752 | 0.668 | 0.8882979 | 0.1117021 |
| 6.5 | 0.668 | 0.501 | 0.7500000 | 0.2500000 |
| 7.5 | 0.501 | 0.460 | 0.9181637 | 0.0818363 |
| 8.5 | 0.460 | 0.409 | 0.8891304 | 0.1108696 |
| 9.5 | 0.409 | 0.000 | 0.0000000 | 1.0000000 |
Y de manera análoga para cuando el catéter se colocó de forma percutánea.
| \(t\) | \(S(u_k)\) | \(S(u_{k-1})\) | \(S(u_k)/S(u_{k-1})\) | \(h(t)\) |
|---|---|---|---|---|
| 0.0 | 1.0000000 | 0.6153846 | 0.6153846 | 0.3846154 |
| 0.5 | 0.6153846 | 0.3692308 | 0.6000001 | 0.3999999 |
| 1.5 | 0.3692308 | 0.3076923 | 0.8333332 | 0.1666668 |
| 2.5 | 0.3076923 | 0.3076923 | 1.0000000 | 0.0000000 |
| 3.5 | 0.3076923 | 0.3076923 | 1.0000000 | 0.0000000 |
| 6.5 | 0.3076923 | 0.0000000 | 0.0000000 | 1.0000000 |
Finalmente, obtener la función de riesgo acumulado es algo muy sencillo
SPCS_risk$Risk_rate %>%
cumsum() %>%
as_tibble() %>% # Una buena estructura de datos
mutate(Time = row_number()) %>% # Añadimos un índice para el tiempo
select(Time, value) %>% # Para un buen orden
rename("$t$" = Time, "$H(t)$" = value) %>%
style_wK()| \(t\) | \(H(t)\) |
|---|---|
| 1 | 0.0690000 |
| 2 | 0.1796337 |
| 3 | 0.2231120 |
| 4 | 0.2736170 |
| 5 | 0.3853192 |
| 6 | 0.6353192 |
| 7 | 0.7171555 |
| 8 | 0.8280251 |
| 9 | 1.8280251 |
Para el caso continuo se verá un ejemplo más adelante cuando se vea lo correspondiente a la sección de modelos paramétricos
Parámetros poblacionales
Como ya se menciono antes, es importante un análisis descriptivo de la información y en esta sección se ven varios de los puntos más importantes para esto. Tomando el ejemplo anterior (sólo el primer procedimiento), y por lo tanto los casos discretos, se tiene la siguiente información sobre la media residual, la esperanza y la varianza.
SPCS_summary <- summary_survival(SPCS$Tiempo, SPCS$Surv)
SPCS_summary$Residual_median %>%
rename("$t$" = Time, "$mr(t)$" = mr) %>%
style_wK()| \(t\) | \(mr(t)\) |
|---|---|
| 0.0 | 6.341000 |
| 2.5 | 6.810956 |
| 3.5 | 7.164251 |
| 4.5 | 6.909091 |
| 5.5 | 6.610372 |
| 6.5 | 6.441617 |
| 7.5 | 7.087824 |
| 8.5 | 5.997826 |
| 9.5 | 4.721271 |
| 26.5 | Inf |
- Esperanza: 6.341
- Var: 17.432719
Y respecto a los cuantiles, creo que la siguiente gráfica será de gran ayuda
library(gganimate)
library(scales)
library(gifski)
animate_quantiles <- data_filled(SPCS, .1) %>% add_row(Tiempo = seq(26.5, 30, .1), Surv = rep(0,36)) %>% #rellenamos los datos faltantes para hacer las transiciones
ggplot(aes(x = Tiempo, y = Surv, yend = Surv)) +
annotate("segment", x = SPCS$Tiempo, xend = c((SPCS$Tiempo %>% lead())[-10], 30),
y = SPCS$Surv, yend = SPCS$Surv) + #Agregamos de manera sintética la gráfica de supervivencia
annotate("segment", x = SPCS$Tiempo[-1], xend = SPCS$Tiempo[-1],
y = SPCS$Surv[-10], yend = (SPCS$Surv %>% lead())[-10]) +
general_theme +
geom_point(color = "#E23220") + #Puntos para la animación
geom_segment(xend = 31, color = "#89160A")+
lims(x = c(0,30))+
scale_x_continuous(breaks = seq(0,30, 5), labels = seq(0,30, 5))+
scale_y_continuous(breaks = seq(0, 1, .05), labels = seq(0, 1, .05),
sec.axis = sec_axis(~ 1- ., name = "Cuantil al",
breaks = seq(0, 1, .1),
labels = label_percent()))+
transition_time(Tiempo) +
ease_aes('linear') +
labs(title = 'Cuantil = {round(frame_time, 2)}', y = "S(t)")
final_animation<-animate(animate_quantiles,100, fps = 15, nframes = 150, duration = 30, res = 150, width = 864, height = 480, renderer = gifski_renderer())
anim_save("Plots/quantile_animate.gif",animation=final_animation)
Modelos Paramétricos
Será difícil encontrar aplicaciones que se ajusten perfectamente a un sólo modelo paramétrico, aunque no es imposible. Para este caso vamos a proponer una serie de problemas y una modelación con funciones paramétricas como las vistas en las sección 5 del Bookdown de referencia.
Vamos a suponer que en tu tiempo libre te dedicas a escribir un blog en internet o que llevas un poco de tiempo creando contenido para alguna plataforma de video. Te haz percatado que en promedio necesitas 10 días para que 500 personas vean tu trabajo, por lo que podríamos decir que 50 personas ven tu contenido por día. Entonces te haces la siguiente pregunta: ¿Cuál es la probabilidad de que en total 50 personas hayan visto mi contenido en menos de una semana? Para contestar esta pregunta veamos la siguiente gráfica.
Tl7 <- pexp(7, rate = 1/50) #Probabilidad acumulada a la semana
exp_example <- tibble(x= c(0,200)) %>%
ggplot(aes(x = x)) +
stat_function(fun = ~dexp(.x, rate = 1/50)) + #Gráficamos la distribución
annotate("segment", x = qexp(Tl7, rate = 1/50), #Agregamos linea
xend = qexp(Tl7, rate = 1/50),
y = 0, yend = dexp(7, rate = 1/50)) +
geom_area(stat = "function", fun = ~dexp(.x, rate = 1/50), #Agregamos relleno
fill = "red", xlim = c(0,7)) +
annotate("text", x = 37, y = dexp(7, rate = 1/50) + 0.001,
label = "$P(T\\leq0.7) = 0.1306418$") +
general_theme +
labs(x = "Tiempo", y = "Probabilidad")
#Función para dar un formato estilo Latex. Esta usa la librería tikzDevice
layoult_latex(exp_example, "Plots/exp_example.tex", height = 5, width = 8.35)
Sólo para aclarar algunas cosas, estamos considerando una distribución exponencial, es decir, \(T\sim exp(\lambda = 1/50)\), ya que utilizamos esta distribución en particular para modelar la cantidad de tiempo de espera para que suceda un evento. Por ejemplo la cantidad de tiempo para que un aparato electrónico falle, el tiempo de espera de llegada del siguiente autobús, etc. Respecto a \(\lambda\), recordemos que esta es una tasa que podemos interpretar como la cantidad de tiempo entre una cierta cantidad de eventos; en este caso tenemos que por una unidad de tiempo (días) suceden 50 eventos (50 visitas).
Dado este modelo, sabemos que \(\mathbb{E}[T] = \frac{1}{\lambda} = 50\) y \(Var[T] = \frac{1}{\lambda^2} = 2500 \implies \sigma = 50\); es decir que en promedio se tendrán 50 visitas las cuales pueden aumentar o disminuir en un 100%, lo cual hace dudar si este es un buen modelo. Además, sabemos que \(h(t) = \lambda = \frac{1}{50}\) en cada instante por lo que nuestra audiencia siempre sería la misma, lo cual no nos ayuda mucho ya que, en general, la probabilidad de que más personas vieran mi contenido aumentaría a medida que uno va ganando reconocimiento.
Una de las características del modelo exponencial es esa tasa de riesgo constante, por lo que un modelo de este estilo serviría para procesos donde se tenga un cierto control o que las condiciones externas no afecten demasiado a las probabilidades de que un evento suceda como los accidentes automovilísticos, accidentes de avión, las llegadas de un nuevo autobús con horarios estrictos, etc. En esta situación se esperaría que aumentará o disminuyera la probabilidad instantánea a medida que avanzan los días y no de manera constante, por lo que un modelo que se ajustaría mejor sería con una distribución Weibull.
Aprovechando este comportamiento, vamos a utilizar algunas funciones de aproximación numérica para calcular los parámetros poblacionales mediante la función de supervivencia, que en este caso tendríamos que \(S(t) = e^{-\frac{1}{50}t}\). Pero primero recordemos lo siguiente para el caso continuo:
- Media: \(\mu = \int_0^{\infty}S(t)dt\)
- Varianza: \(\sigma^2 = 2\int_0^\infty tS(t)dt-(\int_0^\infty S(t)dt)^2\)
- Media Residual: \(mr(x) = \frac{\int_x^\infty S(t)dt}{S(x)}\)
Media:
Surv_blog <- function(t) exp(-t/50)
#Compruebe las propiedades de la función de supervivencia con Surv_blog(0) y Surv_blog(Inf)
integrate(Surv_blog, lower = 0, upper = Inf)50 with absolute error < 0.00032Varianza:
#Primer sumando
FS <- 2*(integrate(f = function(t) t*Surv_blog(t), lower = 0, upper = Inf))$value
#Segundo sumando
SS <- (integrate(f = Surv_blog, lower = 0, upper = Inf))$value^2
#Varianza:
FS-SS[1] 2500Para la media residual vamos a utilizar un data frame y mostrar sólo algunos resultados.
#Ve que sucede si se agregan más datos
tibble(x = 1:10) %>%
#Obtenemos la integral de S(t) desde x para cada x
mutate(Numerator = map(x, ~(integrate(Surv_blog, lower = .x, upper = Inf))$value)) %>%
#Lo anterior regresa una lista por el comportamiento de la función integrate(), así que obtenemos de esas listas los valores que nos interesan
mutate(Numerator = map_dbl(Numerator, ~ `[`(.x, 1))) %>%
#Obtenemos los valores S(x) para cada x
mutate(Denominator = Surv_blog(x)) %>%
mutate(quotien = Numerator/Denominator) %>%
rename("$t$" = x, "$\\int_{t}^\\infty S(u)du$" = Numerator, "$S(t)$" = Denominator,
"$\\frac{\\int_{t}^\\infty S(u)du}{S(t)} = mr(t)$" = quotien) %>% #Para el diseño
style_wK()| \(t\) | \(\int_{t}^\infty S(u)du\) | \(S(t)\) | \(\frac{\int_{t}^\infty S(u)du}{S(t)} = mr(t)\) |
|---|---|---|---|
| 1 | 49.00993 | 0.9801987 | 50 |
| 2 | 48.03947 | 0.9607894 | 50 |
| 3 | 47.08823 | 0.9417645 | 50 |
| 4 | 46.15582 | 0.9231163 | 50 |
| 5 | 45.24187 | 0.9048374 | 50 |
| 6 | 44.34602 | 0.8869204 | 50 |
| 7 | 43.46791 | 0.8693582 | 50 |
| 8 | 42.60719 | 0.8521438 | 50 |
| 9 | 41.76351 | 0.8352702 | 50 |
| 10 | 40.93654 | 0.8187308 | 50 |
- ¿Por qué \(mr(x)\) es constante? Hint: ¿Cómo se define \(mr(x); \forall x>0\)?
Como nuestra variable de interés preferida en este curso es el tiempo, no se tratarán otras cuestiones que podrían resultar interesantes. Por ejemplo, tomando el caso anterior ¿Qué tipo de modelo convendría si ahora deseáramos conocer la información probabilística acerca de la cantidad de visitas entre un número determinado de días? 🤔
Ahora supongamos que se desea modelar el tiempo de vida de un aparato electrónico, como una televisión o las luces de un árbol de navidad. ¿Tiene sentido utilizar una distribución exponencial? ¿Qué significaría la propiedad de pérdida de memoria en este caso? 🤔
Las situaciones anteriores se pueden modelar con una distribución Weibull, así como otros problemas, como el tiempo que tarda en arrancar un automóvil, la vida de las mariposas Monarca y la vida de las estrellas; vamos a tomar este último para una estrella de las características del Sol. De acuerdo a (Wisniewski and Velasco Sotomayor 2001), el tiempo de vida en miles de millones de años de una estrella de este tipo sigue una distribución \(Weibull(\lambda = \frac{1}{11}, \gamma = \frac{8}{3})\). Vamos a ver las funciones de densidad y de riesgo de este comportamiento.
risk_weibull <- function(x, shape, scale) scale*shape*(scale*x)^(shape-1)
DF_sun <- tibble(x = c(0,.3)) %>%
ggplot(aes(x = x)) +
stat_function(fun = ~dweibull(.x, shape= 8/3, scale = 1/11), color = "#0F5148")+
guides(color = guide_legend(title=NULL))+
labs(x = "t", y = TeX('$f(t)$'))+
general_theme
HF_sun <- tibble(x = c(0,20)) %>%
ggplot(aes(x = x)) +
stat_function(fun = ~risk_weibull(.x, shape = 8/3, scale = 1/11), color = "#C13B06")+
guides(color = guide_legend(title=NULL))+
labs(x = "t", y = TeX('$h(t)$'))+
general_theme
DF_sun + HF_sun
Lo más interesante esta en la función de riesgo, la cual indica que a medida que pasa el tiempo de vida de una estrella como el sol, aumenta la probabilidad instantánea de morir aunque no de manera tan drástica, ya que recordemos que estamos hablando de miles de millones de años. Los parámetros poblacionales los obtendremos de manera análoga a lo anterior y se colocarán sólo resultados.
lambda_sun <- 1/11
shape_sun <- 8/3
Surv_sun <- function(t) exp(-(lambda_sun*t)^shape_sun)
PSun <- PParameters(Surv_sun, 100)- Media: 9.7780493
- Varianza: 15.5963184
PSun$mrVéase que en los primeros mil millones de años, la esperanza de vida residual es alta pero con cada unidad de tiempo que pasa, esta va decayendo drásticamente como se puede ver en las siguiente gráficas. La primera corresponde a la media vida residual y la segunda a la función de riesgo acumulado; en esta se observa como va aumentando la tasa de riesgo a través de los millones de años
Sun_plotMR <- PSun$mr %>%
ggplot(aes(x = t, y = rm, color = rm)) +
geom_point()+
scale_color_continuous(high = "#1BC603", low = "#AF1403") +
labs(x = NULL, y = "Media residual") +
guides(color = FALSE) +
general_theme
Sun_plotMR <- ggplotly(Sun_plotMR)
Sun_plotH <- tibble(x = seq(0,100, by = 0.1)) %>%
mutate(H = -log(Surv_sun(x))) %>%
ggplot(aes(x = x, y = H)) +
labs(x = NULL, y ="Riesgo acumulado") +
geom_line(color = "#AF1403")
Sun_plotH <- ggplotly(Sun_plotH)
subplot(Sun_plotMR, Sun_plotH,
titleY = TRUE, margin = 0.05) %>%
add_annotations(text = "Tiempo",
x = 0.5,
y = 0,
yref = "paper",
xref = "paper",
xanchor = "center",
yanchor = "bottom",
yshift = -35,
showarrow = FALSE,
font = list(size = 15)) %>%
layout(margin = list(b = 40))Finalmente, vamos a ver el siguiente ejemplo obtenido de (Wisniewski and Velasco Sotomayor 2001):
“Se sabe que en el aeropuerto de la ciudad de México un turista recién llegado corre un gran riesgo si toma un taxi en la calle, ya que muchos”taxistas" son realmente asaltantes. Para evitar tantos asaltos a los turistas en el aeropuerto, se ha instaurado un sistema de combis colectivas, las cuales dependen del mismo aeropuerto y son más seguras, aunque más caras. Las combis no tienen horario fijo de salida, pero salen tan proton llegan 6 pasajeros independientes (que van todos más o menos por la misma zona). Si los pasajeros independientes van llegando a un ritmo de = 1.5 pasajeros cada 5 minutos […]"
Dado esto, es evidente que la distribución exponencial juego un rol importante, aunque al considerar varios pasajeros para que suceda el evento de interés se necesita otro tipo de modelo. Así como la distribución binomial negativa determina la probabilidad de que suceda el k-ésimo éxito en varios ensayos que se distribuyen bernoulli, la distribución gamma determina la probabilidad del tiempo que se requiere para completar el k-ésimo éxito en ensayos tipo Poisson. Tampoco hay que olvidar que la suma de variables independientes con idéntica distribución exponencial se distribuye gamma.
Por lo que tendríamos el caso de una distribución \(gamma(\beta = 6,\lambda = \frac{3}{2})\) y los siguientes resultados:
Surv_combi <- function(t) 1-pgamma(q = t, shape = 6, rate = 3/2)
PCombi <- PParameters(Surv_combi, 35)- Media: 4
- Varianza: 2.6666667
PSun$mrJunto con las siguientes gráficas
risk_gamma <- function(x, shape, scale) dgamma(x = x, shape = 6, rate = 3/2)/Surv_combi(x)
DF_combi <- tibble(x = c(0,20)) %>%
ggplot(aes(x = x)) +
stat_function(fun = ~dgamma(.x, shape= 6, rate = 3/2), color = "#0F5148")+
guides(color = guide_legend(title=NULL))+
labs(x = "t", y = TeX('$f(t)$'))+
general_theme
SF_combi <- tibble(x = c(0,20)) %>%
ggplot(aes(x = x)) +
stat_function(fun = ~Surv_combi(.x), color = "#0F5148")+
guides(color = guide_legend(title=NULL))+
labs(x = "t", y = TeX('$S(t)$'))+
general_theme
HF_combi <- tibble(x = c(0,20)) %>%
ggplot(aes(x = x)) +
stat_function(fun = ~risk_gamma(.x), color = "#C13B06")+
guides(color = guide_legend(title=NULL))+
labs(x = "t", y = TeX('$h(t)$'))+
general_theme
DF_combi + SF_combi + HF_combi
Más allá de obtener la probabilidad específica del tiempo entre dos transportes, véase que la probabilidad de que el tiempo de llegada entre dos transportes sea mayor a 10 es muy baja, lo cual es un buen indicador para el usuario que desee usar este servicio. Además de que la función de riesgo nos reafirma el hecho de que, dado que el usuario ya ha esperado un cierto tiempo, la probabilidad de llegada del transporte aumenta considerablemente. En general, el usuario tendrá que esperar 4 minutos que, en promedio, este puede aumentar 1.6 minutos o si tiene suerte disminuirlos.
Como bien se menciona en el capítulo 3 del Bookdown de referencia, la función de riesgo nos puede ayudar a tener una primera idea sobre el modelo paramétrico al que podemos ajustar nuestros datos. Hasta ahora con los ejemplos anteriores tenemos que las funciones de riesgo son constantes o monótonas crecientes o decrecientes, pero no se han tratado aquellas que pueden cambiar su tendencia a medida que avanza el tiempo; para este caso véase las siguientes funciones de riesgo para los modelos Log-logístico y Log-normal respectivamente.
En las anteriores gráficas se tienen comportamientos crecientes y decrecientes, aunque no hay que olvidar que se pueden obtener comportamientos monótonos. Ambas distribuciones tienen una gran diversidad de aplicaciones como se mencionan en los artículos Log-logistic distribution y Log-normal distribution donde también se da una “justificación” al porque elegir estas distribuciones. Para el caso de la log-normal se pueden encontrar diversas aplicaciones como los tiempos de resolución de un cubo de rubik en una competición realizada en el 2019 o el tiempo que permanecen los usuarios en linea viendo diversos tipos de artículos. Respecto a la distribución log-logistica se pueden encontrar aplicaciones en hidrología o en el tiempo que transcurre en datos que son recibidos por alguna aplicación desde un ordenador desde la salida hasta su procesamiento en otros ordenadores.
Verosimilitud C&T
Uno de los temas más importantes y utilizados para la estimación paramétrica es la función de verosimilitud y falta ver lo correspondiente a esta con datos censurados y truncados. Para estos fines no se repetirá la teoría vista en el capítulo 6 del Bookdown y se utilizarán los resultados que se plasman en él y en caso contrario se harán las aclaraciones pertinentes.
Con el fin de dar continuidad a los ejemplos anteriores, vamos a retomarlos y calcular los distintos parámetros por medio del estimador máximo verosimil. Comencemos con el ejemplo referente al tiempo de visitas en un blog de internet. Como en la sección anterior sólo se trataron las distribuciones con sus valores exactos, vamos a simular un conjunto de datos junto con un indicador de datos censurados en distintos tipos comenzando con el caso de censura por la derecha tipo 1.
#No olvidemos que es importante la reproducibilidad
set.seed(31)
Tblog_simulated <- rexp(100, rate = 1/50)
#Para la asignación de valores censurados, vamos a crear una variable aleatoria que, con 30% de probabilidad, indique que la observación fue censurada
right_censured_Tb <- sample(c(0,1), size = 100, prob = c(.70, .30), replace = T)Como bien sabemos, la función de verosimilitud en este caso es \(\mathscr{L} = \prod^{n}_{i=1}f(t_{i})^{\delta_{i}}S(t_{i})^{1-\delta_{i}}\) y en particular para la función exponencial es \(\mathscr{L} = \lambda^{r}e^{-\lambda\sum_{i= 1}^{n}t_{i}}\) donde \(r\) es el número de observaciones no censuradas. Por lo que el estimador máximo verosimil para \(\lambda\) es \(\hat{\lambda} = \frac{r}{\sum_{i= 1}^{n}t_{i}}\). Vamos a calcularlo de manera exacta y después con una aproximación numérica.
#Calculamos r
r <- sum(right_censured_Tb == 0)
hat_lambda <- r / sum(Tblog_simulated)
cat(paste("Parámetro real: 1/50 = 0.02", "\nParámetro verosimil: ", hat_lambda))Parámetro real: 1/50 = 0.02
Parámetro verosimil: 0.0107660140471416¿Qué pasará si los reducimos la cantidad de valores censurados?
#Vamos a hacer una función que vaya los distintos valores estimados de lambda para una nueva variable aleatoria con distintas probabilidades para los datos censurados.
sample_censured_blog <- function(prob){
right_censured_Tb <- sample(c(0,1), size = 100, prob = c(1-prob, prob), replace = T)
sum(right_censured_Tb == 0) / sum(Tblog_simulated)
}
#Vamos a hacer una gráfica para evitar tener una mejor representación de estos datos.
tibble(p_censured = seq(.30, 0, by = -.01)) %>%
mutate(lambda = map_dbl(p_censured, ~sample_censured_blog(.x))) %>%
ggplot(aes(x = p_censured, y = lambda)) +
geom_line(color = "royalblue3") +
general_theme +
labs(x = "Probabilidad para los datos censurados simulados", y = TeX("$\\hat{\\lambda}$")) +
scale_x_reverse()
Hay que considerar que estamos trabajando con datos pseudoaleatorios, por lo que todos estos resultados pueden variar en diferentes condiciones, pero sin importar este punto, es clara una tendencia hacia el estimador máximo verosimil ¿Porqué piensas que sucede este comportamiento al reducir la cantidad de valores censurados? Pon a prueba esto corriendo varias veces el anterior código.
Ahora vamos a calcularlos utilizando la paquetería fitdistrplus. Este paquete esta diseñado para el ajuste de distintas distribuciones en un conjunto de datos dados. Como la sección trata con valores censurados o truncados, sólo se verán las funciones que nos permitan hacer un ajuste en este tipo de estructuras y no se abarcarán todas las funciones interesantes que puede tener este paquete.
La función que tomaremos de este paquete es fitdistrplus::fitdistcens(), la cual ajusta una distribución univariada a datos censurados dando los estimadores máximo verosimil. Por detrás, esta función utiliza algunas funciones de optimización como stats::optim con diferentes métodos numéricos; por default fitdistcens() utiliza los métodos Nelder-Mead y BFGS con el objetivo de evitar problemas con las derivadas de distintas funciones.
Para este caso de datos censurados, la función anterior solicita que los datos tengan una estructura especial. Los datos deben estar contenidos en un data frame con dos columnas:
leftNApara datos censurados por la izquierda.- El valor del límite izquierdo en el intervalo para datos censurados por intervalo.
- El valor de las observaciones no censuradas.
rightNApara datos censurados por la derecha.- El valor del límite derecho en el intervalo para datos censurados por intervalo.
- El valor de las observaciones no censuradas.
¿Basta con estas dos columnas para englobar los distintos tipos de censura? Sí. Veamos un ejemplo donde están cada uno de estos casos, supongamos que el estudio tiene como objetivo medir el tiempo de espera en minutos en una recepción y estos fueron medidos desde que abrió la oficina hasta 32 minutos después de que abrió.
#Exammple data for fitdistcens
EDFF <- data.frame(left = c(NA, 32, 22, 30), right = c(0, NA, 24, 30))
row.names(EDFF) <- c("Censura por la izquierda", "Censurad por la derecha", "Censura por intervalos", "No censura")
EDFF %>% style_wK()| left | right | |
|---|---|---|
| Censura por la izquierda | NA | 0 |
| Censurad por la derecha | 32 | NA |
| Censura por intervalos | 22 | 24 |
| No censura | 30 | 30 |
Entonces, para nuestro caso tendríamos lo siguiente:
(data_blog_censored_right <- tibble(left = Tblog_simulated) %>%
mutate(right = if_else(right_censured_Tb == 0, left, NA_real_)))Y ya podemos utilizar la función fitdistcens
library(fitdistrplus)
# La función solicita explicitamente un data frame
fit_blog <- fitdistcens(as.data.frame(data_blog_censored_right), "exp")
fit_blog$estimate rate
0.0107664 Véase que los resultados numéricos y teóricos son muy similares.
Antes de seguir, esta es una gran oportunidad para introducir una de las librerías más utilizadas en este tipo de análisis, es decir el paquete survival. Con esta librería se podrán realizar la mayoría de cálculos de manera automatizada que se verán después en este curso. Por el momento basta con saber que nos permite una interfaz para modelar adecuadamente los datos censurados con un objeto de tipo supervivencia, esto con la función survival::Surv(). Se puede obtener información detallada en la documentación oficial pero, de manera resumida, se tienen los siguientes parámetros para esta función.
time: Tiempo para los datos censurados por la derecha o el tiempo inicial en el caso de censura por intervalos.time2: Tiempo final para el caso de censura por intervalos: (start, end]event: Parámetro para agregar un indicador de censura:- 0 = vivo, 1 = muerto | TRUE = vivo, FALSE = muerto | 1 = vivo, 2 = muerto.
- Para los demás casos:
- 0 = censurado por la derecha.
- 1 = no censura.
- 2 = censurado por la izquierda.
- 3 = censurado por intervalos.
type: Indicador del tipo de censura.
Si se utiliza la función sin colocar el nombre explicito del argumento, se tienen dos casos:
Surv(obj1, obj2):time= obj1 yevent= obj2.Surv(obj1, obj2, obj3):time= obj1,time2= obj2 yevent= obj3.
Entonces, ajustando un modelo con la función fitdistcens utilizando los resultados de un objeto de supervivencia, se tiene lo siguiente para nuestro caso, lo cual es completamente equivalente a lo anterior, sólo que ahora lo hacemos con un objeto diseñado para el análisis de supervivencia.
library(survival)
surv_blog <- Surv(Tblog_simulated, right_censured_Tb)
surv_blog %>% head(15) [1] 2.215791 14.709631+ 15.670172+ 88.394989 8.425897+ 33.023225+
[7] 15.186753+ 211.934514+ 134.970640 38.388905 76.331577+ 6.660786+
[13] 58.328007 36.347459+ 78.324413+#so : from survival object
data_blog_censored_right_fso <- tibble(left = surv_blog[,"time"]) %>%
mutate(right = if_else(surv_blog[,"status"] == 0, left, NA_real_))
fit_blog_fso <- fitdistcens(as.data.frame(data_blog_censored_right_fso), "exp")
fit_blog_fso$estimate rate
0.0107664 Hasta ahorita, en las funciones que se han mencionado, se engloba la censura por la derecha en el tipo 1, ya que, como se vio en el Bookdown, los demás tipos producen funciones de verosimilitud proporcionales a la del tipo 1. Por ejemplo, para la censura tipo 2 y en caso exponencial, desarrollando todas las expresiones necesarias, se tendría que \(\hat{\lambda} = \frac{r}{\left(\sum\limits_{i = 1}^nt_i\right)-1}\). En nuestro ejemplo, suponiendo que se desean observar las primeras 70 veces que 50 personas ven el contenido del blog, tendríamos resultados bastante similares:
# Véase que por el comportamiento del estimador no es necesario simular más datos censurados
r <- 70
hat_lambdaII <- r / (sum(Tblog_simulated)-1)
cat(paste("Parámetro real: 1/50 = 0.02",
"\nParámetro verosimil 1: ", hat_lambda,
"\nParámetro verosimil 2: ", hat_lambdaII))Parámetro real: 1/50 = 0.02
Parámetro verosimil 1: 0.0107660140471416
Parámetro verosimil 2: 0.0115960896455989Como ya vimos, los resultados los podemos obtener de manera más mecanizada con ciertas funciones, así que vamos a ver los ejemplos restantes con los métodos numéricos.
En el segundo ejemplo visto en la sección Modelos Paramétricos se tenia que este es modelado con una distribución \(Weibull(\lambda = \frac{1}{11}, \gamma = \frac{8}{3})\). Vamos a suponer que durante el estudio se descubrió una nueva técnica para calcular de manera más precisa la edad de una estrella, por lo que algunas observaciones quedaron censuradas por la izquierda; entonces se decide calcular los nuevos estimadores mediante el método de máxima verosimilitud.
set.seed(32)
Tsun_simulated <- rweibull(100, scale = 1/11, shape = 8/3)
#Para la asignación de valores censurados, vamos a crear una variable aleatoria que, con 20% de probabilidad, indique que la observación fue censurada
left_censured_Ts <- sample(c(0,1), size = 100, prob = c(.80, .20), replace = T)
surv_sun <- Surv(Tsun_simulated, left_censured_Ts, type = "left")
data_sun_censored <- tibble(right = surv_sun[,"time"]) %>%
mutate(left = if_else(surv_sun[,"status"] == 0, right, NA_real_))
fit_sun <- fitdistcens(as.data.frame(data_sun_censored), "weibull")
fit_sun$estimate shape scale
2.3502016 0.0840386 Los nuevos parámetros establecen una distribución \(Weibull(\hat{\lambda} = 0.0840386, \hat{\gamma} = 2.3502016)\) la cual es muy similar a la distribución sin datos censurados \(Weibull(\lambda = 0.09090909, \gamma = 2.666667)\) como se puede ver en la siguiente gráfica
tibble(x = c(0,.3)) %>%
ggplot(aes(x = x)) +
stat_function(fun = ~dweibull(.x, shape= 8/3, scale = 1/11), color = "#0F5148")+
stat_function(fun = ~dweibull(.x, shape= fit_sun$estimate["shape"],
scale = fit_sun$estimate["scale"]), color = "#C13B06")+
guides(color = guide_legend(title=NULL))+
labs(x = "t", y = TeX('$f(t)$'))+
general_theme
Finalmente, para el caso de censura por intervalos tomaremos nuestro último ejemplo de la sección Modelos paramétricos en el cual se trata un nuevo sistema de transporte. En este caso vamos a suponer que cada 2 minutos se registra el tiempo de llegada de un autobús en base a nuestro modelo \(gamma(\beta = 6, \lambda = \frac{3}{2})\).
set.seed(33)
Tcombi_simulated <- rgamma(100, shape = 6, rate = 3/2)
#Para este caso no vamos a generar ninguna variable aleatoria para agregar censura, si no que vamos a categorizar el vector Tcombi_simulated.
combi_censured <- tibble(time = Tcombi_simulated,
intervals = findInterval(time, seq(0, 10, by = 2))) %>%
#Necesitamos que por cada observación se tengan los límites del intervalo al que pertenece
mutate(lower = map_dbl(intervals, ~switch(.x, 0, 2, 4, 6, 8)),
upper = map_dbl(intervals, ~switch(.x, 2, 4, 6, 8, 10)))
surv_combi <- Surv(time = combi_censured$lower,
time2 = combi_censured$upper, type = "interval2")
surv_combi %>% head(10) [1] [2, 4] [6, 8] [2, 4] [2, 4] [4, 6] [2, 4] [4, 6] [2, 4] [2, 4] [2, 4]data_combi_censored <- tibble(right = surv_combi[,"time2"]) %>%
mutate(left = surv_combi[,"time1"])
fit_combi <- fitdistcens(as.data.frame(data_combi_censored), "gamma")
fit_combi$estimate shape rate
7.197411 1.828527 Para este caso tenemos las siguientes gráficas. Tómese en cuenta que en este ejemplo todas las observaciones son consideradas como censuradas.
tibble(x = c(0,15)) %>%
ggplot(aes(x = x)) +
stat_function(fun = ~dgamma(.x, shape = 6, rate = 3/2), color = "#0F5148")+
stat_function(fun = ~dgamma(.x, shape = fit_combi$estimate["shape"],
rate = fit_combi$estimate["rate"]), color = "#C13B06")+
guides(color = guide_legend(title=NULL))+
labs(x = "t", y = TeX('$f(t)$'))+
general_theme
Modelos No Paramétricos
Hasta el momento se ha trabajado con R de una manera muy tradicional, es decir que hemos tomado ejemplos que se ajustan de una manera muy adecuada al análisis que hasta ahora se ha abarcado. Esto con el único fin de mostrar una opción, utilizando herramientas computacionales, para calcular ciertos resultados teóricos. Ejemplos como los que se vieron en las últimas dos secciones serán difíciles de encontrar o pueden llevar cierto tiempo para determinar si cumplen las condiciones necesarias para realizar un modelo paramétrico.
Como ya se vio en el la sección correspondiente al Bookdown, estamos considerando dos métodos no paramétricos para el ajuste de la función de supervivencia; el método actuarial y el método producto-límite o Kaplan-Meier.
Ya que las operaciones necesarias para crear los siguientes resultados son fáciles de codificar, vamos a utilizar un conjunto de librerías que nos ayuden a realizar todo el cálculo necesario sin necesidad de mostrar aquí una codificación personalizada. Aunque se recomienda y pide al lector crear por ellos mismos las funciones para obtener los resultados que aquí se van a presentar. Comencemos primero con los ejemplos que podemos encontrar en el Bookdown para después ver un ejemplo de un conjunto de datos en la siguiente entrega.
Método Actuarial
Para este caso, vamos a utilizar el paquete KMsurv, el cual en su mayoría es un paquete que contiene una gran cantidad de data sets para el análisis de supervivencia; estos datos, en especial, son aquellos que se utilizan en el libro (Klein and Moeschberger 2006). Finalmente, se tiene la función KMsurv::lifetab() la cual crea una tabla de vida de acuerdo a una cohorte, es decir de acuerdo a un grupo de personas que han experimentado un cierto acontecimiento en un mismo intervalo de tiempo.
Esta función contiene 4 parámetros:
tis(time intervals): Un vector que indica la configuración de los intervalos de tiempo.ninit(number of initial subjects): Número de individuos del estudio.nlost(number of lost subjetcs): Vector que indica la cantidad de individuos “perdidos” por alguna razón del estudio. Es decir, los censuradosnevent(number of events): Vector que indica el número de individuos que experimentaron el evento.
Esta función regresará un data frame con toda la información necesaria para el análisis de supervivencia, entre ellas las mismas columnas que se tienen en el ejemplo del Bookdown y otras que hasta el momento no se han explorado, por lo que estas serán omitidas hasta que sea el momento de estudiarlas. Recordemos que tenemos los siguientes datos
Así que ahora sólo vamos a utilizar la función antes mencionada y darle nuestros datos, pero primero vamos a organizar nuestros datos para que sea más sencillo todo el proceso.
library(KMsurv)
intervals <- c(0, 12, 24, 36, 48, 60, NA)
n_subjects <- dim(myeloma_data)[1]
#myeloma_data_organized es un data frame donde se resume la cantidad de censurados y de muertes por intervalo de los datos de myeloma.
censured_myeloma <- myeloma_data_organized$censured
event_myeloma <- myeloma_data_organized$death
lifetab(tis = intervals, ninit = n_subjects, nlost = censured_myeloma, nevent = event_myeloma) %>%
dplyr::select(nevent, nlost, nsubs, nrisk,surv)Véase que la columna que indica la supervivencia comienza con 1 y en el ejemplo dado en el Bookdown el primer valor de la supervivencia es distinto de 1. Esto es porque esta función da los valores de supervivencia en el inicio del intervalo (se puede consultar el libro citado para verificar dicha información), por lo que basta con agregar un tiempo extra para ver el comportamiento final del último intervalo.
#Vamos a agregar el último valor de tiempo que aparece en los datos
intervals <- c(0, 12, 24, 36, 48, 60, max(myeloma_data$Survival_time), NA)
#Vamos a agregar ceros ya que no hay datos en el último intervalo de la tabla.
censured_myeloma <- c(myeloma_data_organized$censured,0)
event_myeloma <- c(myeloma_data_organized$death,0)
(myeloma_surv_actuar <- lifetab(tis = intervals, ninit = n_subjects,
nlost = censured_myeloma, nevent = event_myeloma) %>%
dplyr::select(nevent, nlost, nsubs, nrisk,surv))Y finalmente sólo falta nuestro formato
myeloma_surv_actuar %>%
rename("$d_j$" = "nevent", "$c_j$" = "nlost", "$n_j$" = "nsubs",
"$n^{*}_j$" = "nrisk", "$S(t)$" = "surv") %>%
style_wK()| \(d_j\) | \(c_j\) | \(n_j\) | \(n^{*}_j\) | \(S(t)\) | |
|---|---|---|---|---|---|
| 0-12 | 16 | 4 | 48 | 46.0 | 1.0000000 |
| 12-24 | 10 | 4 | 28 | 26.0 | 0.6521739 |
| 24-36 | 1 | 0 | 14 | 14.0 | 0.4013378 |
| 36-48 | 3 | 1 | 13 | 12.5 | 0.3726708 |
| 48-60 | 2 | 2 | 9 | 8.0 | 0.2832298 |
| 60-91 | 4 | 1 | 5 | 4.5 | 0.2124224 |
| 91-NA | 0 | 0 | 0 | 0.0 | 0.0236025 |
Ya que la gráfica correspondiente se puede encontrar en el Bookdown, esta se omitirá.
Kaplan-Meier
Bien, para comparar después los métodos, vamos a tomar nuevamente este ejemplo. Para este caso vamos a de la función survival::Surv() que se vio anteriormente y una nueva de la misma paquetería: survival::survfit(), la cual crea una curva de supervivencia estimada mediante el método de Kaplan-Meier. Como estamos considerando sólo observaciones censuradas por la derecha y de tipo I, podemos utilizar la manera más sencilla de la función Surv()
myeloma_surv <- Surv(myeloma_data$Survival_time, myeloma_data$Status)
myeloma_surv_fit <- survfit(myeloma_surv~1)
typeof(myeloma_surv_fit)[1] "list"Véase que lo que regresa la función survfit es una lista con distintos objetos importantes que utilizaremos para crear la tabla que se presenta en el Bookdown.
names(myeloma_surv_fit) [1] "n" "time" "n.risk" "n.event" "n.censor" "surv"
[7] "std.err" "cumhaz" "std.chaz" "type" "logse" "conf.int"
[13] "conf.type" "lower" "upper" "call" Además de esto, el objeto myeloma_surv_fit es un objeto de tipo survfit, por lo que se recomienda probar diferentes códigos como myeloma_surv_fit, summary(myeloma_surv_fit) y plot(myeloma_surv_fit). Otra cosa que es importante aclarar es que la función solicita en sus argumentos una fórmula, las cuales son muy comunes en R y la manera de expresarlas es usando el símbolo ~. Al escribir survfit(myeloma_surv~1) la función entiende que no se desea hacer alguna segregación de los datos, lo cual puede ser de mucha utilidad; más abajo se dará un ejemplo de esto pero antes véase la construcción de la tabla correspondiente. Por esta ocasión ignoraré el formato LaTex para mostrar todos los resultados.
tibble(time = myeloma_surv_fit$time,
dk = myeloma_surv_fit$n.event,
nk = myeloma_surv_fit$n.risk,
ck = myeloma_surv_fit$n.censor,
"S(t)" = myeloma_surv_fit$surv)Ahora veamos la gráfica que se genera con esta curva de supervivencia comparada con la dada en el método actuarial.
#Vamos a modificar un poco la base myeloma_surv_actuar para agregar una variable de tiempo y un renglón más
myeloma_surv_actuarM <- myeloma_surv_actuar %>%
mutate(time = c(na.omit(intervals))) %>%
add_row(time = c(max(myeloma_data$Survival_time),100), surv = 0)
tibble(time = myeloma_surv_fit$time,
S = myeloma_surv_fit$surv) %>%
add_row(time = c(0, 100), S = c(1,0)) %>% #Agregamos dos renglones: 1 para el tiempo 0 y otro por diseño
ggplot(aes(x = time, y = S)) +
geom_step(aes(color = "K-M")) +
geom_step(data = myeloma_surv_actuarM, aes(y = surv, color = "Actuar")) +
labs(x = "Tiempo", y = TeX("$\\hat{S}(t)$"), colour = "Tipo de estimación: ") +
scale_color_manual(values = c( "#c45d16", "#086dbf")) +
theme(legend.position = "top", legend.direction = "horizontal") +
general_theme
¿Qué pasará si modificamos un poco los intervalos?
¿Que podemos comentar sobre estas dos gráficas?
Finalmente, vamos segregar los datos y utilizar el estimador \(K-M\) para comparar las gráficas de supervivencia. Consideremos al sexo del individuo como un estrato para determinar su supervivencia.
myeloma_surv_sex <- survfit(Surv(myeloma_data$Survival_time, myeloma_data$Status)~Sex, data = myeloma_data)
male_strata <- myeloma_surv_sex$strata[1]
female_strata <- myeloma_surv_sex$strata[2]
tibble(time = myeloma_surv_sex$time,
S = myeloma_surv_sex$surv) %>%
mutate(sex = factor(c(rep(1, male_strata), rep(2, female_strata)))) %>%
add_row(time = c(0, 100, 0, 100), S = c(1,0, 1, 0), sex = factor(c(1, 1, 2, 2))) %>% #Agregamos dos por tener dos grupos
ggplot(aes(x = time, y = S, color = sex)) +
geom_step() +
labs(x = "Tiempo", y = TeX("$\\hat{S}(t)$"), colour = "Sexo: ") +
scale_color_manual(values = c( "royalblue2", "mediumpurple1"), labels = c("Hombres", "Mujeres")) +
theme(legend.position = "top", legend.direction = "horizontal") +
general_theme
Habría que formalizarlo con una prueba, pero hay que notar que, por las gráficas, en este tipo de cancer en los primeros 25 meses de tener la enfermedad las probabilidades son bastante similares tanto en hombres como mujeres, pero después de haber pasado ese lapso, las probabilidades de sobrevivir a dicho mal es menor en las mujeres en comparación de los hombres. Aunque habría que considerar que en el estudio se tuvieron 21 hombres y 14 mujeres, por lo que, o bien, en general esta enfermedad, en las condiciones que se dio el estudio, se tiene una estimación más precisa sobre la supervivencia en hombres que en mujeres, lo cual podría producir un sesgo importante que haría creer que la supervivencia en mujeres es más baja.
Estimaciones NP
!NO¡, NP significa No Paramétricas más no estimaciones sobre problemas no computables. Bueno, antes de ver más ejemplos en los que podemos aplicar un análisis de supervivencia, vamos a ver algunos resultados que devuelven las funciones anteriores que no se llegaron a ver debido a que la esta sección no se había estudiado. Con respecto a la función lifetab() no hay mucho que abarcar ya que esta devuelve las estimaciones para la función de riesgo, su desviación estándar al igual que para la de supervivencia. Por lo que vamos a completar esta función para obtener al menos lo mostrado en el Bookdown.
#Vamos a suponer un alpha = 0.05 (es decir, lo usual)
cuantil975 <- qnorm(.975)
(new_myeloma_surv_actuar <- myeloma_surv_actuar %>%
mutate(qj = nevent/nrisk, #Necesitamos las qj para el cálculo de la varianza estimada
varS = surv^2*cumsum(qj/((1-qj)*nrisk))) %>%
mutate(Confidence_5 = surv-cuantil975*sqrt(varS),
Confidence_95 = surv+cuantil975*sqrt(varS)) %>% #Vamos a arreglar los valores por si se tienen valores superiores a 1 o a 0 en el intervalo
mutate(Confidence_5 = if_else(is.na(Confidence_5)|Confidence_5<0, 0, Confidence_5),
Confidence_95 = if_else(Confidence_95>1, 1, Confidence_95)) %>% #Es distinto en este caso
mutate(Confidence_95 = if_else(is.na(Confidence_95), 0, Confidence_95)))Ahora vamos a obtener una mejor gráfica; las líneas con el color degradado representan los intervalos de confianza.
#Gráfica de los datos de myeloma con el método actuarial con intervalos de confianza
myelomaPlot_actuarWC <- new_myeloma_surv_actuar %>%
mutate(time = c(na.omit(intervals))) %>%
add_row(time = c(max(myeloma_data$Survival_time),100), surv = 0, #Añadimos algunos valores extra para el diseño
Confidence_5 = 0, Confidence_95 = 0) %>%
dplyr::select(time, surv, Confidence_5, Confidence_95) %>% #Siempre será lo ideal trabajar con sólo lo necesario
ggplot(aes(x = time, y = surv)) +
geom_step(color = "#c45d16") +
geom_step(aes(y = Confidence_5), color = "#c45d16", alpha = 0.5) +
geom_step(aes(y = Confidence_95), color = "#c45d16", alpha = 0.5) +
general_theme +
labs(x = "Tiempo", y = "Supervivencia estimada")
ggplotly(myelomaPlot_actuarWC)Ahora vamos a ver todo lo que nos falto sobre la función survfit()
myeloma_surv_summary <- summary(survfit(Surv(myeloma_data$Survival_time, myeloma_data$Status)~1))
tibble(time = myeloma_surv_summary$time,
dk = myeloma_surv_summary$n.event,
nk = myeloma_surv_summary$n.risk,
ck = myeloma_surv_summary$n.censor,
S = myeloma_surv_summary$surv,
varS = myeloma_surv_summary$std.err^2) %>%
mutate("lower 95% CI" = myeloma_surv_summary$lower) %>%
mutate("lower to hand" = S-qnorm(0.975)*sqrt(varS)) %>% #Hacemos a mano los cálculos para el lower 95% CI
mutate("uper 95% CI" = myeloma_surv_summary$upper) %>%
mutate("upper to hand" = S+qnorm(0.975)*sqrt(varS)) %>% #Hacemos a mano los cálculos para el upper 95% CI
rename("S(t)" = S) #Más diseño¿Porqué se calcularon los intervalos de manera “manual”? Si el lector es observador, se habrá dado cuenta de que los intervalos obtenidos por la función survfit() y a mano son diferentes. Esto es debido a que es común que aquí los límites de los intervalos se salgan del rango \([0,1]\), por lo que se proponen distintas transformaciones (log y log-log) para corregir esto; en el caso por default, los intervalos son calculados con la transformación logarítmica, para obtener los intervalos de acuerdo al estimador de Greenwood se debe cambiar el argumento conf.type con "plain". Véase en la siguiente gráfica los tres tipos de intervalos que podemos tener con esta función.
myeloma_conf_int_Grewnwood <- summary(survfit(Surv(myeloma_data$Survival_time, myeloma_data$Status)~1,
conf.type = "plain"))
myeloma_conf_int_log_log <- summary(survfit(Surv(myeloma_data$Survival_time, myeloma_data$Status)~1,
conf.type = "log-log"))
tibble(time = myeloma_surv_summary$time,
S = myeloma_surv_summary$surv,
CIlog_lower = myeloma_surv_summary$lower,
CIlog_upper = myeloma_surv_summary$upper,
CIplain_lower = myeloma_conf_int_Grewnwood$lower,
CIplain_upper = myeloma_conf_int_Grewnwood$upper,
CIlog_log_lower = myeloma_conf_int_log_log$lower,
CIlog_log_upper = myeloma_conf_int_log_log$upper) %>%
ggplot(aes(x = time, y = S)) +
geom_step(color = "#086dbf") +
geom_step(aes(y = CIlog_lower, colour = "1"), alpha = 0.5) +
geom_step(aes(y = CIlog_upper, colour = "1"), alpha = 0.5) +
geom_step(aes(y = CIplain_lower, colour = "2"), alpha = 0.5) +
geom_step(aes(y = CIplain_upper, colour = "2"), alpha = 0.5) +
geom_step(aes(y = CIlog_log_lower, colour = "3"), alpha = 0.5) +
geom_step(aes(y = CIlog_log_upper, colour = "3"), alpha = 0.5) +
labs(x = "Tiempo", y = "Supervivencia estimada", colour = "Tipo de intervalos de confianza: ") +
scale_color_manual(labels = c("1" = "log", "2" = "plain", "3" = "log-log"),
values = c("1" = "firebrick1", "2" = "darkorange3", "3" = "chocolate4")) +
general_theme +
theme(legend.position = "bottom")
Ahora, con respecto a la función de riesgo acumulado tenemos los siguientes resultados proporcionados por la función survfit(), los cuales puede comprobarse fácilmente que corresponden al estimador Nelson-Aalen. Además de estos resultados obtenidos de manera automática, se agrega el cálculo para los intervalos de confianza.
tibble(dk = myeloma_surv_summary$n.event,
nk = myeloma_surv_summary$n.risk,
H = myeloma_surv_summary$cumhaz,
stdH = myeloma_surv_summary$std.chaz) %>%
mutate(upper = H + cuantil975*stdH/H,
lower = H - cuantil975*stdH/H)Y mejor vamos a ver todo en una gráfica.
tibble(t = myeloma_surv_summary$time,
H = myeloma_surv_summary$cumhaz,
stdH = myeloma_surv_summary$std.chaz) %>%
mutate(upper = H + cuantil975*stdH/H,
lower = H - cuantil975*stdH/H) %>%
ggplot(aes(x = t, y = H)) +
geom_line(color = "firebrick2") +
geom_line(aes(y = upper), color = "firebrick2", alpha = 0.5) +
geom_line(aes(y = lower), color = "firebrick2", alpha = 0.5) +
general_theme +
labs(x = "Tiempo", y = "Función de Riesgo acumulado")
Ahora, vamos a ver que tan diferentes son las estimaciones de la función de supervivencia; con el método K-M y mediante el estimador Fleming-Harrington
tibble(time = myeloma_surv_summary$time,
S = myeloma_surv_summary$surv,
H = myeloma_surv_summary$cumhaz,
S_FH = exp(-H)) %>%
ggplot(aes(x = time, y = S, color = "1")) +
geom_step() +
geom_step(aes(y = S_FH, color = "2")) +
scale_color_manual(labels = c("1" = "K-M", "2" = "F-H"),
values = c("1" = "#086dbf", "2" = "firebrick2")) +
labs(x = "Tiempo", y = "Función de Superviviencia Estimada", colour = "Tipo de estimación: ") +
general_theme
¿Crees que hay alguna diferencia significativa?
Ahora vamos a ver una propuesta para el código correspondiente a las bandas de confianza. Véase como estas se comparan entre sí y entre los intervalos de confianza.
library(km.ci)
fit_myeloma <- survfit(Surv(myeloma_data$Survival_time, myeloma_data$Status)~1, conf.type = "plain")
EP <- summary(km.ci::km.ci(fit_myeloma, conf.level = 0.95, method = "epband"))
HW <- summary(km.ci::km.ci(fit_myeloma, conf.level = 0.95, method = "hall-wellner"))
tibble(time = myeloma_surv_summary$time,
S = myeloma_surv_summary$surv,
CIplain_lower = myeloma_conf_int_Grewnwood$lower,
CIplain_upper = myeloma_conf_int_Grewnwood$upper) %>%
ggplot(aes(x = time, y = S)) +
geom_step(color = "black") +
geom_step(aes(y = CIplain_lower, colour = "1"), alpha = 0.5) +
geom_step(aes(y = CIplain_upper, colour = "1"), alpha = 0.5) +
geom_step(aes(y = c(EP$lower,NaN), colour = "2"), alpha = 0.5) +
geom_step(aes(y = c(EP$upper, NaN), colour = "2"), alpha = 0.5) +
geom_step(aes(y = c(HW$lower, NaN), colour = "3"), alpha = 0.5) +
geom_step(aes(y = c(HW$upper, NaN), colour = "3"), alpha = 0.5) +
labs(x = "Tiempo", y = "Supervivencia estimada", colour = "Tipo de intervalos/bandas de confianza: ") +
scale_color_manual(labels = c("1" = "plain", "2" = "EPBand", "3" = "Hall.Wellner"),
values = c("1" = "#086dbf", "2" = "navyblue", "3" = "firebrick4")) +
general_theme +
theme(legend.position = "bottom")
¿Qué puedes concluir sobre esto?
Finalmente, vamos a ver un poco de la librería survminer() que puede ser atractiva para algunos
library(survminer)
fit_myeloma <- surv_fit(Surv(Survival_time, Status) ~ 1, data = myeloma_data)
plots_myeloma_survminer <- list()
plots_myeloma_survminer[[1]] <- ggsurvplot(fit_myeloma, data = myeloma_data, ggtheme = general_theme)
plots_myeloma_survminer[[2]] <- ggsurvplot(fit_myeloma, data = myeloma_data, fun = "event", ggtheme = general_theme)
plots_myeloma_survminer[[3]] <- ggsurvplot(fit_myeloma, data = myeloma_data, fun = "cumhaz", ggtheme = general_theme)
arrange_ggsurvplots(plots_myeloma_survminer, nrow = 1, ncol = 3, print = TRUE)
Y también podemos utilizar un estrato para dividir los resultados, obtener una tabla de riesgo y otras propiedades rápidas útiles.
fit_myeloma <- surv_fit(Surv(Survival_time, Status) ~ Sex, data = myeloma_data)
ggsurvplot(fit_myeloma, data = myeloma_data,
ggtheme = general_theme,
risk.table = T,
conf.int = T,
pval = T,
palette = c("darkorchid3", "dodgerblue4"),
legend = "bottom",
legend.title = "Sexo: ",
legend.labs = c("Hombres", "Mujeres"))
Finalmente, se dejan aquí un par de enlaces que pueden ser de gran utilidad al momento de ajustar algún modelo paramétrico en un análisis de supervivencia.
- Examples of focused model comparison: parametric survival models
- flexsurv: A Platform for Parametric Survival Modelling in R
Pruebas de Hipótesis
Preparación de los datos
Ahora sí, para esta sección vamos a tomar otro ejemplo basados en una base de datos que podemos encontrar en diferentes fuentes. Para este caso vamos a tomar los famosos datos de Billboard donde se recopila la información de las canciones más populares de acuerdo a su año y semana. Para este caso vamos a tomar sólo los datos que corresponden a todo el año 2000 junto con unas pocas semanas más; puedes obtener estos datos aquí.
# Primero cargamos los datos
(billboard <- tbl_df(read.csv("Data/billboard.csv", stringsAsFactors = FALSE)))Se puede ver fácilmente que la base necesita mucha limpieza. Este ejemplo aparece en el artículo Tidy Data de Hadley Wickham y aunque en el no se ve el código para dicha limpieza, aquí se proporciona esperando que el entendimiento del código sea el menor de los problemas.
billboard_tidy <- billboard %>%
gather(week, rank, wk1:wk76,na.rm=TRUE) %>% #Pasamos de un formato amplio a largo (algunas columnas crearon una nueva variable)
mutate( #Damos un mejor formato a las fechas
week = extract_numeric(week),
date = as.Date(date.entered) + 7 * (week - 1)) %>%
dplyr::select(-date.entered) #No necesitamos esto
billboard_tidy %>% arrange(artist, track) #Sólo por presentación, ordenamosSe puede ver con calma que ya no se tienen datos perdidos y que tenemos un buen formato para trabajar. Vamos a ver cuantas semanas se tienen en esta base
unique(billboard_tidy$week) [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
[26] 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
[51] 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65¡Genial! Tenemos una gran cantidad de información a través del tiempo. Veamos rápidamente cuales son esas canciones que duraron más de 60 semanas.
billboard_tidy %>% filter(week >= 60)Que interesante es todo esto, en este link puedes ver un Dashboard muy interesante con esta información. Ahora vamos a pensar en como podemos usar esta base para un análisis de supervivencia. Vamos a definir nuestro estudio como lo siguiente: Se desea estudiar cuantas semanas puede permanecer en la lista una canción que este entre un determinado rango de rankeo. Vamos a suponer que los rangos son de 10 en 10 (El rankeo más pequeño es 1). Teniendo este análisis podríamos responder preguntas del siguiente estilo:
- Si una canción no baja de un cierto rankeo, ¿Cuál es probabilidad de seguir en la lista?
- ¿En que grupo de rankeo puede estar más segura una canción para no salir de la lista?
- ¿Cuál es el peor grupo de rankeo en el que unca canción puede permanecer?
Entonces, vamos a preparar la base
#Función para colocar el grupo de rankeo
billboard_tidy %>%
mutate(group_rank = floor(rank/10)) %>%
dplyr::select(-date, -time, -year, -rank)Hasta aquí todo bien pero la variable de interés es la semana en que se sale de la lista (no en el momento en que se cambia de rankeo). Así que vamos a determinar la última semana de permanencia de cada canción, para esto debemos encontrar dicha semana. Si una canción pasa a otro grupo es una observación censurada más no “muerta”.
billboard_tidy %>%
mutate(group_rank = floor(rank/10)) %>%
dplyr::select(-date, -time, -year, -rank) %>%
group_by(track) %>%
mutate(week_max = max(week)) %>% #Obtenemos la última semana en que esta la canción
mutate(max_group = if_else(week == week_max, 1, 0)) %>% #Obtenemos en que grupo se quedo la canción antes de salir
group_by(group_rank, track) %>% # max_group no da un valor por grupo como deseamos, así que agrupamos de mejor manera
mutate(status = sum(max_group)) #Ahora si se tiene el último grupo al que perteneció la canción antes de salir de la lista y esta será nuestra variable para indicar el eventoEntonces ya tenemos nuestra variable que indicaría cuando una canción salio y además por grupo de rankeo. Por lo que sólo falta aplicar el estimador K-M. Vamos a hacerlo de una manera más adecuada a nuestro nivel de programación.
(models_billboard <- billboard_tidy %>%
mutate(group_rank = floor(rank/10)) %>%
dplyr::select(-date, -time, -year, -rank) %>%
group_by(track) %>%
mutate(week_max = max(week)) %>% #Obtenemos la última semana en que esta la canción
mutate(max_group = if_else(week == week_max, 1, 0)) %>% #Obtenemos en que grupo se quedo la canción antes de salir
group_by(group_rank, track) %>% # max_group no da un valor por grupo como deseamos, así que agrupamos de mejor manera
mutate(status = sum(max_group)) %>% #Ahora si se tiene el último grupo al que perteneció la canción antes de salir de la lista y esta será nuestra variable para indicar el evento
group_by(group_rank) %>% #Ahora solo agrupamos por el grupo de rankeo.
dplyr::select(-week_max) %>% #Ya no necesitamos esta variable
nest() %>% #Anidamos en data frames independientes por grupos para hacer más rápido el análisis
mutate(surv_model = map(
data, function(df) survfit(Surv(week, status)~1, data = df)
)))Sólo para aclarar, la variable surv_model contiene en cada uno de sus elementos un modelo de supervivencia basado en el estimador K-M hecho con la función survfit() y Surv(). Por ejemplo, vamos a ver el correspondiente al grupo más alto de rankeo.
survBillboar_1rank <- (models_billboard %>% arrange(group_rank))$surv_model[[11]]
#La siguiente función simplemente condensa el código previamente visto para crear este tipo de graficas de acuerdo a un objeto de supervivencia.
surv_graph(survBillboar_1rank, intervals = TRUE, color = "#086dbf")
Ahora sería interesante ver como se comportan cada una de las curvas de supervivencia
survBillboar_plots <- surv_graph(survBillboar_1rank, intervals = FALSE, colour = TRUE, color = "a")
#Agregamos múltiples capas con diferentes curvas de supervivencia.
for(i in 2:dim(models_billboard)[1]){
survBillboar_plots <- survBillboar_plots +
surv_layer(model = (models_billboard %>% arrange(group_rank))$surv_model[[i]],
intervals = FALSE, colour = TRUE, color = letters[i])
}
#Vamos a usar una paquetería divertida:
library(pjocolors)
colors_Billboar <- c(pjocolors::pjo_palettes$BatOfTheLabyrinth, pjo_palettes$TitansCurse)
survBillboar_plots +
scale_colour_manual(labels = 0:10,
values = c("black", colors_Billboar)) +
labs(colour = "Grupo de rankeo: ") +
theme(legend.position = "bottom") +
guides(colour = guide_legend(nrow = 1))
¿Qué podemos concluir hasta ahora con esta gráfica?
Si uno es observador, habrá notado que grupos 0, 1 y 2 no tienen una gran varianza, esto es debido a que, por el comportamiento de los datos, en esos grupos, mientras una canción permanezcan en dicha categoría, esta no saldrá de la lista, aunque eso no garantiza que la canción no saldrá de la lista. Vamos a retirar aquellos grupos que nos parezcan sospechosos a simple vista por el efecto que no tiene considerar la salida de un grupo de rankeo.
survBillboar_plots <- surv_graph((models_billboard %>% arrange(group_rank))$surv_model[[4]],
intervals = FALSE, colour = TRUE, color = "a")
#Agregamos múltiples capas con diferentes curvas de supervivencia.
for(i in c(5:11)){
survBillboar_plots <- survBillboar_plots +
surv_layer(model = (models_billboard %>% arrange(group_rank))$surv_model[[i]],
intervals = FALSE, colour = TRUE, color = letters[i])
}
survBillboar_plots +
scale_colour_manual(labels = 3:10,
values = colors_Billboar[3:10]) +
labs(colour = "Grupo de rankeo: ") +
theme(legend.position = "bottom") +
guides(colour = guide_legend(nrow = 1))
- ¿Entre que grupos de rankeo existe una diferencia más grande de tiempo?
- ¿Es indistinto estar en los grupos 6, 7 y 8?
- ¿Se ajustará algún modelo paramétrico a estas distribuciones estimadas?
1 población
Bien, ahora que ya tenemos un análisis sobre la base de datos ya aplicamos un ajuste K-M por estratos, vamos a hacer algunas pruebas de hipótesis. Para nuestra primer prueba de 1 población, vamos a tratar de ver si alguna de las funciones estimadas con las que contamos siguen alguna distribución, ya sea exponencial, weibull o gamma.
# Creamos una función que, de acuerdo a los datos que tenemos de billboar, cada uno de los subsets en data tengan un formato correcto la función fitdistcens
to_fitdistcens <- function(data){
tibble(left = data$week) %>%
mutate(right = if_else(data$status == 0, left, NA_real_))
}
(parModels_billboard <- models_billboard %>%
mutate(to_fit = map(data, function(x) to_fitdistcens(x))) %>%
mutate(exp_fit = map(to_fit, function(x) fitdistcens(as.data.frame(x), "exp"))) %>%
mutate(weibull_fit = map(to_fit, function(x) fitdistcens(as.data.frame(x), "weibull"))) %>%
mutate(gamma_fit = map(to_fit, function(x) fitdistcens(as.data.frame(x), "gamma"))))Si se lee el vingette Frequently Asked Questions del paquete fitdistrplus, se tiene el siguiente mensaje si se desea usar la función gofstat() en los anteriores objetos fitdistcens.
Function gofstat is not yet proposed in our package for fits on censored data but to develop one is among one of our objectives in the future. Published works on goodness-of-fit statistics based on the empirical distribution function for censored data generally focused on data containing only one type of censoring (e.g. right censored data in survival data). Build such statistics in the general case, with data containing in the same time (right, left and interval censoring), remains tricky. Nevertheless, it is possible for any type of censored data, to use information criteria (AIC and BIC given in the summary of an object of class fitdistcens) to compare the fits of various distributions to a same data set.
Por lo que vamos a obtener el modelo con el AIC más pequeño.
#Una función para obtener el modelo con el AIC más pequeño
aic_smaller_model <- function(exp, weibull, gamma){
smallerAIC <- min(exp$aic, weibull$aic, gamma$aic)
if(exp$aic == smallerAIC){
exp
}else{
if(weibull$aic){
weibull
}else{
gamma
}
}
}
(parModels_billboard <- parModels_billboard %>% ungroup() %>%
mutate(betterModel = pmap(list(exp_fit, weibull_fit, gamma_fit), function(exp, weibull, gamma) aic_smaller_model(exp, weibull, gamma))))Veamos el primer caso
cat("distribución ", parModels_billboard$betterModel[[1]]$distname, "con parámetros: ", parModels_billboard$betterModel[[1]]$estimate)distribución weibull con parámetros: 1.156293 10.25227Por lo que la función de riesgo para esta distribución es
\[ h(t) = \lambda\gamma(\lambda t)^{\gamma-1} = 0.09753939*1.156293(0.09753939t)^{1.156293-1} = 0.1127841(0.09753939\cdot t)^{0.156293} \]
Véase que si este hubiera resultado en un modelo exponencial, los cálculos se simplifican mucho, ya que la tasa de riesgo es \(\lambda\) y la esperanza de esta distribución es \(1/\lambda\) por lo que es sencillo obtener los esperados y la varianza de \(Z\), pero en este caso tendremos que calcular estos dos valores. Vamos a realizar nuestra prueba:
parameters <- as.vector(parModels_billboard$betterModel[[1]]$estimate)
observados <- sum(parModels_billboard$surv_model[[1]]$n.event/parModels_billboard$surv_model[[1]]$n.risk)
esperados <- sum(sapply(parModels_billboard$surv_model[[1]]$time, function(x) 0.1127841*(0.09753939*x)^0.156293))
z <- observados-esperados
var_z <- sum(sapply(parModels_billboard$surv_model[[1]]$time,
function(x) (0.1127841*(0.09753939*x)^0.156293))/parModels_billboard$surv_model[[1]]$n.risk)
statistic <- z^2/var_z
cat("Estadístico T=", statistic)Estadístico T= 0.8439997Finalmente vamos a agregar el cuantil teórico con el cual comparamos a nuestro estadístico y el correspondiente \(p-value\).
cat("Estadístico T =", statistic,
"\nCuantil al 95%: ", qchisq(.95, df = 1),
"\nP-value: ", 1- pchisq(statistic, df = 1))Estadístico T = 0.8439997
Cuantil al 95%: 3.841459
P-value: 0.3582553¿Qué podemos concluir entonces?
Vamos a hacer un poco más general todo. En el siguiente chunk se encuentra una función que aplicará la prueba de manera automática a un modelo de tipo fitdistcens dado.
general_1_sample <- function(fit_model, surv_model, alpha = 0.05){
observados <- sum(surv_model$n.event/surv_model$n.risk)
z <- 0
var_z <- 0
parameters <- as.vector(fit_model$estimate)
if(fit_model$distname == "weibull"){
esperados <- sum(sapply(surv_model$time,
function(x) (1/parameters[2])*parameters[1]*((1/parameters[2])*x)^(parameters[1]-1)))
z <- observados-esperados
var_z <- sum(sapply(surv_model$time,
function(x) ((1/parameters[2])*parameters[1]*
((1/parameters[2])*x)^(parameters[1]-1)))/surv_model$n.risk)
}else{
if(fit_model$distname == "exp"){
#Ayuda mucho a que la tasa de riesgo sea constante
esperados <- length(surv_model$time)*parameters[1]
z <- observados-esperados
var_z <- sum(parameters[1]/surv_model$n.risk)
}else{ #Gamma
esperados <- sum(sapply(surv_model$time,
function(x) (dgamma(.x, parameters[1],
rate = parameters[2])/(1-pgamma(.x, parameters[1],
rate = parameters[2])))))
z <- observados-esperados
var_z <- sum(sapply(surv_model$time,
function(x) (dgamma(.x, parameters[1],
rate = parameters[2])/(1-pgamma(.x, parameters[1],
rate = parameters[2]))))/surv_model$n.risk)
}
}
statistic <- z^2/var_z
quantil <- qchisq(1-alpha, df = 1)
list(model = fit_model$distname,
parameters = fit_model$estimate,
statistic = statistic,
result = if_else(statistic>quantil, "RejectH0", "CannotRejectH0"),
quantil = quantil,
pvalue = 1- pchisq(statistic, df = 1))
}Y aquí lo probamos con el ejemplo anterior.
general_1_sample(parModels_billboard$betterModel[[1]],surv_model= parModels_billboard$surv_model[[1]], 0.05)$model
[1] "weibull"
$parameters
shape scale
1.156293 10.252268
$statistic
[1] 0.8439991
$result
[1] "CannotRejectH0"
$quantil
[1] 3.841459
$pvalue
[1] 0.3582555Entonces, vamos a aplicarlo a todos los “mejores modelos” y vamos a obtener información interesante.
parModels_billboard %>%
#Aplicamos la función a cada uno de los mejores modelos
mutate(test = map2(betterModel, surv_model, ~general_1_sample(.x, .y, 0.05))) %>%
#Obtenemos si se rechazo o no la prueba y el estadístico
mutate(dist = map_chr(test, ~'[['(.x, "model")),
stat = map_dbl(test, ~'[['(.x, "statistic")),
reject = map_chr(test, ~'[['(.x, "result"))) %>%
#Ordenamos por
arrange(stat) %>%
#Solo tomamos las variables que nos interesan
dplyr::select(group_rank, dist, stat, reject)Véase que sólo para los grupos 4 y 8, el mejor modelo basado en la estadística \(AIC\) se “ajusta” a los modelos propuestos, en este caso weibull. Y para el modelo que menos se ajusta a esta propuesta es el modelo exponencial en el grupo de rankeo 10.
Log-Rank
Ahora, hagamos lo correspondiente para la prueba Log-Rank para dos curvas de supervivencia. Por ejemplo las curvas correspondientes a los grupos de rankeo 8 y 9
#Vamos a ordenarlos
parModels_billboard <- parModels_billboard %>% arrange(group_rank)
#Hay grupo de rankeo 0
#Obtenemos los modelos 8 y 9
S1 <- parModels_billboard$surv_model[[9]]
S2 <- parModels_billboard$surv_model[[10]]
#Calculamos los fallos totales al tiempo ti
di <- S1$n.event+S2$n.event
#Calculamos los individuos en riesgo al tiempo ti
ni <- S1$n.risk+S2$n.risk
#Calculamos los esperados y Var(di) de acuerdo al modelo hipergeométrico
e1i <- S1$n.risk*di/ni
vard1 <- S1$n.risk*S2$n.risk*di*(ni-di)/((ni^2)*(ni-1))
#Calculamos el estadístico
L <- (sum(S1$n.event-e1i))^2/(sum(vard1))
cat("Estadístico: ", L, "\ncuantil al 95%: ", qchisq(0.95, df = 1), "\np-value: ", 1 - pchisq(L, df = 1))Estadístico: 194.4549
cuantil al 95%: 3.841459
p-value: 0Donde recordemos que se rechaza la hipótesis nula cuando \(L^2>J_{1-\alpha}\). Por lo que hay una muy fuerte evidencia estadística para rechazar que \(S_8(t) = S_9(t)\), lo cual es evidente al ver las gráficas de estas dos poblaciones.
surv_graph((parModels_billboard$surv_model)[[9]],
intervals = FALSE, colour = TRUE, color = "a") +
surv_layer(model = (parModels_billboard$surv_model)[[10]],
intervals = FALSE, colour = TRUE, color = "b") +
scale_colour_manual(labels = 8:9, values = colors_Billboar[8:9]) +
labs(colour = "Grupo de rankeo: ") +
theme(legend.position = "bottom") +
guides(colour = guide_legend(nrow = 1))
Como se ve, es sencillo realizar estas pruebas de hipótesis, por lo que es factible pensar en generalizar esta prueba como la anterior, aplicar esta entre las distintas combinaciones de modelos y determinar cuales son los pares de modelos que son aprobadas por la prueba, es decir aquellas funciones de supervivencia que son estadísticamente iguales o diferentes. Para que este contenido no sea aburrido, se deja como ejercicio al lector. Hay que tener cuidado con esta implementación, ya que se tendrían \(n \choose 2\), en este caso se tendrían ¡55!
Wilcoxon
Para la prueba de Wilcoxon tenemos una implementación bastante similar (solo se agrega la ponderación), de hecho, sólo se agrega la implementación.
#Calculamos el estadístico
W <- (sum(ni*(S1$n.event-e1i)))^2/(sum(ni^2*vard1))
cat("Estadístico: ", W, "\ncuantil al 95%: ", qchisq(0.95, df = 1), "\np-value: ", 1 - pchisq(L, df = 1))Estadístico: 166.9611
cuantil al 95%: 3.841459
p-value: 0Y para agregar un poco más de contenido, aquí se deja un data frame con las pruebas aplicadas a todas las combinaciones de 2 en 2 de todos los modelos.
Y finalmente, vamos a ver dichas graficas.
surv_graph((parModels_billboard$surv_model)[[1]],
intervals = FALSE, colour = TRUE, color = "a") +
surv_layer(model = (parModels_billboard$surv_model)[[2]],
intervals = FALSE, colour = TRUE, color = "b") +
surv_layer(model = (parModels_billboard$surv_model)[[3]],
intervals = FALSE, colour = TRUE, color = "c") +
surv_layer(model = (parModels_billboard$surv_model)[[4]],
intervals = FALSE, colour = TRUE, color = "d") +
surv_layer(model = (parModels_billboard$surv_model)[[6]],
intervals = FALSE, colour = TRUE, color = "e") +
surv_layer(model = (parModels_billboard$surv_model)[[7]],
intervals = FALSE, colour = TRUE, color = "f") +
surv_layer(model = (parModels_billboard$surv_model)[[8]],
intervals = FALSE, colour = TRUE, color = "g") +
surv_layer(model = (parModels_billboard$surv_model)[[9]],
intervals = FALSE, colour = TRUE, color = "h") +
scale_colour_manual(labels = c(0,1,2,3,5,6,7,8), values = colors_Billboar[c(1,2,3,4,6,7,8,9)]) +
labs(colour = "Grupo de rankeo: ") +
theme(legend.position = "bottom") +
guides(colour = guide_legend(nrow = 1))
Proportional Hazar Models
Finalmente, sólo resta tratar la sección de modelos proporcionales como una propuesta de un modelo lineal sobre distintas covariables y las variables de interés con las que ya hemos tratado. Para este caso, la función que estaremos utilizando es survival::coxph() la cual, entre todos sus argumentos, nos interesan los siguientes:
formula. Ya hemos tratado con este tipo de expresiones, por ejemplo en la funciónsurvfit:survfit(Surv(time, Status)~strata.method. En este parámetro se coloca el nombre del tipo de estimador que se utilizará para los distintos parámetros del modelo de Cox. Las opciones sonexactpara utilizar los estimadores derivados del modelo de máxima verosimilitud.efronpara el modelo de verosimilitud parcial. Este es el valor por defecto en el parámetromethod.breslowpara el uso del estimador de Breslow.
Para nuestro ejemplo anterior será un poco complicado utilizar las covariables restantes ya que estas tienen una gran diversidad entre sus valores, por lo que vamos a utilizar otros datos para este análisis, estos nos serán proporcionados por el paquete JM y serán los datos pbc2, los cuales contienen el seguimiento de 312 pacientes aleatorizados con cirrosis biliar primaria, una enfermedad hepática autoinmune rara, en Mayo Clinic. Toda la información al respecto de las variables se puede encontrar en la correspondiente documentación.
library(JM)
(pbc <- JM::pbc2)Rápidamente, se puede observar que las variables que nos interesan son status2, la cual indica si un paciente murió (1) o si sobrevivío o tuvo un trasplante (0), además de years, que es el número de años entre el registro y el tiempo de fallecimiento, trasplante o análisis del estudio, lo que ocurra primero.
Entonces, véase que sucede al aplicar un modelo de Cox basado en las covariables drug, sex y age, es decir, sobre el siguiente modelo:
\[ h(t|X(t)) = h_0 \cdot exp \{ \beta_1 X_1 +\beta_2X_2(t) \} = h_0(t) \cdot exp \{ \beta_1\cdot drug+\beta_2\cdot sex +\beta_3\cdot age\} \]
model_cox_1 <- coxph(Surv(years, status2) ~ drug + sex+ age, data = pbc)
summary(model_cox_1)Call:
coxph(formula = Surv(years, status2) ~ drug + sex + age, data = pbc)
n= 1945, number of events= 725
coef exp(coef) se(coef) z Pr(>|z|)
drugD-penicil -0.130201 0.877919 0.075958 -1.714 0.0865 .
sexfemale -0.547998 0.578106 0.093038 -5.890 3.86e-09 ***
age 0.042168 1.043070 0.003966 10.634 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
drugD-penicil 0.8779 1.1391 0.7565 1.0188
sexfemale 0.5781 1.7298 0.4817 0.6937
age 1.0431 0.9587 1.0350 1.0512
Concordance= 0.602 (se = 0.012 )
Likelihood ratio test= 151 on 3 df, p=<2e-16
Wald test = 157.8 on 3 df, p=<2e-16
Score (logrank) test = 157.8 on 3 df, p=<2e-16Esta función nos devuelve mucha información, por ejemplo, tenemos directamente los valores de \(e^{\beta_i}\) para así tener el siguiente modelo:
\[ h(t|X(t)) = h_0(t) \cdot exp \{-0.130201\cdot drug-0.547998\cdot sex +0.042168\cdot age\} \]
También tenemos directamente las pruebas que son de nuestro interés como la significancia de los parámetros a través de la prueba de Wald. En este caso tenemos que se rechaza la hipótesis de que para algún coeficiente \(\theta_i = \theta_{H_0}\), por lo que nuestros coeficientes sí otorgan información relevante al modelo; además se tiene de manera individual los \(p-values\) para cada uno de los parámetros como si de una regresión lineal se tratara. Si somos estrictos, todos los coeficientes son significativos, aunque la variable drug podría ser reconsiderada. Además de esta, tenemos los resultados para las siguientes pruebas
- Concordancia.
- Prueba de puntajes.
- De acuerdo al modelo que se proponga, se da el estadístico \(R^2\).
- Likehood ratio test. Este último es uno de los más conocidos y utilizados en la toma de decisiones así como el \(BIC\) y \(AIC\).
Otro punto interesante que hay que notar es que la variable sex (la cual es un factor con niveles male y female) fue codificada como sexfemale dando a entender que esta variable se codifico de manea binaria, tal como se muestra en nuestro libro de referencia. Algo similar sucedió con la variable drugD. En caso de que uno quiera crear variables dummies, se puede utilizar la función I(), la cual es ampliamente utilizada en modelos y fórmulas en R. Aquí se deja una publicación con un gran resumen de funciones basadas en supervivencia junto con un ejemplo de esto.
Véase que la variable numérica age no fue codificada como las anteriores, pero esta fue tomada de manera literal, bien se pudo haber utilizado alguna transformación como un logaritmo para reducir su varianza o alguna otra función lineal o no lineal; por ejemplo, podemos colocar log(age) en lugar de la edad o hacer un filtro:
model_cox_2 <- coxph(Surv(years, status2) ~ drug + sex + I(age>40), data = pbc)
summary(model_cox_2)Call:
coxph(formula = Surv(years, status2) ~ drug + sex + I(age > 40),
data = pbc)
n= 1945, number of events= 725
coef exp(coef) se(coef) z Pr(>|z|)
drugD-penicil -0.005658 0.994358 0.074380 -0.076 0.939366
sexfemale -0.559153 0.571693 0.093553 -5.977 2.27e-09 ***
I(age > 40)TRUE 0.391450 1.479123 0.116668 3.355 0.000793 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
drugD-penicil 0.9944 1.0057 0.8595 1.1504
sexfemale 0.5717 1.7492 0.4759 0.6867
I(age > 40)TRUE 1.4791 0.6761 1.1768 1.8591
Concordance= 0.556 (se = 0.011 )
Likelihood ratio test= 50.32 on 3 df, p=7e-11
Wald test = 53.17 on 3 df, p=2e-11
Score (logrank) test = 54.94 on 3 df, p=7e-12Es interesante ver que parace que necesitamos de las observaciones con edades menores a 40 años para que la variable drug no se considere como no significativa. Aquí también utilizaremos la función I() para evitar problemas y realizar libremente estas transformaciones.
model_cox_3 <- coxph(Surv(years, status2) ~ drug + sex+ I(age^2), data = pbc)
summary(model_cox_1)Call:
coxph(formula = Surv(years, status2) ~ drug + sex + age, data = pbc)
n= 1945, number of events= 725
coef exp(coef) se(coef) z Pr(>|z|)
drugD-penicil -0.130201 0.877919 0.075958 -1.714 0.0865 .
sexfemale -0.547998 0.578106 0.093038 -5.890 3.86e-09 ***
age 0.042168 1.043070 0.003966 10.634 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
drugD-penicil 0.8779 1.1391 0.7565 1.0188
sexfemale 0.5781 1.7298 0.4817 0.6937
age 1.0431 0.9587 1.0350 1.0512
Concordance= 0.602 (se = 0.012 )
Likelihood ratio test= 151 on 3 df, p=<2e-16
Wald test = 157.8 on 3 df, p=<2e-16
Score (logrank) test = 157.8 on 3 df, p=<2e-16Otra cosa importante que hay que considerar es que pueden haber interacciones entre los datos, es decir que dos variables o más se multipliquen. Ejemplos de estos se pueden encontrar, nuevamente, en este post.
Todo esto esta muy interesante, pero no hay que olvidar que hay supuestos que tenemos que verificar, y en este caso tenemos que ver la proporcionalidad entre los riesgos de los individuos. Esto es sencillo verlo con las gráficas que podemos obtener de las funciones de supervivencia de Kaplan-Meir realizando una estratificación de acuerdo a la variable que estemos analizando. Además de esto tenemos otras opciones como las pruebas y gráficos basados en los residuos de Schoenfeld. Para ver más acerca de esto se deja el siguiente enlace y una mención al paquete survminer, ya que en este se encuentran muchas funciones especializadas en todo esto. Aquí se deja el enlace a la Sheet Cheat para tenerla a la mano.
Se deja al lector comprobar si los supuestos son validados en los anteriores modelos o con otra combinación de variables. Para determinar si un modelo se ajusta mejor a los datos que otro podemos hacer la prueba del cociente de verosimiltudes para ver si entre los dos modelos existe una diferencia significante. Hay que tomar en cuenta que estos modelos deben diferir muy poco, como un interacción adicional o una variable extra considerada. Ya con dichos modelos, se puede utilizar la función anova(fit_null, fit_alt) donde fit_null será el modelo que se esta considerando para la hipótesis nula y se desea comprobar que este es distinto al modelo fit_alt.
Además de todo lo anterior habría que considerar otras cosas que podrían afectar nuestros resultados para bien o para mal
- Limpieza de datos en general
- Considerar más covariables o transformaciones no lineales entre ellas.
- El logaritmo de una variable reduce su varianza.
- Puede existir colinealidad. Para este caso véase más acerca del argumento
singular.oken la funcióncoxph()
Aquí se dejan algunos enlaces interesantes como un resumen a todo lo que hemos visto o sobre este último tema
- Diagnostics for the Cox model.
- Survival Analysis in R Companion.
- Análisis de Supervivencia en R..
- Survival Analysis in R For Beginners
Y para concluir, aquí se dejen algunas cosas que no se vieron en este trabajo.
- La prueba Peto-Peto logrank test es un test derivado de la prueba Logrank. Aquí se deja un poco de información al respecto.
- En el siguiente artículo se comparan los diferentes métodos de estimación sobre los parámetros de Cox.
- ATF: Accelerated failure time model. Otro modelo propuesto que considera covariables al modelo como una alternativa al modelo de Cox. Aquí otro archivo.
- El modelo de Cox tiene una relación con la regresión Poisson. En el siguiente artículo se deja una aplicación de esta interesante relación: Approximating the Cox Model.
Ejercicios
Escribe el código para obtener cada uno de los siguientes puntos
Obtener las funciones de densidad mediante la función de supervivencia de los ejemplos vistos en la sección Análisis de Supervivencia.
Obtener las funciones de densidad y supervivencia mediante la función \(h(t)\) de los ejemplos vistos en la sección Análisis de Supervivencia.
Obtén la función de riesgo acumulado para ambas funciones de supervivencia de los ejemplos vistos en la sección Análisis de Supervivencia pero con la forma alternativa de calcular \(H(t)\).
Crea una función que, dada la información correspondiente a la función de supervivencia en el caso discreto y el tiempo correspondiente, se obtenga la esperanza, la varianza y un marco de datos con la información correspondiente a la media residual.
Crea una función que, dada la información de una función de supervivencia discreta y su tiempo, devuelve algún o algunos cuantiles que determine el usuario.
Toma alguno de los ejemplos vistos en la sección Modelos paramétricos y obtén la función de supervivencia de manera discreta para ver como se aproximan los parámetros poblaciones con tus datos a los teóricos.
Con los ejemplos vistos en la sección Modelos paramétricos, comprueba las distintas igualdades que se tienen entre las funciones \(f(t)\), \(S(t)\), \(h(t)\) y \(H(t)\).
Crea un simple programa que, al menos con las distribuciones vistas en las sección Modelos paramétricos, se obtengan los resultados correspondientes a \(S(t)\), \(h(t)\) y \(H(t)\), además de sus parámetros poblacionales.
Utilizando los resultados teóricos, comprobar los resultados referentes a los parámetros poblacionales de la sección Modelos paramétricos obtenidos con aproximaciones numéricas.
Calcula la asimetría (Skewness) con un conjunto de datos simulados de la distribución gamma vista en la sección Modelos paramétricos.
Con el ejemplo de la distribución exponencial de la sección Verosimilitud C&T, crea una función donde puedas obtener los resultados del estimador máximo verosimil sin utilizar métodos numéricos, con dos parámetros: uno que indique la probabilidad de obtener datos censurados en un vector como
right_censured_Tby otro para determinar un tamaño de muestra \(n\).Con base al ejercicio anterior, crea una gráfica de lo que sucede al variar el tamaño de muestra con el estimador máximo verosimil.
Utilizar las ecuaciones vistas en el capítulo 6 del Bookdown referentes al caso de una distribución \(Weibull\) para datos censurados por la derecha tipo 1 para calcular los estimadores máximo verosimil de esta distribución con valores simulados. Calcule estos mediante métodos numéricos también.
Calcule la distancia Kullback Leibler para las distintas funciones de distribución que se vieron en la sección Verosimilitud C&T.
Obtén mediante el método de \(K-M\) la tabla de ejemplo referente a este método del Bookdown.
Crea una animación donde, con el ejemplo de los datos de myeloma vistos en la sección Modelos no paramétricos, se vea el comportamiento de la estimación de \(S(t)\) por medio del método actuarial a medida que se va disminuyendo el tamaño de los intervalos.
Realiza más segregaciones con los datos del myeloma y determina si existe algún factor relevante en la base de datos que modifique drásticamente la función de supervivencia.
Investiga la paquetería
survminery utiliza alguna de las funciones para obtener una gráfica interactiva de los datos de myeloma por sexo.Utiliza alguno de los datos contenidos en el paquete
KMsurvy realiza un análisis de supervivencia mediante el método actuarial y \(K-M\).Crea una función que, dada la función de supervivencia estimada por el método \(K-M\) y el método actuarial, se calculen todos los resultados que se vieron en la sección Análisis de supervivencia y Parámetros poblaciones.
Simula algunos datos de la sección Modelos Paramétricos de acuerdo a cada problema y utiliza la función
flexsurv::flexsurvregpara crear un ajuste a cada uno de los modelos correspondientes.Crea una función que determine la curva de supervivencia más parecida para cada uno de los estratos de acuerdo a la prueba \(Log-Rank\) con los datos de Billboard. ***
Klein, John P, and Melvin L Moeschberger. 2006. Survival Analysis: Techniques for Censored and Truncated Data. Second. Springer Science & Business Media.
Wisniewski, Piotr Marian, and Gabriel Velasco Sotomayor. 2001. Problemario de Probabilidad. 519.2076 W5P7.
A work by Carlos Vásquez
carlosfvasquez@ciencias.unam.mx