El Mejor Ajuste Paramétrico a una Función de Supervivencia Estimada por K-M

Buscaremos el mejor ajuste para una función de supervivencia

Obs: En este caso sólo usaremos 3 distribuciones para comparar con la función de supervivencia: 1. Exponencial 2. Weibull 3. Log-Normal Sin embargo, existe una gama amplia de distribuciones que se pueden utilizar para encontrar el mejor ajuste. Lo primero que debemos hacer es obtener la función de supervivencia estimada de alguna base de datos que nos interese.

Iniciamos cargando las librerías

library(survival)
library(KMsurv)
library(survMisc)
library(survminer)
## Loading required package: ggplot2
## 
## Attaching package: 'ggplot2'
## The following object is masked from 'package:survMisc':
## 
##     autoplot
## Loading required package: ggpubr
## Loading required package: magrittr
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

Ahora carrgamos la base de datos

getwd()
## [1] "C:/Users/bec-rene3/Documents/Lizbeth/Trabajos/R"
datos <- read.csv(file = "churn.csv", header = T)
head(datos)
##   churn accountlength internationalplan voicemailplan numbervmailmessages
## 1    No           128                no           yes                  25
## 2    No           107                no           yes                  26
## 3    No           137                no            no                   0
## 4    No            84               yes            no                   0
## 5    No            75               yes            no                   0
## 6    No           118               yes            no                   0
##   totaldayminutes totaldaycalls totaldaycharge totaleveminutes
## 1           265.1           110          45.07           197.4
## 2           161.6           123          27.47           195.5
## 3           243.4           114          41.38           121.2
## 4           299.4            71          50.90            61.9
## 5           166.7           113          28.34           148.3
## 6           223.4            98          37.98           220.6
##   totalevecalls totalevecharge totalnightminutes totalnightcalls
## 1            99          16.78             244.7              91
## 2           103          16.62             254.4             103
## 3           110          10.30             162.6             104
## 4            88           5.26             196.9              89
## 5           122          12.61             186.9             121
## 6           101          18.75             203.9             118
##   totalnightcharge totalintlminutes totalintlcalls totalintlcharge
## 1            11.01             10.0              3            2.70
## 2            11.45             13.7              3            3.70
## 3             7.32             12.2              5            3.29
## 4             8.86              6.6              7            1.78
## 5             8.41             10.1              3            2.73
## 6             9.18              6.3              6            1.70
##   numbercustomerservicecalls
## 1                          1
## 2                          1
## 3                          0
## 4                          2
## 5                          3
## 6                          0

Recodificamos la variable que nos interesa tomar como censura

datos$churn <- ifelse(datos$churn == "Yes", 1, 0)

Ahora realizamos el ajuste de la función de supervivencia (para esta función sólo podemos analizar el mejor ajuste con un objeto del tipo surv_fit() sin estratificación), sin embargo se puede desarrollar una función que calcule la mejor aproximación para cada estrato [Tal vez un futuro proyecto :D ]

ajuste<- surv_fit(Surv(accountlength, churn) ~ 1, data = datos)

Creamos la función para obtener el mejor modelo parámétrico que se ajuste a nuestra función de supervivencia estimada.

mejorajuste<-function(ajuste){
  
Dist<-c("exp","weibull","lnorm")
a1<- lm(log(ajuste$surv)~ ajuste$time)
a2<-lm(log(-log(ajuste$surv))~log(ajuste$time))
a3<-lm(log(ajuste$time)~qnorm(1-ajuste$surv))
r1<-summary(a1)
r2<-summary(a2)
r3<-summary(a3)
vrsqr<-c(r1$adj.r.squared,r2$adj.r.squared,r3$adj.r.squared)
majs<-Dist[which.max(vrsqr)]
majs
plot(ajuste$time,log(ajuste$surv), main="Comparación con distribución Exponencial", col= "lightcoral")
abline (a1, lwd=1, col ="black" )
plot(log(ajuste$time),log(-log(ajuste$surv)), main="Comparación con distribución Weibull", col= "cyan4")
abline (a2, lwd=1, col ="black" )
plot(qnorm(1-ajuste$surv),log(ajuste$time), main="Comparación con distribución Log-normal", col= "limegreen")
abline (a3, lwd=1, col ="black" )
p<-print(c("La disribución que más se ajusta a esta función de supervivencia es",majs))
}

Ahora probemos a nuestra bebé:

mejorajuste(ajuste)

## [1] "La disribución que más se ajusta a esta función de supervivencia es"
## [2] "weibull"

como podemos observar la distribución que más se asemeja a nuestros datos es una Weibull, para tomar esta decisión se usó como criterio el valor de la R cuadrada resultante de aplicar una regresión lineal a nuestros datos transformados tanto del tiempo como de la función de supervivencia ajustada. Sin embargo se debe de tener cuidado al momento de la elección.

Bueno, eso es todo amiguitos :)