library(dplyr)
library(ggplot2)
library(latex2exp)
library(knitr)
library(kableExtra)

El objetivo de esta práctica es generar desde cero los métodos no paramétricos para la función de Supervivencia, los métodos son:

  • Método Actuarial
  • Estimador Producto-Límite(Kaplan-Meier)

Método Actuarial:

Se divide el tiempo en intervalos y se calcula \(S(t)\) dentro de cada intervalo, el proceso es el siguiente:

Se divide el eje de tiempo en \(k+1\) intervalos, en este caso se toma como limite inferior el \(0\) y como límite superior \(\lceil \text{max}\{T_i\}\rceil\) donde \(T_i\) son los tiempos de fallo observados.

Cada intervalo es de la forma: \[l_j=(a_{j-1},a_j]\] Se calculan:

  • \(n_j\): número de individuos en riesgo(vivos o no censurados) al tiempo \(a_{j-1}\)
  • \(d_j\): número de fallas en \(l_j\)
  • \(d_j\): número de censuras en \(l_j\)
  • \(c_j^*\): número ajustado en riesgo \[n_j^*=n_j-\frac{c_j}{2}\]
  • \(q_j\): probabilidad de fallo dentro del intervalo \(l_j\), donde:

\[\hat q_j=1-\frac{d_j}{n_j-\frac{c_j}{2}}\]
* \(p_j=1-q_j\) es la probabilidad de sobrevivir hasta \(a_j\) despues de haver sobrevivido hasta \(a_{j-1}\).
Si no hay censuras en \(l_j\), entonces \(n_j^*=n_j\)

Finalmente se calcula \(\hat S(t)\) como:

\[\hat S_j=\prod_{i=1}^jp_j\] Se cargan los datos

myeloma<-read.csv("Myeloma_data.csv")
d<-data.frame(myeloma)

options(digits=4)

La siguiente función sigue el proceso descrito anteriormente y devuelve una tabla que contiene la informacion relevante sobre el método Actuarial.

#Funcion que calcula la supervivencia mediante el metodo actuarial
#d:conjunto de datos(dataframe)
#k:numero de intervalos
#La funcion devuelve una tabla($tabla) y un vector de tiempos($t)
metodo_actuarial<-function(d,k){
  #longitud de cada intervalo
  l=ceiling(max(d$Survival_time))/k
  
  #generamos los vectores que contienen los intervalos aj-1 y aj,
  #donde [-(k+1)] se usa para eliminar el ultimo elemento ya que no nos sirve
  aj_1=seq(0,l*k,l)[-(k+1)]
  aj=lead(seq(0,l*k,l))[-(k+1)]
  
  #obtenemos lo siguiente: 
  #nj: cantidad de supervivientes a tiempo aj-1
  #nj ajustado: ajuste sobre nj
  #dj: numero de fallas en el intervalo (aj-1,aj]
  #cj: numero de censuras en el intervalo (aj-1,aj]
  #pj: probabilidad de sobrevivir hasta aj, despues de haber sobrevivido hasta aj-1
  nj<-rep(0,k+1)
  nj_ajustado<-rep(0,k+1)
  dj<-rep(0,k+1)
  cj<-rep(0,k+1)
  pj<-rep(0,k+1)
  Sj<-rep(1,k+1)
  intervalo<-seq(1,k)
  for(i in 1:k){
    nj[i] = nrow(subset(d, Survival_time>aj_1[i]))
    dj[i] = nrow(subset(d, Survival_time>aj_1[i] & Survival_time<=aj[i] & Status==1))
    cj[i] = nrow(subset(d, Survival_time>aj_1[i] & Survival_time<=aj[i] & Status==0))
    nj_ajustado[i]=nj[i]-cj[i]/2
    pj[i] = 1-(dj[i]/nj_ajustado[i])
    Sj[i+1] = prod(pj[1:i])
    intervalo[i] = paste("(",round(aj_1[i],3),",",round(aj[i],3),"]",sep = "")
  }
  intervalo[k+1] = paste(">",aj[k],sep = "")

  tabla<-data.frame(intervalo,dj,cj,nj,nj_ajustado,Sj)
  #la funcion devuelve la tabla y un vector con los tiempos que sirve para graficar
  return(list(tabla=tabla,t=seq(0,l*k,l)[-(k)]))
}

