Cargamos los datos de Myeloma y para visualizar los datos imprimimos las primeras entradas de estos mismos. Los paquetes que usaremos en esta practica son: stats, dplyr, openxlsx, ggplot2, scales, ggthemes.

data <- read.csv("Myeloma_data.csv", header = TRUE, sep = ",")
head(data)
##   Patient_number Survival_time Status Age Sex Bun   Ca   Hb Pcells Protein
## 1              1            12      1  66   1  25 14.6 14.6     18       1
## 2              2            52      0  66   1  13 12.0 12.0    100       0
## 3              3             6      1  53   2  15 11.4 11.4     33       1
## 4              4            40      1  69   1  10 10.2 10.2     30       1
## 5              5            10      1  65   1  20 13.2 13.2     66       0
## 6              6             7      0  57   2  12  9.9  9.9     45       0

Método Actuarial

Sabemos que el Método Actuarial divide el tiempo en intervalos y calcula la supervivencia en cada intervalo, la longitud del invertvalo depende de la frecuencia con que ocurre el suceso de interés.

Ahora bien, en este caso el intervalo que escogemos fue el de \(l_{j}\) que va de año en año (Valores en meses).

interv <- seq(1, 8, 1)
l_j <- seq(0, 100, 12)
data <- data%>%
  cbind(l_k = cut(data$Survival_time, breaks = l_j,labels = interv, right =FALSE))%>%
  mutate(l_k = as.numeric(l_k))

Ahora bien, necesitamos hacer desde cero las funciones para los calculos pedidos.

Necesitamos checar tanto el intervalo en el que estamos como el status del mismo (Ya que en este caso hay censura por la derecha y es importante hacer la distinción).

c_j <- data.frame(n = 0, l_k =0)

Esta función nos ayuda a poder filtrar los Status para ver los datos censurados

funcion_c_j <- function(x){
   u <- data%>%
     filter(l_k == x,
            Status == 0)%>%
     count()%>%
     mutate(l_k = x)
   
   c_j <- c_j %>%
     rbind(u)
     }
for (i in 1:8) {
  c_j <-funcion_c_j(i)
}

Para poder determinar el número de pacientes por en cada intervalo se hace la siguiente función

n_j <- data.frame(n = 0, l_k =0)
funcion_n_j <- function(x){
   u <- data%>%
     filter(l_k >= x)%>%
     count()%>%
     mutate(l_k = x)
   
   n_j <- n_j %>%
     rbind(u)
     }
for (i in 1:8) {
  n_j <- funcion_n_j(i)
}

Ya teniendo los pacientes de cada intervalo necesitamos ver las defunciones, hacemos la siguiente funcion que calcula las defunciones.

d_j <- data.frame(n = 0, l_k =0)
funcion_d_j <- function(x){
   u <- data%>%
     filter(l_k == x,
            Status == 1)%>%
     count()%>%
     mutate(l_k = x)
   
   d_j <- d_j %>%
     rbind(u)
     }
for (i in 1:8) {
  d_j <- funcion_d_j(i)
}

Ahora simplemenente almacenamos la funcion de supervivencia entera

actuarial <- data.frame(l_k = seq(1,8,1),
                    t = seq(0,84,12))%>%
  left_join(c_j%>%
              rename(c_j = n), by ="l_k")%>%
  left_join(d_j%>%
              rename(d_j = n), by ="l_k")%>%
  left_join(n_j%>%
              rename(n_j = n), by ="l_k")


actuarial <- actuarial%>%
  filter(l_k !=0)%>%
  rename(intervalo = l_k)%>%
  mutate(n_j_1 = round(n_j - (c_j/2), digits = 0),
         qj = d_j/n_j_1,
         pj = 1-qj)

Por ultimos solo falta obtener la funcion de supervivencia \(S(t)\), esto lo hacemos con los estimadores ya obtenidos

t <- seq(1,8,1)

m <- actuarial[1,8]
st <- data.frame(1, m)
for (i in 2:8) {
  m <- actuarial[i,8] * m
  
  j <-c(i, m)
  st <- st%>%
    rbind(j)
  }

actuarial <- actuarial%>%
  left_join(st%>%
              rename(intervalo = X1,
                     st = m), by = "intervalo")

Método Kaplan - Meier

El método de Kaplan-Meier calcula la supervivencia cada vez que un paciente muere. Da proporciones exactas de supervivencia debido a que utiliza tiempos de supervivencia precisos.

y <- data%>%
  arrange(Survival_time)

y <-unique(y$Survival_time)

Tenemos que hacer una función para obtener el numero de pacientes en riesgo dado un tiempo \(t\) para poder usarla para calcular la función final.

nk <- data.frame(n = 0, t =0)
nk_fun <- function(x){
   u <- data%>%
     filter(Survival_time >= y[x])%>%
     count()%>%
     mutate(t = y[x])
   
   nk <- nk %>%
     rbind(u)
     }
for (i in 1:28) {
  nk <- nk_fun(i)
}

Para calcular los fallos que han ocurrido en tiempo \(t\) tenemos la siguiente función

d_k <- data.frame(n = 0, t =0)

funcion_d_k <- function(x){
   u <- data%>%
     filter(Survival_time == y[x],
            Status == 1)%>%
     count()%>%
     mutate(t = y[x])
   
   d_k <- d_k %>%
     rbind(u)
     }
for (i in 1:28) {
  d_k <- funcion_d_k(i)
}

Juntamos en un dataframe los datos necesarios para poder realizar el calculo de la función final

kaplan_meier <- data.frame(y)%>%
  rename(t = y)%>%
  
  left_join(d_k%>%
              rename(d_k =n), by = "t")%>%
  left_join(nk%>%
              rename(nk =n), by = "t")%>%
  mutate(hk = 1-(d_k/nk))

Ahora bien, se hace una función que calcula la función de supervivencia segun Kaplan Meier y queda de la siguiente forma.

m <- kaplan_meier[1,4]
st_km <- data.frame(y[1], m)

for (i in 2:28) {
  m <- kaplan_meier[i,4] * m
  
  j <-c(y[i], m)
  st_km<- st_km%>%
    rbind(j)
  }

kaplan_meier <- kaplan_meier%>%
  left_join(st_km%>%
              rename(t = 1,
                     st_km = m), by = "t")

Vemos en la siguiente fráfica que se ve de una manera clara la diferencia entre ambos métodos, siendo estos

dat <- data%>%
  rename(intervalo = l_k)%>%
  left_join(actuarial, by = "intervalo")%>%
  left_join(kaplan_meier%>%
              rename(Survival_time= t), by = "Survival_time")%>%
  select(Survival_time,
         st,
         st_km)%>%
  arrange(Survival_time)

ggplot(dat) +
  geom_line(aes(Survival_time,st), color= "violet")+
  geom_point(aes(Survival_time,st),size =2,shape= 20, fill ="white")+
  geom_line(aes(Survival_time, st_km), color= "orange")+
  geom_point(aes(Survival_time, st_km),size =2,shape=20, fill ="white")+
  theme_economist()+
  scale_y_continuous(breaks = breaks_pretty(10))+
  ylab("S(t)")+
  xlab("Tiempo")+
  ggtitle("Tabla Actuariales")