data(veteran)

La finalidad de esta práctica es crear una función que nos de la distancia máxima entre dos curvas de supervivencia. Primero vamos a estimar los datos de las bases con la función survfit y los resultados obtenidos aquí los usaremos en la función.

km.veteran <- survfit(Surv(time,status)~trt,veteran)

Hay que tener algo en cuenta: como estamos estimando, no es una función de supervivencia “real”, por lo que puede contradecir algunas propiedades. Por esto, será conveniente definir un límite en el eje Y de nuestra gráfica (la probabilidad de supervivencia) hasta donde tomaremos en cuenta los valores. Este límite es arbitario y para darnos una idea de cual podría ser, primero usaremos la función ggsurvplot(), esto solo como ayuda visual

ggsurvplot(km.veteran)

Después de ver esta gráfica, parace que un límite superior adecuado podría ser 0.1. Acontinuación, mostramos la función, en esta metemos la información obtenida al estimar nuestra base, además del límite mencionado anteriormente. Usamos una función de R para estimar las funciones de supervivencia “inversas” (el eje X sería las probabilidades de supervivencia) y generamos un vector de probabilidades, que irá del 0 hasta el valor máximo de las probabilidades de supervivencia estimadas, pero usamos la función seq() para conseguir un vector de 1001 entradas, esto con la finalidad de hacer “un poco más continua nuestra función”. Como tenemos aproximaciones para los valores adentro del intervalo, no debe de causar mucho problema este enfoque.

Es importante mencionar que esta función está hecha para una estimación con dos estratos, por lo que hay que tener en cuenta esto a la hora de intentar utilizarla.

En general, lo demás no causa mucho problema. Ya hemos trabajado con ggplot() y es con esta función que hacemos la gráfica correspondiente.

d.max <- function(km.surv,km.time,km.strata,surv.min){
  s<-km.strata
  t.surv1 <- data.frame(tiempo=km.time[1:s[1]], surv=km.surv[1:s[1]])
  t.surv2 <- data.frame(tiempo=km.time[(s[1]+1):sum(s)], surv=km.surv[(s[1]+1):sum(s)])
  
  aprox_1 <- approxfun(t.surv1[,2],t.surv1[,1])
  aprox_2 <- approxfun(t.surv2[,2],t.surv2[,1])
  
  m <- max(t.surv1$surv[1],t.surv2$surv[1]) #el maximo entre los dos conjuntos de supervivencia
  
  p.eval <- seq(0,m,m/1000) #puntos que vamos a evaluar en nuestras aproximaciones
  
  t1.aprox <- aprox_1(p.eval)
  t2.aprox <- aprox_2(p.eval)
  
  dif <- abs(t2.aprox-t1.aprox)
  
  tabla <- data.frame(surv=p.eval, tiempo1=t1.aprox, tiempo2=t2.aprox, diferencia=dif)
  
  tabla <- tabla[!is.na(tabla$diferencia),]
  tabla <- tabla[tabla$surv>surv.min,]
  
  dif.max <- max(tabla$diferencia)
  
  grafica <- ggplot(tabla) + geom_line(aes(surv,tiempo1,colour="blue"))+
    geom_line(aes(surv,tiempo2,colour="red"))+ 
    ggtitle(paste("Distancia maxima:",round(dif.max,4)))+
    geom_vline(xintercept = tabla$surv[tabla$diferencia==dif.max], linetype="dashed", color = "black")+
    geom_hline(yintercept = tabla$tiempo1[tabla$diferencia==dif.max], linetype="dashed", color = "black")+
    geom_hline(yintercept = tabla$tiempo2[tabla$diferencia==dif.max], linetype="dashed", color = "black") +
    labs(x = "Probabilidad de supervivencia", y = "Tiempo",colour="Estratos")+coord_flip() 
  
  return(grafica)
  
}

A continuación, un ejemplo de la función

d.max(km.veteran$surv,km.veteran$time,km.veteran$strata,0.1)
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values