Generemos la tabla con \(15\) intervalos.

k=15
metodo_act<-metodo_actuarial(d,k)

#le damos formato a la tabla
kable(metodo_act$tabla, caption = "Método Actuarial"
             , align = c('c', 'c', 'c', 'c', 'c', 'c'), 
      col.names=c("Intervalo","$d_j$", "$c_j$", "$n_j$", "$n^{*}_j$", "$S(t)$")) %>% 
  kable_styling(full_width = F, bootstrap_options = c("striped", "hover", "condensed"))
Método Actuarial
Intervalo \(d_j\) \(c_j\) \(n_j\) \(n^{*}_j\) \(S(t)\)
(0,6.067] 11 1 48 47.5 1.0000
(6.067,12.133] 7 4 36 34.0 0.7684
(12.133,18.2] 7 3 25 23.5 0.6102
(18.2,24.267] 2 0 15 15.0 0.4285
(24.267,30.333] 0 0 13 13.0 0.3713
(30.333,36.4] 1 0 13 13.0 0.3713
(36.4,42.467] 2 1 12 11.5 0.3428
(42.467,48.533] 0 0 9 9.0 0.2831
(48.533,54.6] 2 1 9 8.5 0.2831
(54.6,60.667] 0 1 6 5.5 0.2165
(60.667,66.733] 2 0 5 5.0 0.2165
(66.733,72.8] 0 0 3 3.0 0.1299
(72.8,78.867] 0 1 3 2.5 0.1299
(78.867,84.933] 0 0 2 2.0 0.1299
(84.933,91] 2 0 2 2.0 0.1299
>91 0 0 0 0.0 0.0000

La gráfica es la siguiente:

ggplot()+ 
  geom_step(aes(x=metodo_act$t,y=metodo_act$tabla$Sj[1:k], color = "Método-Actuarial"),size=1)+
  theme_bw()+
  labs(title="Función de Supervivencia", x = "Tiempo", y = TeX("$\\hat{S}(t)$")) +
  scale_color_manual(breaks=c("Método-Actuarial","K-M") ,values = c("Método-Actuarial"="steelblue","K-M"="#c45d16")) +
  theme(legend.position = "top", legend.direction = "horizontal", plot.title = element_text(hjust = 0.5))+
  labs(color = "Estimador")

Estimador Producto-Límite(Kaplan-Meier)

Se define el soporte \(\{u_1,u_2,...\}\) de los tiempos de fallo \(T_1,T_2,...,T_n\), dado que es posible obtener observacions repetidas.

Se calculan:

  • \(d_k\): número de tiempos de fallo ocurridos en \(T=u_k\)
  • \(c_k\): número de tiempos de fallo ocurridos en \(T=u_k\)
  • \(n_k\): número de individuos en riesgo a tiempo \(u_k\)

Se calcula \(\hat S(t)\) como:

\[\hat S(t)=\prod_{k:u_k\leq t}\left(1-\frac{d_k}{n_k}\right)\]

La siguiente función sigue el proceso descrito anteriormente y devuelve una tabla que contiene la informacion relevante sobre el método K-M.

#Funcion que calcula la supervivencia mediante el metodo Kaplan-Meier
#d:conjunto de datos(dataframe)
K_M<-function(d){
  #determinamos el soporte de los tiempos(uk)
  t<-sort(unique(d$Survival_time))
  l<-length(t)
  #se calculan:
  #dk:numero de tiempos de fallo iguales a uk
  #nk:numero de individuos en riesgo a tiempo uk
  #ck:numero de censuras a tiempo uk
  #Sk:supervivencia
  dk<-rep(0,l)
  nk<-rep(0,l)
  ck<-rep(0,l)
  Sk<-rep(1,l)
  pk<-rep(0,l)
  for(i in 1:(l)){
    dk[i] = nrow(subset(d, Survival_time==t[i] & Status==1))
    ck[i] = nrow(subset(d, Survival_time==t[i] & Status==0))
    nk[i] = nrow(subset(d, Survival_time>=t[i]))
    pk[i] = 1-dk[i]/nk[i]
    Sk[i] = prod(pk[1:i])
  }
  tabla<-data.frame(t,dk,nk,ck,Sk)  
  return(tabla)
}
metodo_KM<-K_M(d)

