Introducción

El análisis estadístico más que tratarse de variables aleatorias y predicciones se trata de comparar, en este caso se ha generado una función que de manera gráfica ilustra la distancia máxima en el tiempo entre dos curvas de supervivencia y con esto darse una idea de a lo más que tanto se mejora (empeora) la situación entre ambos grupos.

dcurvG

Primero que nada esta función necesita de que la base de datos con la que se vaya a trabajar haya sido previamente tratada y ajustada a un modelo de supervivencia de su agrado. Después, se van a utilizar los vectores del modelo de supervivencia de los datos. Es altamente recomendado el no trabajar con la base de datos de manera directa.

Forma de uso

dcurvG()

dcurvG(t1,t2,s1,s2,...)

#Uso por default
dcurvG(t1,t2,s1,s2,et1="Grupo 1", et2= "Grupo 2", ns=0)

Argumentos

Los argumentos son los siguientes:

  • t1, t2 son los vectores de tiempo del estrato 1 y del estrato 2 respectivamente
  • s1, s2 son los vectores de supervivencia del estrato 1 y del estrato 2 respectivamente
  • et1 y et2 son las etiquetas que se le pondrán en la gráfica al primer y segundo estrato respectivamente.
  • ns se trata del nivel de significancia de las distancias, es decir, a partir de qué probabilidad de supervivencia (es decir \(0\leq{ns}<1\) ) nos interesa comparar las curvas.

Código

dcurvG <- function(t1,t2,s1,s2,et1="grupo 1",et2= "grupo 2",ns=0){
  
  #Primero... si diferentes tiempos tienen una misma probabilidad en la cola de la función de supervivencia 
  ultimos <- c(s1[length(s1)], s2[length(s2)])
  if(1 < sum(s1==ultimos[1])){
    suma <- sum(s1==ultimos[1])
    s1 <- s1[-((length(s1)-suma+2):length(s1))]
    t1 <- t1[-((length(t1)-suma+2):length(t1))]
  }
  if(1 < sum(s2==ultimos[2])){
    suma <- sum(s2==ultimos[2])
    s2 <- s2[-((length(s2)-suma+2):length(s2))]
    t2 <- t2[-((length(t2)-suma+2):length(t2))]
  }
  
  
  # funcion que aproxima el tiempo de supervivencia del estrato 1
  at1 <- approxfun(s1,t1)
  # funcion que aproxima el tiempo de supervivencia del estrato 2
  at2 <- approxfun(s2,t2)

  smerge <- unique(sort(c(s1,s2), decreasing = TRUE)) #survival vector merge
  t1m <- at1(smerge[smerge>=ns])
  t2m <- at2(smerge[smerge>=ns])
  d <- abs(t1m-t2m) #d = distancia

  dmax <- max(d, na.rm = TRUE)
  indice <- 0

  for(i in 1:length(smerge)){
    if(is.na(d[i])!=T){
      if(d[i] == dmax){
        indice <- i
        break
      }}
  }

plot(at1(smerge),smerge, main = paste("La maxima distancia entre curvas es ",format(round(dmax, 2), nsmall = 2)), 
     type="l", xlim = c(0,max(at1(smerge[length(smerge)]),at2(smerge[length(smerge)]),na.rm = TRUE)), 
     xlab="Tiempo", ylab="Probabilidad de Supervivencia")
lines(at2(smerge),smerge, col="red")
legend(max(at1(smerge[length(smerge)]),at2(smerge[length(smerge)]),na.rm = TRUE),1,c(et1, et2), 
       col = c("black", "red"), seg.len = 1, bty = "n",lty = 1,xjust = 1, yjust = 1)
lines(c(at1(smerge[indice]),at1(smerge[indice])),c(-2,1), lty="dotted", col="bisque4")
lines(c(at2(smerge[indice]),at2(smerge[indice])),c(-2,1), lty="dotted", col="bisque4")
lines(c(-100,max(at1(smerge[length(smerge)]),at2(smerge[length(smerge)]),na.rm = TRUE)),
      c(smerge[indice],smerge[indice]), lty="dotted", col="bisque4")

}

Ejemplos

Para los ejemplos se requiere usar las siguientes paqueterías

require(survival)
require(survminer)

Para construir el modelo de supervivencia se usará la función Surv de la paquetería survival y para hacer la función de supervivenvcia a partir del modelo se usará la función surv_fit de la paquetería survminer.

1. Estudio de cáncer de pulmón de veteranos

Es una muestra de pruebas de dos tipos de tratamientos de cáncer de pulmón incluída en la paquetería survival. Se desea comparar los resultados de los tratamientos realizados a los veteranos, para esto se ajusta el modelo de supervivencia con el tiempo (time) y el estatus de censura (status) agrupado de acuerdo a los tratamientos (trt).

