Introducción

Existen funciones en algunas paqueterías de R que ajustan observaciones empíricas a un modelo de distribución paramétrica. Pero la mayoría de las veces (si no es que todas) no tenemos ni una pista a qué función paramétrica se ajuste más nuestro modelo. Para este caso se creó la función pparam que nos da una noción gráfica de a qué posibles distribuciones paramétricas se ajuste.

pparam

Esta función compara la relación entre función de supervivencia y el tiempo del modelo que se esté trabajando de manera equivalente a como se comportaría si fuera parte de los modelos: “Exponencial”, “Weibull”, “Gompertz”, “Lognormal” y “Log-logística”.

Forma de uso

pparam(vsurv, vtime)

Argumentos

El modelo debió haber sido ajustado previamente. De manera preferencial con las funciones surv de la paquetería survival y surv_fit de la paquetería KMsurv. Los argumentos son:

  • vsurv vector de supervivencia S(t).
  • vtime vector de tiempo t.

ambos del modelo ajustado, no usar directo de la base de datos.

Código

pparam <- function(vsurv,vtime){
  #tratamiento previo
  par(mfrow=c(3,2)) #divide la gráfica en una matriz de 3 filas y 2 columnas
  #Exponencial
  plot(log(vsurv),vtime, main = "Diagnóstico Exponencial", xlab = expression(hat(S)(t)), ylab="t")
  mtext("Se espera una línea recta con pendiente negativa", side=3, adj = 0, cex = 0.8) #info
  mtext("Hecho por: Itzel Escutia", side=1, line=4, adj = 1, cex=0.5) #etiqueta
  grid()
  #Weibull
  plot(log(-log(vsurv)),log(vtime), main = "Diagnóstico Weibull", xlab = expression(hat(S)(t)), ylab="t")
  mtext("Se espera una línea recta con pendiente positiva", side=3, adj = 0, cex = 0.8) #info
  mtext("Hecho por: Itzel Escutia", side=1, line=4, adj = 1, cex=0.5) #etiqueta
  grid()
  #Gompertz
  plot(log(1-log(vsurv)),vtime, main = "Diagnóstico Gompertz", xlab = expression(hat(S)(t)), ylab="t")
  mtext("Se espera una línea recta con pendiente positiva", side=3, adj = 0, cex = 0.8) #info
  mtext("Hecho por: Itzel Escutia", side=1, line=4, adj = 1, cex=0.5) #etiqueta
  grid()
  #Lognormal
  plot(qnorm(1-vsurv),log(vtime), main = "Diagnóstico Lognormal", xlab = expression(hat(S)(t)), ylab="t")
  mtext("Se espera una línea recta con pendiente positiva", side=3, adj = 0, cex = 0.8) #info
  mtext("Hecho por: Itzel Escutia", side=1, line=4, adj = 1, cex=0.5) #etiqueta
  grid()
  #Log-logistic
  plot(log(1/vsurv - 1),log(vtime), main = "Diagnóstico Log-logistic", xlab = expression(hat(S)(t)), ylab="t")
  mtext("Se espera una línea recta con pendiente positiva", side=3, adj = 0, cex = 0.8) #info
  mtext("Hecho por: Itzel Escutia", side=1, line=4, adj = 1, cex=0.5) #etiqueta
  grid()
}

Básicamente no se necesita ninguna paquetería para su uso.

Ejemplos

Para los ejemplos se requiere usar las siguientes paqueterías

require(survival)
require(flexsurv)
require(KMsurv)

1. 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. Sea datos el data frame.

datos <- read.csv("HR_comma_sep.csv") #se lee la base de datos desde un archivo
km_model <- surv_fit(Surv(time_spend_company,left)~1, datos) #Se ajusta el modelo de acuerdo al data frame
pparam(km_model$surv,km_model$time) #se usa la función

Entre el título y la gráfica se observa cual es el criterio de comparación para decir si se ajusta nuestro modelo a la respectiva distribución paramétrica. En este caso no se ve ninguna recta, por lo que no se ajusta a ninguna de las distribuciones mencionadas.

2. Pulmón

Se tiene la base de datos lung incluida en la paqueteria survival y hace referencia a la supervivencia de pacientes con cáncer pulmonar avanzado provenientes del North Central Cancer Treatment Group. De sus variables nos fijamos en el tiempo (time) y el estatus (status), para hacer el modelo de supervivencia. Para esta base de datos el estatus es 1 si es censurado y 2 si murió, por lo que se tienen que convertir a unos y ceros respectivamente para ajustarlos. Sea pulmon el data frame.

pulmon <- lung %>% mutate(status = status-1) #se guarda el data frame con el nuevo status en "pulmon"
km_pulmon <- surv_fit(Surv(time,status)~1,pulmon) #Se hace el modelo de supervivencia
pparam(km_pulmon$surv,km_pulmon$time)

En este se puede observar que se forma casi una recta de pendiente positiva cuando se compara como la distribución paramétrica de Gompertz, por lo que podemos especular que el tiempo se distribuye de esa manera.

3. 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).

km_kidney <- surv_fit(Surv(time,status)~1,kidney)
pparam(km_kidney$surv,km_kidney$time)

Esta muestra tampoco da indicios de distribuirse en alguna de las familias paramétricas presentadas en la función.