#le damos formato a la tabla resultante
kable(metodo_KM, caption = "Método K-M"
      , align = c('c', 'c', 'c', 'c', 'c'), 
      col.names=c("Tiempo","$d_k$", "$n_k$", "$c_k$", "$S(t)$")) %>% 
  kable_styling(full_width = F, bootstrap_options = c("striped", "hover", "condensed"))
Método K-M
Tiempo \(d_k\) \(n_k\) \(c_k\) \(S(t)\)
1 3 48 0 0.9375
3 0 45 1 0.9375
4 2 44 0 0.8949
5 4 42 0 0.8097
6 2 38 0 0.7670
7 0 36 1 0.7670
8 1 35 0 0.7451
10 4 34 1 0.6575
11 0 29 1 0.6575
12 2 28 1 0.6105
14 1 25 0 0.5861
15 1 24 1 0.5617
16 2 22 0 0.5106
17 1 20 0 0.4851
18 2 19 2 0.4340
23 1 15 0 0.4051
24 1 14 0 0.3761
36 1 13 0 0.3472
40 2 12 1 0.2893
50 1 9 0 0.2572
51 1 8 0 0.2250
52 0 7 1 0.2250
56 0 6 1 0.2250
65 1 5 0 0.1800
66 1 4 0 0.1350
76 0 3 1 0.1350
88 1 2 0 0.0675
91 1 1 0 0.0000

La gráfica es la siguiente:

ggplot()+ 
  geom_step(aes(x=metodo_KM$t,y=metodo_KM$Sk, color = "K-M"),size=1)+
  theme_bw()+
  labs(title="Función de Supervivencia", x = "Tiempo", y = TeX("$\\hat{S}(t)$")) +
  scale_color_manual(breaks=c("Método-Actuarial","K-M") ,values = c("Método-Actuarial"="steelblue","K-M"="#c45d16")) +
  theme(legend.position = "top", legend.direction = "horizontal", plot.title = element_text(hjust = 0.5))+
  labs(color = "Estimador")

Comparación

Realizamos una comparación gráfica de ambos modelos

ggplot()+ 
  geom_step(aes(x=metodo_act$t,y=metodo_act$tabla$Sj[1:k], color = "Método-Actuarial"),size=1)+
  geom_step(aes(x=metodo_KM$t,y=metodo_KM$Sk, color = "K-M"),size=1)+
  theme_bw()+
  labs(title="Función de Supervivencia", x = "Tiempo", y = TeX("$\\hat{S}(t)$")) +
  scale_color_manual(breaks=c("Método-Actuarial","K-M") ,values = c("Método-Actuarial"="steelblue","K-M"="#c45d16")) +
  theme(legend.position = "top", legend.direction = "horizontal", plot.title = element_text(hjust = 0.5))+
  labs(color = "Estimador")

Si usamos una mayor cantidad de intervalos para el método actuarial, la grafica de \(S(t)\) se “suavizará” más. Para \(k=30\) se obtiene la siguiente gráfica.

k=30
metodo_act<-metodo_actuarial(d,k)

ggplot()+ 
  geom_step(aes(x=metodo_act$t,y=metodo_act$tabla$Sj[1:k], color = "Método-Actuarial"),size=1)+
  geom_step(aes(x=metodo_KM$t,y=metodo_KM$Sk, color = "K-M"),size=1)+
  theme_bw()+
  labs(title="Función de Supervivencia", x = "Tiempo", y = TeX("$\\hat{S}(t)$")) +
  scale_color_manual(breaks=c("Método-Actuarial","K-M") ,values = c("Método-Actuarial"="steelblue","K-M"="#c45d16")) +
  theme(legend.position = "top", legend.direction = "horizontal", plot.title = element_text(hjust = 0.5))+
  labs(color = "Estimador")