df <- veteran
km_model_trt <- surv_fit(Surv(time, status)~trt, df)
ggsurvplot
km_model_trt$strata # nos dice cuantos datos corresponden a la curva 1 y a la curva 2
#en este caso en el primer estrato (trt=1) hay 61 obs
#y en el estrato 2 (trt=2) hay 53

#Entonces la curva de supervivencia del strato 1 está en
t1 <- km_model_trt$time[1:61]
s1 <- km_model_trt$surv[1:61]

#La segunda curva de supervivencia está en
t2 <- km_model_trt$time[(61+1):(61+53)] #53 son las observaciones del estrato 2 (trt=2)
s2 <- km_model_trt$surv[(61+1):(61+53)] 

#Al final se evalúa con la función
dcurvG(t1,t2,s1,s2)

En la gráfica se muestra que aquellos que optaron por el tratamiento 2 tienen una mayor supervivencia pero con una probabilidad mínima (menos del 5%) así que como la curva de supervivencia del grupo 1 es mejor en probabilidades arriba del 20% se usará ese valor(0.2) para el nivel de significancia (ns). Además se le agregarán etiquetas para denotar de qué tratamiento se trata.

dcurvG(t1,t2,s1,s2,ns=0.2,"trt 1","trt 2")
#nótese que da el mismo resultado si se escribe:
#dcurvG(t1,t2,s1,s2,"trt 1","trt 2",0.2)

La diferencia máxima en tiempos tignificativos es de 50.13 unidades de tiempo al casi 50%, de esta manera se puede notar que el tratamiento 1 da mejores resultados que el tratamiento 2.

2. Catéter en riñones

Se tiene la base de datos kidney incluída en la paquetería survival. Esta base de datos representa el tiempo de reaparición de una infección al ponerle un catéter a los pacientes con problemas de riñones que usan equipo de diálisis. El catéter puede ser removido por causas ajenas a la infección, en tal caso se censuran los datos. Nos fijamos en las variables: tiempo antes de remover el cateter (time) y si fue removido por la infección o por otras causas (status) para hacer el modelo de supervivencia al final se agrupa respecto a la variable sexo (sex). Se quiere ver que género tiene mejores probabilides de supervivencia después de remover el catéter.

km_model_sex <- surv_fit(Surv(time,status)~sex, kidney)
km_model_sex$strata
# En el grupo sex=1 (male) hay 18 observaciones y en el segundo grupo (female) hay 50

#Entonces la curva de supervivencia del strato 1 está en
t1 <- km_model_sex$time[1:18]
s1 <- km_model_sex$surv[1:18]

#La segunda curva de supervivencia está en
t2 <- km_model_sex$time[(18+1):(18+50)] #50 son las observaciones del estrato 2
s2 <- km_model_sex$surv[(18+1):(18+50)]

#Al final se evalúa en la curva
dcurvG(t1,t2,s1,s2,"Hombres","Mujeres")

Aquí se observa que la supervivencia de las mujeres siempre está por encima por la de los hombres por lo que no es necesario pedir un nivel de significancia más que para ver la distancia máxima arriba de esa medida.

3. Empleados

Se tiene una base de datos proveniente de recursos humanos donde se ve la información de los trabajadores como nivel de satisfacción, horas promedio mensuales, etcétera. Por el momento nos enfocaremos en tiempo trabajado en la compañía (time_spend_company) y si se fue (left) que es una variable de tiempo y la de estatus respectivamente para hacer el modelo de supervivencia y se ajustará de acuerdo a la variable salario (salary). Sea empleados el data frame.

empleados <- read.csv("HR_comma_sep.csv")
km_model_salary <- surv_fit(Surv(time_spend_company,left)~salary, empleados)
km_model_salary$strata
# Hay 8 observaciones en cada grupo, el orden es salario alto, bajo y medio 
# Nos fijaremos en los últimos dos estratos

#Entonces la curva de supervivencia de los trabajadores con salario medio son
t1 <- km_model_salary$time[(8+8+1):(8+8+8)] #Cada 8 es uno de los estratos basados en el salario
s1 <- km_model_salary$surv[(8+8+1):(8+8+8)]

#Entonces la curva de supervivencia de los trabajadores con salario bajo son
t2 <- km_model_salary$time[(8+1):(8+8)] 
s2 <- km_model_salary$surv[(8+1):(8+8)]

dcurvG(t1,t2,s1,s2,"salario medio","salario bajo") #Se evalúa la función

Se puede observar que la aproximación de la curva de supervivencia de salarios medios es mayor a la de bajos por un máximo de 1.35 unidades de tiempo.

Nota: la cola de los datos queda recta para diferentes tiempos cuando se hace el ajuste del modelo como se muestra en la siguiente gráfica.

La función se encarga de borrar esas observaciones y hacer la aproximación con probabilidades no repetidas.