5 Estimación con datos completos y censurados
La teoría que veremos en esta sección, usualmente se denomina “Análisis de sobrevivencia”. La aplicación típica es conocer el tiempo de vida de un grupo de personas, bajo ciertas condiciones de riesgo. Por ejemplo, la esperanza de vida de un grupo de fumadores vs no fumadores. Sin embargo, los mismo modelos se pueden usar cuando se tienen montos de pagos en un seguro.
Usaremos ambas aplicaciones en esta sección.
5.1 Datos completos
Primero recuerde que la función de distribución se define como,
\[\begin{equation*} F_n(x) = \frac{1}{n}\sum_{i=1}^{n} 1_{X_i\leq x} \end{equation*}\]
y la función de sobrevivencia es \(S_n(x) = 1 - F_n(x)\). La tasa acumulativa de riesgo es
\[\begin{equation*} H(x)=-\ln S(x) \end{equation*}\]
Si \(S(x)\) es diferenciable, entonces \[\begin{equation*} H^{\prime}(x)=-\frac{S^{\prime}(x)}{S(x)}=\frac{f(x)}{S(x)}=h(x), \end{equation*}\] \[\begin{equation*} H(x)=\int_{-\infty}^x h(y) d y . \end{equation*}\]
Por lo tanto obtenemos la relacione \(F(x)=1-S(x)=1-e^{-H(x)}\)
Empíricamente para estimar \(F_n(x)\) de \(X_1, \dots, X_n\) tome \(y_1<y_2<\cdots<y_k\) los cuales son \(k\leq n\) valores únicos en la muestra. Sea \(s_j\) el número de veces que la observación \(y_j\) de modo que \(\sum_{j=1}^k s_j=n\). Adicional, defina el conjunto de riesgo \(r_j=\sum_{i=j}^{k}s_j\) como el número de observaciones igual o más grandes que \(y_j\). Entonces,
\[\begin{equation*} F_n(x)=\left\{\begin{array}{ll} 0, & x<y_1 \\ 1-\frac{r_j}{n}, & y_{j-1} \leq x<y_j, j=2, \ldots, k, \\ 1, & x \geq y_k . \end{array}\right. \end{equation*}\]
Podemos además definir empíricamente el estimador de Nelson-Aalen (NelsonTheory1972?; AalenNonparametric1978?) como
\[\begin{equation*} \hat{H}(x)= \begin{cases} 0, & x<y_1 \\ \displaystyle \sum_{i=1}^{j-1} \frac{s_i}{r_i}, & y_{j-1} \leq x<y_j, j=2, \ldots, k \\ \displaystyle \sum_{i=1}^k \frac{s_i}{r_i}, & x \geq y_k \end{cases} \end{equation*}\]
Para el caso de datos agrupados, se tiene la fórmula \[\begin{equation*} F_n(x)=\frac{c_j-x}{c_j-c_{j-1}} F_n\left(c_{j-1}\right)+\frac{x-c_{j-1}}{c_j-c_{j-1}} F_n\left(c_j\right), c_{j-1} \leq x \leq c_j . \end{equation*}\] donde los \(c_j\)’s son las fronteras de cada grupo. Esta fórmula conecta cada frontera con una línea recta.
Para el caso de la densidad se tiene que
\[\begin{equation*} f_n(x)=\frac{F_n\left(c_j\right)-F_n\left(c_{j-1}\right)}{c_j-c_{j-1}}=\frac{n_j}{n\left(c_j-c_{j-1}\right)}, c_{j-1} \leq x<c_j . \end{equation*}\]
Considere la siguiente tabla
tabla_c <- read.csv(text = "a,b,payments
0,7500,99
7500,17500,42
17500,32500,29
32500,67500,28
67500,125000,17
125000,300000,9
300000,Inf,3")
knitr::kable(tabla_c)| a | b | payments |
|---|---|---|
| 0 | 7500 | 99 |
| 7500 | 17500 | 42 |
| 17500 | 32500 | 29 |
| 32500 | 67500 | 28 |
| 67500 | 125000 | 17 |
| 125000 | 300000 | 9 |
| 300000 | Inf | 3 |
(x <- grouped.data(Group = c(0, tabla_c$b), Frequency = tabla_c$payments)) Group Frequency
1 (0, 7500] 99
2 (7500, 17500] 42
3 (17500, 32500] 29
4 (32500, 67500] 28
5 (67500, 125000] 17
6 (125000, 300000] 9
7 (300000, Inf] 3(Fo <- ogive(x))Ogive for grouped data
Call: ogive(x = x)
x = 0, 7500, 17500, ..., 3e+05, Inf
F(x) = 0, 0.43612, 0.62115, ..., 0.98678, 1plot(Fo)
5.2 Datos censurados
Usaremos las siguientes definiciones de ahora para las observaciones:
- Truncada desde abajo o a la izquierda en \(d\): Si cuando está en o por debajo de \(d\) no se registra, pero cuando está por encima de \(d\) se registra en su valor observado.
- Truncada desde arriba o a la derecha en \(u\): Si cuando está en o por encima de \(u\) no se registra, pero cuando está debajo de \(u\) se registra en su valor observado.
- Censurada desde abajo o por la izquierda en \(d\): Si cuando está en o por debajo de \(d\) se registra como igual a \(d\), pero cuando está por encima de \(d\) se registra en su valor observado.
- Censura desde arriba o por la derecha) en \(u\): Si cuando está en o por encima de \(u\) se registra como igual a \(u\), pero cuando está por debajo de \(u\) se registra en su valor observado.
Estos son ejemplos aplicados al curso:
- Un deducible ordinario de \(d\) es un ejemplo de dato truncado por la izquierda. Cuando el asegurado tiene una pérdida menor que \(d\), está no se reporta a la compañía.
- Un límite de póliza de \(u\), es un ejemplo de dato censurado por la derecha. Cuando el monto de la pérdida es igual o superior a \(u\), se registra solo el valor de \(u\).
- En estudios clínicos la personas entran a cierta edad, por lo que los eventos antes de la edad \(d\) son truncados por la izquierda. Cuando el estudio termina, todos los eventos son censurados a la derecha a la edad \(u\).
Para construir el estimador empírico, necesitamos calcular el grupo de riesgo en cada momento \(j\). Para esto defina
- \(d_j\) como el valor de truncamiento para la observación \(j\). Si no hay truncamiento entonces \(d_j=0\).
- \(y_1<y_2<\cdots<y_k\) valores únicos de los \(x_j\)’s que aparece en la muestra. Aquí \(k\) debe ser menor o igual al número de datos no censurados.
- \(s_j\) es el número de veces que la observación \(y_j\) aparece en la muestra.
- \(r_j\) es el número de datos “vivos” en la muestra incluyendo los censurados. En otras palabras: \[\begin{equation*} r_j= (\text {número de } d_i <y_j)- (\text {número de } x_i <y_j)-(\text{ número de } u_i <y_j) \end{equation*}\]
En otras palabras, son todos los individuos truncados antes de \(y_j\), menos todos los individuos que tuvieron un evento observado antes de \(y_j\), menos las censuras hasta \(y_j\).
\[\begin{equation*} \begin{aligned} r_j=& r_{j-1}+\left(\text { number of } d_i \text { s between } y_{j-1} \text { and } y_j\right) \\ &-\left(\text { number of } x_i \text { s equal to } y_{j-1}\right) \\ &-\left(\text { number of } u_i \text { s between } y_{j-1} \text { and } y_j\right) \end{aligned} \end{equation*}\]
Suponga que se tiene una póliza de vida a 5 años plazo. Se tienen las condiciones de muerte \(d\), cancelación \(s\) o expiración \(e\).
Usemos el gráfico anterior para calcular los dos primeros conjuntos de riesgo:
\(r_1\): Tenemos \(y_1=0.8\). Se observa que el número de \(d_i\)’s menores a 0.8 son 32. Luego, el número de eventos \(x_i\)’s menor estrictos a 0.8 es 0, y finalmente el número de censuras \(u_i\)’s menores estrictas a 0.8 es 2. Por lo tanto las personas en riesgo son \(r_1=32-0-2 =30\). En el momento \(y_1\) ocurrió exactamente 1 evento \(s_1\).
\(r_2\) Para \(y_2=2.9\), 3 nuevas pólizas ingresaron por los que el número de \(d_i\) es 35. Ahora, los eventos \(x_i\) pasados son 1, y las censuras pasadas son 8. Entonces \(r_2=35-1-8= 26\). A ese momento ocurrieron exactamente 2 eventos \(s_2\).
5.3 Estimador de Kaplan-Meier y Nelson Allen
El estimador de Kaplan-Meier se define como, \(S(0)=1\). En el punto \(y_1\), habían \(r_1\) personas con posibilidades de morir, de las cuales \(s_1\) lo hicieron. Así que tenemos una tasa de \(\frac{r_1-s_1}{r_1}\). Continuado así con cada punto obtenemos lo siguiente
\[\begin{equation*} S_n(t)= \begin{cases}1, & 0 \leq t<y_1, \\ \prod_{i=1}^{j-1}\left(\frac{r_i-s_i}{r_i}\right), & y_{j-1} \leq t<y_j, j=2, \ldots, k, \\ \prod_{i=1}^k\left(\frac{r_i-s_i}{r_i}\right) \text { or } 0, & t \geq y_k .\end{cases} \end{equation*}\]
Una versión alternativa al estimador de Kaplan-Meier es el de Nelson-Allen. En este caso estima directamente \(H(t)\) en lugar de \(S(t)\).
\[\begin{equation*} s(t)=\int_0^t r(u) h(u) d u . \end{equation*}\] \[\begin{equation*} d s(t)=r(t) h(t) d t . \end{equation*}\] \[\begin{equation*} \frac{d s(t)}{r(t)}=h(t) d t . \end{equation*}\] \[\begin{equation*} \int_0^t \frac{d s(u)}{r(u)}=\int_0^t h(u) d u=H(t) . \end{equation*}\]
Ahora reemplace \(s(t)\) por su versión empírica \(\hat{s}(t)\).
\[\begin{equation*} \sum_{t_i \leq t} \frac{s_i}{r_i} \end{equation*}\]
\[\begin{equation*} \hat{H}(t)= \begin{cases} 0, & 0 \leq t<y_1 \\ \displaystyle \sum_{i=1}^{j-1} \frac{s_i}{r_i}, & y_{j-1} \leq t<y_j, j=2, \ldots, k \\ \displaystyle \sum_{i=1}^k \frac{s_i}{r_i}, & t \geq y_k\end{cases} \end{equation*}\]
library(tidyverse)
library(actuar)
library(survival)
library(insuranceData)
library(ggfortify)
library(survminer)
library(goftest)Usaremos la tabla en el seiguiente linkg
tabla_d <- read.csv(file = "./data/tabla_d.csv")knitr::kable(tabla_d)| policy | d | obs | event |
|---|---|---|---|
| 1 | 0.0 | 0.1 | s |
| 2 | 0.0 | 0.5 | s |
| 3 | 0.0 | 0.8 | s |
| 4 | 0.0 | 0.8 | d |
| 5 | 0.0 | 1.8 | s |
| 6 | 0.0 | 1.8 | s |
| 7 | 0.0 | 2.1 | s |
| 8 | 0.0 | 2.5 | s |
| 9 | 0.0 | 2.8 | s |
| 10 | 0.0 | 2.9 | d |
| 11 | 0.0 | 2.9 | d |
| 12 | 0.0 | 3.9 | s |
| 13 | 0.0 | 4.0 | d |
| 14 | 0.0 | 4.0 | s |
| 15 | 0.0 | 4.1 | s |
| 16 | 0.0 | 4.8 | d |
| 17 | 0.0 | 4.8 | s |
| 18 | 0.0 | 4.8 | s |
| 19 | 0.0 | 5.0 | e |
| 20 | 0.0 | 5.0 | e |
| 21 | 0.0 | 5.0 | e |
| 22 | 0.0 | 5.0 | e |
| 23 | 0.0 | 5.0 | e |
| 24 | 0.0 | 5.0 | e |
| 25 | 0.0 | 5.0 | e |
| 26 | 0.0 | 5.0 | e |
| 27 | 0.0 | 5.0 | e |
| 28 | 0.0 | 5.0 | e |
| 29 | 0.0 | 5.0 | e |
| 30 | 0.0 | 5.0 | e |
| 31 | 0.3 | 5.0 | e |
| 32 | 0.7 | 5.0 | e |
| 33 | 1.0 | 4.1 | d |
| 34 | 1.8 | 3.1 | d |
| 35 | 2.1 | 3.9 | s |
| 36 | 2.9 | 5.0 | e |
| 37 | 2.9 | 4.8 | s |
| 38 | 3.2 | 4.0 | d |
| 39 | 3.4 | 5.0 | e |
| 40 | 3.9 | 5.0 | e |
tabla_d <- tabla_d %>%
mutate(delta = case_when(
event == "d" ~ 1,
event %in% c("s", "e") ~ 0
))El estimador de Kaplan Meier se construye de la siguiente forma
surv1 <-
Surv(
time = tabla_d$d,
time2 = tabla_d$obs,
event = tabla_d$delta
)
surv1 [1] (0.0,0.1+] (0.0,0.5+] (0.0,0.8+] (0.0,0.8] (0.0,1.8+] (0.0,1.8+]
[7] (0.0,2.1+] (0.0,2.5+] (0.0,2.8+] (0.0,2.9] (0.0,2.9] (0.0,3.9+]
[13] (0.0,4.0] (0.0,4.0+] (0.0,4.1+] (0.0,4.8] (0.0,4.8+] (0.0,4.8+]
[19] (0.0,5.0+] (0.0,5.0+] (0.0,5.0+] (0.0,5.0+] (0.0,5.0+] (0.0,5.0+]
[25] (0.0,5.0+] (0.0,5.0+] (0.0,5.0+] (0.0,5.0+] (0.0,5.0+] (0.0,5.0+]
[31] (0.3,5.0+] (0.7,5.0+] (1.0,4.1] (1.8,3.1] (2.1,3.9+] (2.9,5.0+]
[37] (2.9,4.8+] (3.2,4.0] (3.4,5.0+] (3.9,5.0+]kaplan_meier <- survfit(surv1 ~ 1)
summary(kaplan_meier)Call: survfit(formula = surv1 ~ 1)
time n.risk n.event censored survival std.err lower 95% CI upper 95% CI
0.8 30 1 3 0.967 0.0328 0.905 1.000
2.9 26 2 5 0.892 0.0589 0.784 1.000
3.1 26 1 0 0.858 0.0659 0.738 0.997
4.0 26 2 3 0.792 0.0755 0.657 0.955
4.1 23 1 1 0.758 0.0797 0.616 0.931
4.8 21 1 3 0.721 0.0837 0.575 0.906autoplot(kaplan_meier)
ggsurvplot(kaplan_meier, data = tabla_d, conf.int = TRUE)
ggsurvplot(kaplan_meier, data = tabla_d, conf.int = TRUE, risk.table = TRUE)
El de Nelson-Allen
nelson_allen <- survfit(surv1 ~ 1, type = "fleming-harrington")
summary(nelson_allen)Call: survfit(formula = surv1 ~ 1, type = "fleming-harrington")
time n.risk n.event censored survival std.err lower 95% CI upper 95% CI
0.8 30 1 3 0.967 0.0322 0.906 1.000
2.9 26 2 5 0.896 0.0571 0.790 1.000
3.1 26 1 0 0.862 0.0642 0.745 0.997
4.0 26 2 3 0.798 0.0736 0.666 0.956
4.1 23 1 1 0.764 0.0779 0.626 0.933
4.8 21 1 3 0.729 0.0820 0.584 0.908autoplot(nelson_allen)
ggsurvplot(nelson_allen, data = tabla_d, conf.int = TRUE)
ggsurvplot(nelson_allen, data = tabla_d, conf.int = TRUE, risk.table = TRUE)
Combinemos los dos resultados para comparar
ggsurvplot_combine(
list(kaplan_meier, nelson_allen),
data = tabla_d,
legend.labs = c("Kaplan Meier", "Nelson Allen")
)Warning: `select_()` was deprecated in dplyr 0.7.0.
Please use `select()` instead.
ggsurvplot_combine(
list(kaplan_meier, nelson_allen),
data = tabla_d,
legend.labs = c("Kaplan Meier", "Nelson Allen"),
fun = "cumhaz"
)
5.4 Media, varianza e intervalos de confianza
De los cursos pasados de estadística, no es difícil probar que
\[\begin{align*} \mathrm{E}\left[S_n(x)\right]&=\mathrm{E}\left(\frac{Y}{n}\right)=\frac{n S(x)}{n}=S(x) \\ \operatorname{Var}\left[S_n(x)\right]&=\operatorname{Var}\left(\frac{Y}{n}\right)=\frac{S(x)[1-S(x)]}{n} \\ \widehat{\operatorname{Var}}\left[S_n(x)\right]&=\frac{S_n(x)\left[1-S_n(x)\right]}{n} \end{align*}\]
Para el caso de Kaplan-Meier, se puede probar (no lo haré por tema de tiempo) que
\[\begin{equation} \mathrm{E}\left[\hat{S}\left(y_j\right)\right] =\frac{S\left(y_j\right)}{S\left(y_0\right)} \end{equation}\]
Usando aproximaciones de la varianza, se puede encontrar la aproximación de Greenwood.
\[\begin{equation} \widehat{\operatorname{Var}}\left[S_n\left(y_j\right)\right] \doteq S_n\left(y_j\right)^2 \sum_{i=1}^j \frac{s_i}{r_i\left(r_i-s_i\right)} \end{equation}\]
Otra opción es usar el método delta con la función \(g(x)=\ln(-\ln(x))\). La varianza de \(S_n(x)\) se puede aproximar con
\[\begin{equation} \left\{g^{\prime}\left[\mathbf{E}\left(S_n(t)\right)\right]\right\}^2 \operatorname{Var}\left[S_n(t)\right]=\frac{\operatorname{Var}\left[S_n(t)\right]}{\left[S_n(t) \ln S_n(t)\right]^2} \end{equation}\]
y el intervalo de confianza al 95% para \(\theta=\ln(-\ln(S(t)))\)
\[\begin{equation} \ln \left[-\ln S_n(t)\right] \pm 1.96 \frac{\sqrt{\widehat{\operatorname{Var}}\left[S_n(t)\right]}}{S_n(t) \ln S_n(t)} \end{equation}\]
Sin embargo este intervalo debe ser vuelto a transformar con \(S(t)=\exp(-e^{\theta})\).
La aproximación de Greenwood y la basada en el método delta también funciona para el estimador de Nelson-Allen,
Para el caso de la aproximación de Greenwood, se tiene que \[\begin{equation} \widehat{\operatorname{Var}}\left[\hat{H}\left(y_j\right)\right]\doteq \sum_{i=1}^j \frac{s_i}{r_i^2} \end{equation}\]
\[\begin{equation} \hat{H}(t) \pm z_{\alpha / 2} \sqrt{\widehat{\operatorname{Var}}\left[\hat{H}\left(y_j\right)\right]} \end{equation}\]
Para el método delta, el intervalo de confiaza quedaría
\[\begin{equation} \hat{H}(t) U \text {, where } U=\exp \left[\pm \frac{z_{\alpha / 2} \sqrt{\widehat{\operatorname{Var}}\left[\hat{H}\left(y_j\right)\right]}}{\hat{H}(t)}\right] \text {. } \end{equation}\]
```{r}kaplan_meier_plain <- survfit(surv1 ~ 1, conf.type = “plain”) kaplan_meier_loglog <- survfit(surv1 ~ 1, conf.type = “log-log”)
ggsurvplot_combine( list(kaplan_meier_plain, kaplan_meier_loglog), data = tabla_d, legend.labs = c(“IC Normal”, “IC Log-Log”), conf.int = TRUE ) ```