#Descargar paquetes

library(openxlsx)
## Warning: package 'openxlsx' was built under R version 4.0.3
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(stats)

Descargamos los datos

datos <- read.csv("Myeloma_data.csv",
                  header = TRUE,
                  sep = ",")

#Metodo actuarial

El metodo actuarial consiste en dividir en intervalos de tiempo y modelar la supervivencia. El ejemplo que vamos a modelar cuenta con 48 participantes, el tiempo esta dado en meses.

El intervalo de tiempo escogido será de 12 meses en 12 meses. Ahora bien, vamos a colocar los intervalos.

lj <- seq(0,100, 12)
lab <- seq(1,8, 1
           )
'Nos apoyaremos de la funcion cut para poner los intervalos'
## [1] "Nos apoyaremos de la funcion cut para poner los intervalos"
datos <- datos%>%
  cbind(lk = cut(datos$Survival_time, breaks = lj, 
                 labels = lab, right = FALSE))%>%
  mutate(lk = as.numeric(lk))

De acuerdo con las indicaciones de la práctica, no podemos utilizar paquetes especializados, por lo que modelamos las funciones que nos ayuden al cálculo de los mismos.

Para el metodo actuarial es importante fijarse en los datos censurados, este experimento cuenta con una censura por la derecha debido a los pacientes que murieron despues de que concluyera el estudio. En el “status” los que se encuentran con 0 son los censurados.

cj <- data.frame(n = 0, lk =0)
'esta funcion nos ayuda a filtrar el status 0 y el intervalo en el que nos estamos colocando'
## [1] "esta funcion nos ayuda a filtrar el status 0 y el intervalo en el que nos estamos colocando"
cj_fun <- function(x){
   u <- datos%>%
     filter(lk == x,
            Status == 0)%>%
     count()%>%
     mutate(lk = x)
   
   cj <- cj %>%
     rbind(u)
     }
for (i in 1:8) {
  cj <-cj_fun(i)
}

Ahora bien vamos a determinar el número de pacientes en riesgo en el intervalo t.

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

Para este caso, vamos a determinar los pacientes que murieron en el intervalo t.

dj <- data.frame(n = 0, lk =0)
'la siguiente funcion ayuda a filtrar el tiempo t y el status 1'
## [1] "la siguiente funcion ayuda a filtrar el tiempo t y el status 1"
dj_fun <- function(x){
   u <- datos%>%
     filter(lk == x,
            Status == 1)%>%
     count()%>%
     mutate(lk = x)
   
   dj <- dj %>%
     rbind(u)
     }
for (i in 1:8) {
  dj <- dj_fun(i)
}

Ahora obtendremos el data frame que contiene todos los calculos para obtener la funcion de supervivencia.

m_act <- data.frame(lk = seq(1,8,1),
                    t = seq(0,84,12))%>%
  left_join(cj%>%
              rename(cj = n), by ="lk")%>%
  left_join(dj%>%
              rename(dj = n), by ="lk")%>%
  left_join(nj%>%
              rename(nj = n), by ="lk")


m_act <- m_act%>%
  filter(lk !=0)%>%
  rename(intervalo = lk)%>%
  mutate(nj_1 = round(nj - (cj/2), digits = 0),
         qj = dj/nj_1,
         pj = 1-qj)

Ahora obtendremos S(t) con los estimadores, la funcion modelada ayuda a realizar la multiplicacion iterada.

t <- seq(1,8,1)

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

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

#Método por estimador Kaplan - Meier

El metodo Kaplan-Meier calcula la supervivencia cada que un paciente muere. A diferencia del metodo actuarial, esta no utiliza intervalos, sino tiempos exactos de supervivencia.

Por lo cual ahora obtendremos todos los tiempos que aparecen en nuestros datos.

y <- datos%>%
  arrange(Survival_time)

y <-unique(y$Survival_time)

Obtendremos el numero de pacientes en riesgo al tiempo t.

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

Ahora bien, es importante calcular los fallos al tiempo t. Para este punto los que tienen Status 1 al tiempo t.

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

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

Obtenemos el data frame que contiene los tiempos y estimadores necesarios para el calculo de la funcion de supervivencia.

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

La funcion de supervivencia por le Metodo Kaplan Meier la obtenemos con una multiplicacion iterada de los hj, modelamos una funcion para esto filtrando por el tiempo y multiplicando con el valor anterior.

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

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

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

Obtenemos la grafica comparativa de ambos metodos.

Para realizarla vamos a juntar ambas funciones de supervivencia en el data frame original para obtener el tiempo t exacto por paciente.

En la grafica es clara la diferencia entre ambas, para el metodo actuarial observamos saltos importantes debido que utiliza intervalos, mientras que la Kaplan-Meier utiliza tiempos exactos de supervivencia.

En rojo observamos la S(t) por el Metodo Kaplan-Meier y en azul por el Metodo actuarial.

library(ggplot2)
library(scales)

dat <- datos%>%
  rename(intervalo = lk)%>%
  left_join(m_act, by = "intervalo")%>%
  left_join(m_km%>%
              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= "blue")+
  geom_point(aes(Survival_time,st),size =2,shape= 20, fill ="white")+
  geom_line(aes(Survival_time, st_km), color= "red")+
  geom_point(aes(Survival_time, st_km),size =2,shape=20, fill ="white")+
  theme_minimal()+
  scale_y_continuous(breaks = breaks_pretty(10))+
  ylab("S(t)")+
  xlab("tiempo")