El objetivo de esta práctica es crear una función que nos haga gráficos diagnóstico para determinar si a una curva de supervivencia estimada con KM le subyace algún modelo paramétrico. Nuestra función recibe los siguientes parámetros:\ - df: Corresponde al DataFrame a estudiar con datos de supervivencia. - n: El indice de la columna que representa el tiempo en df\ - m: El indice de la columna que representa las fallas en df\ Así, se creará una curva de supervivencia estimada y posteriormente se comparará gráficamente con los modelos exponencial, gengamma, lognormal y gompertz. Para visualizar que tan buenas son nuestras comparaciones, creamos un gráfico de dispersión entre F la función de distribución según el modelo con el parametro de maxima verosimilitus y 1-S_{KM}, la curva estimada por el método de Kaplan Meier. Si la aproximación es buena, nos gustaría que este ultimo gráfico se vea como una identidad.

require(plotly)
require(survival)
require(KMsurv)
require(dplyr)
require(ggplot2)
require(DataExplorer)
require(gganimate)
require(flexsurv)
require(actuar)
require(survminer)
require(latex2exp)

datos<-read.csv("C:/Users/Sol/Downloads/HR_comma_sep.csv")
#Estructura de los datos

diag<-function(df,n,m){
  pm1<-flexsurvreg(Surv(df[,n], df[,m])~1,dist= "exp")
  pm2<-flexsurvreg(Surv(df[,n], df[,m])~1,dist= "lognormal")
  pm3<-flexsurvreg(Surv(df[,n], df[,m])~1, dist= "gompertz")
  pm4<-flexsurvreg(Surv(df[,n], df[,m])~1, dist= "gengamma")
  
  km_model<-surv_fit(Surv(df[,n], df[,m])~1, df, type="kaplan-meier")
  par(mfrow=c(2,2))
  plot(pm1,type = "survival", ci = F, main = "Ajuste exponencial", 
       conf.int = F, col = "#FF9900", col.obs = "black", xlab = "Tiempo", ylab = "Probabilidad de Supervivencia")
legend(1, .5, legend=c( "Kaplan Meier","exponencial"),
       col=c("black", "#FF9900"), lty=c(1,1), cex=0.8)

  plot(pm2, type = "survival", ci = F, main = "Ajuste Lognormal", 
       conf.int = F, col = "#99FF33", col.obs = "black", xlab = "Tiempo", ylab = "Probabilidad de Supervivencia")
  legend(1, .5, legend=c( "Kaplan Meier","Lognormal"),
       col=c("black", "#99FF33"), lty=c(1,1), cex=0.8)
  plot(pm3, type = "survival", ci = F, main = "Ajuste gompertz", 
       conf.int = F, col="#00FFFF", col.obs = "black", xlab = "Tiempo", ylab = "Probabilidad de Supervivencia")
 legend(1, .5, legend=c( "Kaplan Meier","Gompertz"),
       col=c("black", "#00FFFF"), lty=c(1,1), cex=0.8)
  plot(pm4, type = "survival", ci = F, main = "Ajuste gengamma", 
       conf.int = F, col = "#FF00FF", col.obs = "black", xlab = "Tiempo", ylab = "Probabilidad de Supervivencia")
  legend(1, .5, legend=c( "Kaplan Meier","Gengamma"),
       col=c("black", "#FF00FF"), lty=c(1,1), cex=0.8)
  print("A continuacion se muestran graficos de bondad de ajuste")
 
  F_est1 <- 1 - as.data.frame(summary(pm1))[, "est"]
  F_est2 <- 1 - as.data.frame(summary(pm2))[, "est"]
  F_est3 <- 1 - as.data.frame(summary(pm3))[, "est"]
  F_est4 <- 1 - as.data.frame(summary(pm4))[, "est"]
  KM_est <- 1 - km_model$surv
  par(mfrow=c(2,2))
  plot(F_est1, KM_est,col="#FF9900", main="Bondad de ajuste: Modelo exponencial", ylab = TeX('$1-S_{KM}$'), xlab=TeX('$F_{EMV}$'))
  plot(F_est2, KM_est,col="#99FF33", main="Bondad de ajuste: Modelo Log-normal", ylab = TeX('$1-S_{KM}$'), xlab=TeX('$F_{EMV}$'))
  plot(F_est3, KM_est,col="#00FFFF", main="Bondad de ajuste: Modelo gompertz", ylab = TeX('$1-S_{KM}$'), xlab=TeX('$F_{EMV}$'))
  plot(F_est4, KM_est,col="#FF00FF", main="Bondad de ajuste: Modelo gamma generalizada", ylab = TeX('$1-S_{KM}$'), xlab=TeX('$F_{EMV}$'))
}

Ejemplo 1. Dataframe que muestra en la columna 5 el tiempo gastado en la compania y 7 si el individuó la dejo o no.

diag(datos,4,6)

## [1] "A continuacion se muestran graficos de bondad de ajuste"

Ejemplo 2. Data frame que representa un experimento con ratas, tenemos que la columna 3 corresponde al tiempo que ocurrio la falla y la columna 6 corresponde a una columna binaria dependiendo si ocurrio el evento o no(muerte).

diag(rats,3,4)

## [1] "A continuacion se muestran graficos de bondad de ajuste"

Ejemplo 3. Data frame que estudia cancer en la vejiga, tenemos que la columna 5 corresponde al tiempo ocurrió la falla y la columna 6 corresponde a una columna binaria dependiendo si ocurrió el evento o no.

diag(bladder,5,6 )

## [1] "A continuacion se muestran graficos de bondad de ajuste"