#### DIGITALIZACI?N ####

GUYOT <- function(ensayo,
                  farmaco,
                  brazo,
                  triski_fco,
                  Nriski_fco,
                  total.events_fco,
                  totalpts,
                  datos_digitalizados,
                  cortes_FP) #if doing a range lookup, should the smaller or larger key be used?)
{
  
  fco_csv <- read_delim(datos_digitalizados, ";", escape_double = FALSE, col_names = FALSE,
                        locale = locale(decimal_mark = "."), trim_ws = TRUE, col_types = cols())
  fco_csv <- data.frame(fco_csv)
  trisk <- triski_fco
  nrisk <- Nriski_fco
  pre <- preprocess(dat=fco_csv, trisk=trisk, nrisk=nrisk,maxy=1)
  est <- getIPD(prep=pre,armID=brazo,tot.events=total.events_fco)
  KM.est_fco <- survfit(Surv(time, status) ~ treat, data=est$IPD, type="kaplan-meier")
  tiempos <- seq(from=0, to=max(trisk), by=cortes_FP)
  summaryKM_fco <- summary(KM.est_fco, tiempos)
  write.xlsx(est$IPD, file=paste0("./RESULTADOS/IPD_",farmaco,".xlsx"))
  output <- list("ENSAYO" = ensayo, "FARMACO"= farmaco, "BRAZO"=brazo, "KAPLAN" = KM.est_fco, "IPD"=est$IPD)
  return(output)}


#### COX ####

COX2 <- function(fco1, fco2) #if doing a range lookup, should the smaller or larger key be used?)
{
    datos_fco1 <- fco1
    datos_fco2 <- fco2
    ####
    data_fco <- rbind(datos_fco1$IPD, datos_fco2$IPD)
    KM.est_combinado <- survfit(Surv(time, status) ~ treat, data=data_fco, type="kaplan-meier")
    data.survdiff <- survdiff(Surv(time, status) ~ treat, data=data_fco)
    p.val = 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1)
    HR = (data.survdiff$obs[2]/data.survdiff$exp[2])/(data.survdiff$obs[1]/data.survdiff$exp[1])
    up95 = exp(log(HR) + qnorm(0.975)*sqrt(1/data.survdiff$exp[2]+1/data.survdiff$exp[1]))
    low95 = exp(log(HR) - qnorm(0.975)*sqrt(1/data.survdiff$exp[2]+1/data.survdiff$exp[1]))
    resultados_cox_data_1 <- cbind(1/HR,1/up95,1/low95)
    colnames(resultados_cox_data_1) <- c("HR", "IC95MIN", "IC95MAX")
    ####
    #### GRAFICO
    data_fco_2 <<- rbind(datos_fco1$IPD, datos_fco2$IPD)
    KM.est_combinado_2 <<- survfit(Surv(time, status) ~ treat, data=data_fco_2, type="kaplan-meier")
    grafico <- ggsurvplot(
      KM.est_combinado_2,
      size = 1,                 # change line size
      palette =
        c("#E7B800", "#2E9FDF"),# custom color palettes
      conf.int = TRUE,          # Add confidence interval
      pval = FALSE,              # Add p-value
      risk.table = TRUE,        # Add risk table
      risk.table.col = "strata",# Risk table color by groups
      legend.labs =
        c(datos_fco1$FARMACO, datos_fco2$FARMACO),    # Change legend labels
      risk.table.height = 0.25, # Useful to change when you have multiple groups
      ggtheme = theme_bw(),
      break.time.by = 3# Change ggplot2 theme
    )
    #### OUTPUT
    output <-list("KAPLAN"=KM.est_combinado,
                  "COX_1"=resultados_cox_data_1,
                  "PLOT"=grafico)
    names(output)[2] <- paste("HR", fco1$FARMACO, "vs", fco2$FARMACO, sep=" ")
    return(output)}


COX3 <- function(fco1, fco2, fco3) #if doing a range lookup, should the smaller or larger key be used?)
{
    datos_fco1 <- fco1
    datos_fco2 <- fco2
    datos_fco3 <- fco3
    ###
    data_fco <- rbind(datos_fco1$IPD, datos_fco2$IPD)
    KM.est_combinado <- survfit(Surv(time, status) ~ treat, data=data_fco, type="kaplan-meier")
    data.survdiff <- survdiff(Surv(time, status) ~ treat, data=data_fco)
    p.val = 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1)
    HR = (data.survdiff$obs[2]/data.survdiff$exp[2])/(data.survdiff$obs[1]/data.survdiff$exp[1])
    up95 = exp(log(HR) + qnorm(0.975)*sqrt(1/data.survdiff$exp[2]+1/data.survdiff$exp[1]))
    low95 = exp(log(HR) - qnorm(0.975)*sqrt(1/data.survdiff$exp[2]+1/data.survdiff$exp[1]))
    resultados_cox_data_1 <- cbind(1/HR,1/up95,1/low95)
    colnames(resultados_cox_data_1) <- c("HR", "IC95MIN", "IC95MAX")
    ###
    data_fco <- rbind(datos_fco1$IPD, datos_fco3$IPD)
    KM.est_combinado <- survfit(Surv(time, status) ~ treat, data=data_fco, type="kaplan-meier")
    data.survdiff <- survdiff(Surv(time, status) ~ treat, data=data_fco)
    p.val = 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1)
    HR = (data.survdiff$obs[2]/data.survdiff$exp[2])/(data.survdiff$obs[1]/data.survdiff$exp[1])
    up95 = exp(log(HR) + qnorm(0.975)*sqrt(1/data.survdiff$exp[2]+1/data.survdiff$exp[1]))
    low95 = exp(log(HR) - qnorm(0.975)*sqrt(1/data.survdiff$exp[2]+1/data.survdiff$exp[1]))
    resultados_cox_data_2 <- cbind(1/HR,1/up95,1/low95)
    colnames(resultados_cox_data_2) <- c("HR", "IC95MIN", "IC95MAX")
    ###
    data_fco <- rbind(datos_fco2$IPD, datos_fco3$IPD)
    KM.est_combinado <- survfit(Surv(time, status) ~ treat, data=data_fco, type="kaplan-meier")
    data.survdiff <- survdiff(Surv(time, status) ~ treat, data=data_fco)
    p.val = 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1)
    HR = (data.survdiff$obs[2]/data.survdiff$exp[2])/(data.survdiff$obs[1]/data.survdiff$exp[1])
    up95 = exp(log(HR) + qnorm(0.975)*sqrt(1/data.survdiff$exp[2]+1/data.survdiff$exp[1]))
    low95 = exp(log(HR) - qnorm(0.975)*sqrt(1/data.survdiff$exp[2]+1/data.survdiff$exp[1]))
    resultados_cox_data_3 <- cbind(1/HR,1/up95,1/low95)
    colnames(resultados_cox_data_3) <- c("HR", "IC95MIN", "IC95MAX")
    #### GRAFICO
    data_fco_3 <<- rbind(datos_fco1$IPD, datos_fco2$IPD, datos_fco3$IPD)
    KM.est_combinado_3 <<- survfit(Surv(time, status) ~ treat, data=data_fco_3, type="kaplan-meier")
    grafico <- ggsurvplot(
      KM.est_combinado_3,
      size = 1,                 # change line size
      palette =
        c("#E7B800", "#2E9FDF","red" ),# custom color palettes
      conf.int = TRUE,          # Add confidence interval
      pval = FALSE,              # Add p-value
      risk.table = TRUE,        # Add risk table
      risk.table.col = "strata",# Risk table color by groups
      legend.labs =
        c(datos_fco1$FARMACO, datos_fco2$FARMACO, datos_fco3$FARMACO),    # Change legend labels
      risk.table.height = 0.25, # Useful to change when you have multiple groups
      ggtheme = theme_bw(),
      break.time.by = 3# Change ggplot2 theme
    )
    #### OUTPUT
    data_fco <- rbind(datos_fco1$IPD, datos_fco2$IPD, datos_fco3$IPD)
    KM.est_combinado <- survfit(Surv(time, status) ~ treat, data=data_fco, type="kaplan-meier")
    output <-list("KAPLAN"=KM.est_combinado,
                  "COX_1"=resultados_cox_data_1,
                  "COX_2"=resultados_cox_data_2,
                  "COX_3"=resultados_cox_data_3,
                  "PLOT"=grafico)
    names(output)[2] <- paste("HR", fco1$FARMACO, "vs", fco2$FARMACO, sep=" ")
    names(output)[3] <- paste("HR", fco1$FARMACO, "vs", fco3$FARMACO, sep=" ")
    names(output)[4] <- paste("HR", fco2$FARMACO, "vs", fco3$FARMACO, sep=" ")
    return(output)}

COX <- function(fco1,
                fco2,
                fco3="NULO"){
  datos_fco1 <- fco1
  datos_fco2 <- fco2
  datos_fco3 <- fco3
  if(datos_fco3 == "NULO") {
    ####
    data_fco <- rbind(datos_fco1$IPD, datos_fco2$IPD)
    KM.est_combinado <- survfit(Surv(time, status) ~ treat, data=data_fco, type="kaplan-meier")
    data.survdiff <- survdiff(Surv(time, status) ~ treat, data=data_fco)
    p.val = 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1)
    HR = (data.survdiff$obs[2]/data.survdiff$exp[2])/(data.survdiff$obs[1]/data.survdiff$exp[1])
    up95 = exp(log(HR) + qnorm(0.975)*sqrt(1/data.survdiff$exp[2]+1/data.survdiff$exp[1]))
    low95 = exp(log(HR) - qnorm(0.975)*sqrt(1/data.survdiff$exp[2]+1/data.survdiff$exp[1]))
    resultados_cox_data_1 <- cbind(1/HR,1/up95,1/low95)
    colnames(resultados_cox_data_1) <- c("HR", "IC95MIN", "IC95MAX")
    ####
    #### GRAFICO
    data_fco_2 <<- rbind(datos_fco1$IPD, datos_fco2$IPD)
    KM.est_combinado_2 <<- survfit(Surv(time, status) ~ treat, data=data_fco_2, type="kaplan-meier")
    grafico <- ggsurvplot(
      KM.est_combinado_2,
      size = 1,                 # change line size
      palette =
        c("#E7B800", "#2E9FDF"),# custom color palettes
      conf.int = TRUE,          # Add confidence interval
      pval = FALSE,              # Add p-value
      risk.table = TRUE,        # Add risk table
      risk.table.col = "strata",# Risk table color by groups
      legend.labs =
        c(datos_fco1$FARMACO, datos_fco2$FARMACO),    # Change legend labels
      risk.table.height = 0.25, # Useful to change when you have multiple groups
      ggtheme = theme_bw(),
      break.time.by = 3# Change ggplot2 theme
    )
    #### OUTPUT
    output <-list("KAPLAN"=KM.est_combinado,
                  "COX_1"=resultados_cox_data_1,
                  "PLOT"=grafico)
    names(output)[2] <- paste("HR", fco1$FARMACO, "vs", fco2$FARMACO, sep=" ")
    return(output)
  }
  else {###
    data_fco <- rbind(datos_fco1$IPD, datos_fco2$IPD)
    KM.est_combinado <- survfit(Surv(time, status) ~ treat, data=data_fco, type="kaplan-meier")
    data.survdiff <- survdiff(Surv(time, status) ~ treat, data=data_fco)
    p.val = 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1)
    HR = (data.survdiff$obs[2]/data.survdiff$exp[2])/(data.survdiff$obs[1]/data.survdiff$exp[1])
    up95 = exp(log(HR) + qnorm(0.975)*sqrt(1/data.survdiff$exp[2]+1/data.survdiff$exp[1]))
    low95 = exp(log(HR) - qnorm(0.975)*sqrt(1/data.survdiff$exp[2]+1/data.survdiff$exp[1]))
    resultados_cox_data_1 <- cbind(1/HR,1/up95,1/low95)
    colnames(resultados_cox_data_1) <- c("HR", "IC95MIN", "IC95MAX")
    ###
    data_fco <- rbind(datos_fco1$IPD, datos_fco3$IPD)
    KM.est_combinado <- survfit(Surv(time, status) ~ treat, data=data_fco, type="kaplan-meier")
    data.survdiff <- survdiff(Surv(time, status) ~ treat, data=data_fco)
    p.val = 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1)
    HR = (data.survdiff$obs[2]/data.survdiff$exp[2])/(data.survdiff$obs[1]/data.survdiff$exp[1])
    up95 = exp(log(HR) + qnorm(0.975)*sqrt(1/data.survdiff$exp[2]+1/data.survdiff$exp[1]))
    low95 = exp(log(HR) - qnorm(0.975)*sqrt(1/data.survdiff$exp[2]+1/data.survdiff$exp[1]))
    resultados_cox_data_2 <- cbind(1/HR,1/up95,1/low95)
    colnames(resultados_cox_data_2) <- c("HR", "IC95MIN", "IC95MAX")
    ###
    data_fco <- rbind(datos_fco2$IPD, datos_fco3$IPD)
    KM.est_combinado <- survfit(Surv(time, status) ~ treat, data=data_fco, type="kaplan-meier")
    data.survdiff <- survdiff(Surv(time, status) ~ treat, data=data_fco)
    p.val = 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1)
    HR = (data.survdiff$obs[2]/data.survdiff$exp[2])/(data.survdiff$obs[1]/data.survdiff$exp[1])
    up95 = exp(log(HR) + qnorm(0.975)*sqrt(1/data.survdiff$exp[2]+1/data.survdiff$exp[1]))
    low95 = exp(log(HR) - qnorm(0.975)*sqrt(1/data.survdiff$exp[2]+1/data.survdiff$exp[1]))
    resultados_cox_data_3 <- cbind(1/HR,1/up95,1/low95)
    colnames(resultados_cox_data_3) <- c("HR", "IC95MIN", "IC95MAX")
    #### GRAFICO
    data_fco_3 <<- rbind(datos_fco1$IPD, datos_fco2$IPD, datos_fco3$IPD)
    KM.est_combinado_3 <<- survfit(Surv(time, status) ~ treat, data=data_fco_3, type="kaplan-meier")
    grafico <- ggsurvplot(
      KM.est_combinado_3,
      size = 1,                 # change line size
      palette =
        c("#E7B800", "#2E9FDF","red" ),# custom color palettes
      conf.int = TRUE,          # Add confidence interval
      pval = FALSE,              # Add p-value
      risk.table = TRUE,        # Add risk table
      risk.table.col = "strata",# Risk table color by groups
      legend.labs =
        c(datos_fco1$FARMACO, datos_fco2$FARMACO, datos_fco3$FARMACO),    # Change legend labels
      risk.table.height = 0.25, # Useful to change when you have multiple groups
      ggtheme = theme_bw(),
      break.time.by = 3# Change ggplot2 theme
    )
    #### OUTPUT
    data_fco <- rbind(datos_fco1$IPD, datos_fco2$IPD, datos_fco3$IPD)
    KM.est_combinado <- survfit(Surv(time, status) ~ treat, data=data_fco, type="kaplan-meier")
    output <-list("KAPLAN"=KM.est_combinado,
                  "COX_1"=resultados_cox_data_1,
                  "COX_2"=resultados_cox_data_2,
                  "COX_3"=resultados_cox_data_3,
                  "PLOT"=grafico)
    names(output)[2] <- paste("HR", fco1$FARMACO, "vs", fco2$FARMACO, sep=" ")
    names(output)[3] <- paste("HR", fco1$FARMACO, "vs", fco3$FARMACO, sep=" ")
    names(output)[4] <- paste("HR", fco2$FARMACO, "vs", fco3$FARMACO, sep=" ")
    return(output)}
}

#### MODELIZACI?N ####
# method = c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN","Nelder-Mead")

MODELIZACION <- function(fco)
{
datos_fco <- fco
KAPLAN <- datos_fco$KAPLAN
IPD <- datos_fco$IPD
FARMACO <- datos_fco$FARMACO
sGamma_fco <- tryCatch(expr =  flexsurvreg(formula = Surv(time, status) ~ 1, data = datos_fco$IPD,  dist="gamma"),
                        error = function(e){flexsurvreg(formula = Surv(time, status) ~ 1, data = datos_fco$IPD,  dist="gamma",
                                                        control=list(fnscale = 2500), method="SANN")})
sLlogis_fco <- tryCatch(expr =  flexsurvreg(formula = Surv(time, status) ~ 1, data = datos_fco$IPD,  dist="llogis"),
                        error = function(e){flexsurvreg(formula = Surv(time, status) ~ 1, data = datos_fco$IPD,  dist="llogis",
                                                        control=list(fnscale = 2500), method="SANN")})
sWei_fco <- tryCatch(expr =  flexsurvreg(formula = Surv(time, status) ~ 1, data = datos_fco$IPD,  dist="weibull"),
                     error = function(e){flexsurvreg(formula = Surv(time, status) ~ 1, data = datos_fco$IPD,  dist="weibull",
                                                     control=list(fnscale = 2500), method="SANN")})
sLn_fco <- tryCatch(expr =  flexsurvreg(formula = Surv(time, status) ~ 1, data = datos_fco$IPD,  dist="lnorm"),
                    error = function(e){flexsurvreg(formula = Surv(time, status) ~ 1, data = datos_fco$IPD,  dist="lnorm",
                                                   control=list(fnscale = 2500), method="SANN")})
sGo_fco <- tryCatch(expr =  flexsurvreg(formula = Surv(time, status) ~ 1, data = datos_fco$IPD,  dist="gompertz"),
                    error = function(e){flexsurvreg(formula = Surv(time, status) ~ 1, data = datos_fco$IPD,  dist="gompertz",
                                                    control=list(fnscale = 2500), method="SANN")})
sGen_fco <- tryCatch(expr =  flexsurvreg(formula = Surv(time, status) ~ 1, data = datos_fco$IPD,  dist="gengamma"),
                     error = function(e){flexsurvreg(formula = Surv(time, status) ~ 1, data = datos_fco$IPD,  dist="gengamma",
                                                     control=list(fnscale = 2500), method="SANN")})
sGf_fco <- tryCatch(expr =  flexsurvreg(formula = Surv(time, status) ~ 1, data = datos_fco$IPD,  dist="genf"),
                    error = function(e){flexsurvreg(formula = Surv(time, status) ~ 1, data = datos_fco$IPD,  dist="genf",
                                                    control=list(fnscale = 2500), method="SANN")})
sExp_fco <- tryCatch(expr =  flexsurvreg(formula = Surv(time, status) ~ 1, data = datos_fco$IPD,  dist="exponential"),
                     error = function(e){flexsurvreg(formula = Surv(time, status) ~ 1, data = datos_fco$IPD,  dist="exponential",
                                                     control=list(fnscale = 2500), method="SANN")})
# spline1_hazard_fco <- tryCatch(expr = flexsurvspline(Surv(time, status) ~ 1, data=datos_fco$IPD, k=1, scale="hazard"),
#                                error = function(e){flexsurvspline(Surv(time, status) ~ 1, data=datos_fco$IPD, k=1, scale="hazard", 
#                                                                   control=list(fnscale = 2500), method="SANN")})
# spline2_hazard_fco <- tryCatch(expr = flexsurvspline(Surv(time, status) ~ 1, data=datos_fco$IPD, k=2, scale="hazard"),
#                                error = function(e){flexsurvspline(Surv(time, status) ~ 1, data=datos_fco$IPD, k=2, scale="hazard", 
#                                                                   control=list(fnscale = 2500), method="SANN")})
# spline3_hazard_fco <- tryCatch(expr = flexsurvspline(Surv(time, status) ~ 1, data=datos_fco$IPD, k=3, scale="hazard"),
#                                error = function(e){flexsurvspline(Surv(time, status) ~ 1, data=datos_fco$IPD, k=3, scale="hazard", 
#                                                                   control=list(fnscale = 2500), method="SANN")})
# spline4_hazard_fco <- tryCatch(expr = flexsurvspline(Surv(time, status) ~ 1, data=datos_fco$IPD, k=4, scale="hazard"),
#                                error = function(e){flexsurvspline(Surv(time, status) ~ 1, data=datos_fco$IPD, k=4, scale="hazard", 
#                                                                   control=list(fnscale = 2500), method="SANN")})
# spline5_hazard_fco <- tryCatch(expr = flexsurvspline(Surv(time, status) ~ 1, data=datos_fco$IPD, k=5, scale="hazard"),
#                                error = function(e){flexsurvspline(Surv(time, status) ~ 1, data=datos_fco$IPD, k=5, scale="hazard", 
#                                                                   control=list(fnscale = 2500), method="SANN")})
# spline6_hazard_fco <- tryCatch(expr = flexsurvspline(Surv(time, status) ~ 1, data=datos_fco$IPD, k=6, scale="hazard"),
#                                error = function(e){flexsurvspline(Surv(time, status) ~ 1, data=datos_fco$IPD, k=6, scale="hazard", 
#                                                                   control=list(fnscale = 2500), method="SANN")})
# spline1_odds_fco <- tryCatch(expr = flexsurvspline(Surv(time, status) ~ 1, data=datos_fco$IPD, k=1, scale="odds"),
#                              error = function(e){flexsurvspline(Surv(time, status) ~ 1, data=datos_fco$IPD, k=1, scale="odds", 
#                                                                 control=list(fnscale = 2500), method="SANN")})
# spline2_odds_fco <- tryCatch(expr = flexsurvspline(Surv(time, status) ~ 1, data=datos_fco$IPD, k=2, scale="odds"),
#                              error = function(e){flexsurvspline(Surv(time, status) ~ 1, data=datos_fco$IPD, k=2, scale="odds", 
#                                                                 control=list(fnscale = 2500), method="SANN")})
# spline3_odds_fco <- tryCatch(expr = flexsurvspline(Surv(time, status) ~ 1, data=datos_fco$IPD, k=3, scale="odds"),
#                              error = function(e){flexsurvspline(Surv(time, status) ~ 1, data=datos_fco$IPD, k=3, scale="odds", 
#                                                                 control=list(fnscale = 2500), method="SANN")})
# spline4_odds_fco <- tryCatch(expr = flexsurvspline(Surv(time, status) ~ 1, data=datos_fco$IPD, k=4, scale="odds"),
#                              error = function(e){flexsurvspline(Surv(time, status) ~ 1, data=datos_fco$IPD, k=4, scale="odds", 
#                                                                 control=list(fnscale = 2500), method="SANN")})
# spline5_odds_fco <- tryCatch(expr = flexsurvspline(Surv(time, status) ~ 1, data=datos_fco$IPD, k=5, scale="odds"),
#                              error = function(e){flexsurvspline(Surv(time, status) ~ 1, data=datos_fco$IPD, k=5, scale="odds", 
#                                                                 control=list(fnscale = 2500), method="SANN")})
# spline6_odds_fco <- tryCatch(expr = flexsurvspline(Surv(time, status) ~ 1, data=datos_fco$IPD, k=6, scale="odds"),
#                              error = function(e){flexsurvspline(Surv(time, status) ~ 1, data=datos_fco$IPD, k=6, scale="odds", 
#                                                                 control=list(fnscale = 2500), method="SANN")})
# spline1_normal_fco <- tryCatch(expr = flexsurvspline(Surv(time, status) ~ 1, data=datos_fco$IPD, k=1, scale="normal"),
#                                error = function(e){flexsurvspline(Surv(time, status) ~ 1, data=datos_fco$IPD, k=1, scale="normal", 
#                                                                   control=list(fnscale = 2500), method="SANN")})
# spline2_normal_fco <- tryCatch(expr = flexsurvspline(Surv(time, status) ~ 1, data=datos_fco$IPD, k=2, scale="normal"),
#                                error = function(e){flexsurvspline(Surv(time, status) ~ 1, data=datos_fco$IPD, k=2, scale="normal", 
#                                                                   control=list(fnscale = 2500), method="SANN")})
# spline3_normal_fco <- tryCatch(expr = flexsurvspline(Surv(time, status) ~ 1, data=datos_fco$IPD, k=3, scale="normal"),
#                                error = function(e){flexsurvspline(Surv(time, status) ~ 1, data=datos_fco$IPD, k=3, scale="normal", 
#                                                                   control=list(fnscale = 2500), method="SANN")})
# spline4_normal_fco <- tryCatch(expr = flexsurvspline(Surv(time, status) ~ 1, data=datos_fco$IPD, k=4, scale="normal"),
#                                error = function(e){flexsurvspline(Surv(time, status) ~ 1, data=datos_fco$IPD, k=4, scale="normal", 
#                                                                   control=list(fnscale = 2500), method="SANN")})
# spline5_normal_fco <- tryCatch(expr = flexsurvspline(Surv(time, status) ~ 1, data=datos_fco$IPD, k=5, scale="normal"),
#                                error = function(e){flexsurvspline(Surv(time, status) ~ 1, data=datos_fco$IPD, k=5, scale="normal", 
#                                                                   control=list(fnscale = 2500), method="SANN")})
# spline6_normal_fco <- tryCatch(expr = flexsurvspline(Surv(time, status) ~ 1, data=datos_fco$IPD, k=6, scale="normal"),
#                                error = function(e){flexsurvspline(Surv(time, status) ~ 1, data=datos_fco$IPD, k=6, scale="normal", 
#                                                       control=list(fnscale = 2500), method="SANN")})
AIC <- rbind.data.frame(
  cbind("Exponential", round(sExp_fco$AIC,1)),
  cbind("Gamma", round(sGamma_fco$AIC,1)),
  cbind("Lognormal", round(sLn_fco$AIC,1)),
  cbind("Loglogistic", round(sLlogis_fco$AIC,1)),
  cbind("Weibull", round(sWei_fco$AIC,1)),
  cbind("Gompertz", round(sGo_fco$AIC,1)),
  cbind("Generalized_gamma", round(sGen_fco$AIC,1)),
  #cbind("Lognormal", round(sLn_fco$AIC,1)),
  cbind("Generalized_F", round(sGf_fco$AIC,1))
  # cbind("spline1_hazard", round(spline1_hazard_fco$AIC,1)),
  # cbind("spline2_hazard", round(spline2_hazard_fco$AIC,1)),
  # cbind("spline3_hazard", round(spline3_hazard_fco$AIC,1)),
  # cbind("spline4_hazard", round(spline4_hazard_fco$AIC,1)),
  # cbind("spline5_hazard", round(spline5_hazard_fco$AIC,1)),
  # cbind("spline6_hazard", round(spline6_hazard_fco$AIC,1)),
  # cbind("spline1_odds", round(spline1_odds_fco$AIC,1)),
  # cbind("spline2_odds", round(spline2_odds_fco$AIC,1)),
  # cbind("spline3_odds", round(spline3_odds_fco$AIC,1)),
  # cbind("spline4_odds", round(spline4_odds_fco$AIC,1)),
  # cbind("spline5_odds", round(spline5_odds_fco$AIC,1)),
  # cbind("spline6_odds", round(spline6_odds_fco$AIC,1)),
  # cbind("spline1_normal", round(spline1_normal_fco$AIC,1)),
  # cbind("spline2_normal", round(spline2_normal_fco$AIC,1)),
  # cbind("spline3_normal", round(spline3_normal_fco$AIC,1)),
  # cbind("spline4_normal", round(spline4_normal_fco$AIC,1)),
  # cbind("spline5_normal", round(spline5_normal_fco$AIC,1)),
  # cbind("spline6_normal", round(spline6_normal_fco$AIC,1))
)
colnames(AIC) <- c("Modelo","AIC")
AIC <- AIC[order(AIC$AIC),]
output <- list(FARMACO = FARMACO,
               KAPLAN = KAPLAN,
               IPD = IPD,
               modelos = list("Exponential" = sExp_fco,
                              "Gamma" = sGamma_fco,
                              "Loglogistic"= sLlogis_fco,
                              "Weibull"=sWei_fco,
                              "Lognormal" = sLn_fco,
                              "Gompertz" = sGo_fco,
                              "Generalized_gamma" = sGen_fco,
                              "Generalized_F" = sGf_fco),
                              # "spline1_hazard" = spline1_hazard_fco,
                              # "spline2_hazard" = spline2_hazard_fco,
                              # "spline3_hazard" = spline3_hazard_fco,
                              # "spline4_hazard" = spline4_hazard_fco,
                              # "spline5_hazard" = spline5_hazard_fco,
                              # "spline6_hazard" = spline6_hazard_fco,
                              # "spline1_odds" = spline1_odds_fco,
                              # "spline2_odds" = spline2_odds_fco,
                              # "spline3_odds" = spline3_odds_fco,
                              # "spline4_odds" = spline4_odds_fco,
                              # "spline5_odds" = spline5_odds_fco,
                              # "spline6_odds" = spline6_odds_fco,
                              # "spline1_normal" = spline1_normal_fco,
                              # "spline2_normal" = spline2_normal_fco,
                              # "spline3_normal" = spline3_normal_fco,
                              # "spline4_normal" = spline4_normal_fco,
                              # "spline5_normal" = spline5_normal_fco,
                              # "spline6_normal" = spline6_normal_fco),
              NOMBRE_MODELO = c("Exponential",
                                "Gamma",
                                "Loglogistic",
                                   "Weibull",
                                   "Lognormal",
                                   "Gompertz",
                                   "Generalized_gamma",
                                   "Generalized_F"),
                                   # "spline1_hazard",
                                   # "spline2_hazard",
                                   # "spline3_hazard",
                                   # "spline4_hazard",
                                   # "spline5_hazard",
                                   # "spline6_hazard",
                                   # "spline1_odds",
                                   # "spline2_odds",
                                   # "spline3_odds",
                                   # "spline4_odds",
                                   # "spline5_odds",
                                   # "spline6_odds",
                                   # "spline1_normal",
                                   # "spline2_normal",
                                   # "spline3_normal",
                                   # "spline4_normal",
                                   # "spline5_normal",
                                   # "spline6_normal"),
              AIC=AIC)
return(output)}

#### CURVAS ####

CURVAS <- function(fco_mod,
                   tipo = "PFS",
                   tmax = 100,
                   cortes = 1)
{
  puntos <- c()
  for (i in 1:length(fco_mod$modelos)) {
    add <- data.frame(summary(fco_mod$modelos[[i]], type="survival",t=seq(0,tmax,by=cortes)))
    puntos <- c(puntos, list(add))
    names(puntos)[i] <- fco_mod$NOMBRE_MODELO[[i]]
  }
  wb_fco <- createWorkbook()
  for (i in 1:length(fco_mod$modelos)) {
    addWorksheet(wb_fco, fco_mod$NOMBRE_MODELO[[i]])
    puntos_curva <- data.frame(summary(fco_mod$modelos[[i]], type="survival",t=seq(0,tmax,by=cortes)))
    writeData(wb_fco, sheet = fco_mod$NOMBRE_MODELO[[i]], x = puntos_curva)
    saveWorkbook(wb_fco, file=paste0("./RESULTADOS/curvas_",tipo,"_",fco_mod$FARMACO,".xlsx"), overwrite = TRUE)
  }
  for (i in 1:length(fco_mod$modelos)) {
    PLOT_KAPLAN <- ggsurvplot(fco_mod$KAPLAN, data=fco_mod$IPD, xlim = c(0, tmax),
                              break.x.by = cortes, palette= "#E7B800")
    plot <- PLOT_KAPLAN$plot +
      geom_line(data = puntos[[i]], aes(x=time, y=est), colour = "#2E9FDF", size =0.5, alpha=0.8) +
      geom_line(data = puntos[[i]], aes(x=time, y=ucl), colour = "#2E9FDF", size =0.2, linetype="dotted")+
      geom_line(data = puntos[[i]], aes(x=time, y=lcl), colour = "#2E9FDF", size =0.2,
                linetype="dotted")+
      theme(legend.position = "none")+
      labs(title = fco_mod$NOMBRE_MODELO[[i]],
           subtitle = round(fco_mod$modelos[[i]]$AIC,2))
    tryCatch(expr = print(plot),
             error = function(e){
               PLOT_KAPLAN <- ggsurvplot(fco_mod$KAPLAN, data=fco_mod$IPD, xlim = c(0, tmax),
                 break.x.by = cortes, palette= "#E7B800")
               plot <- PLOT_KAPLAN$plot +
                 geom_line(data = puntos[[i]], aes(x=time, y=est), colour = "#2E9FDF", size =0.5, alpha=0.8) +
                 theme(legend.position = "none")+
                 labs(title = fco_mod$NOMBRE_MODELO[[i]],
                   subtitle = round(fco_mod$modelos[[i]]$AIC,2))
               print(plot)})
  ggsave(filename=paste0("./RESULTADOS/",tipo,"_",fco_mod$FARMACO,"_",fco_mod$NOMBRE_MODELO[[i]],".pdf"),
           plot = last_plot(), dpi=2400)
  }
  output <- puntos
  return(list(puntos=output))}
#### MEDIANAS ####

MEDIANAS <- function(fco_mod)
{
  fco_mod <- fco_mod
  mediana_sLn <- qlnorm(0.5, fco_mod$modelos$Lognormal$res[1,1], fco_mod$modelos$Lognormal$res[2,1])
  mediana_sLlogis <- qllogis(0.5, fco_mod$modelos$Loglogistic$res[1,1], fco_mod$modelos$Loglogistic$res[2,1])
  mediana_sGamma <- qgamma(0.5, fco_mod$modelos$Gamma$res[1,1],  fco_mod$modelos$Gamma$res[2,1])
  mediana_sWei <- qweibull(0.5, fco_mod$modelos$Weibull$res[1,1], fco_mod$modelos$Weibull$res[2,1])
  mediana_sGo <- qgompertz(0.5, fco_mod$modelos$Gompertz$res[1,1], fco_mod$modelos$Gompertz$res[2,1])
  mediana_sExp <- qexp(0.5, fco_mod$modelos$Exponential$res[1,1])
  mediana_sGen <- qgengamma(0.5, fco_mod$modelos$Generalized_gamma$res[1,1],
                            fco_mod$modelos$Generalized_gamma$res[2,1],
                            fco_mod$modelos$Generalized_gamma$res[3,1])
  mediana_sGf <- qgenf(0.5, fco_mod$modelos$Generalized_F$res[1,1],fco_mod$modelos$Generalized_F$res[2,1], fco_mod$modelos$Generalized_F$res[3,1], fco_mod$modelos$Generalized_F$res[4,1])

  # mediana_spline1_hazard <- (qsurvspline(0.5, fco_mod$modelos$spline1_hazard$coefficients,
  #                                        beta=0, X=0, scale="hazard", fco_mod$modelos$spline1_hazard$knots,
  #                                        timescale="log", offset=0, lower.tail=TRUE, log.p=FALSE))
  # 
  # mediana_spline2_hazard <- (qsurvspline(0.5, fco_mod$modelos$spline2_hazard$coefficients,
  #                                        beta=0, X=0, scale="hazard", fco_mod$modelos$spline2_hazard$knots,
  #                                        timescale="log", offset=0, lower.tail=TRUE, log.p=FALSE))
  # 
  # mediana_spline3_hazard <- (qsurvspline(0.5, fco_mod$modelos$spline3_hazard$coefficients,
  #                                        beta=0, X=0, scale="hazard", fco_mod$modelos$spline3_hazard$knots,
  #                                        timescale="log", offset=0, lower.tail=TRUE, log.p=FALSE))
  # 
  # mediana_spline4_hazard <- (qsurvspline(0.5, fco_mod$modelos$spline4_hazard$coefficients,
  #                                        beta=0, X=0, scale="hazard", fco_mod$modelos$spline4_hazard$knots,
  #                                        timescale="log", offset=0, lower.tail=TRUE, log.p=FALSE))
  # 
  # mediana_spline5_hazard <- (qsurvspline(0.5, fco_mod$modelos$spline5_hazard$coefficients,
  #                                        beta=0, X=0, scale="hazard", fco_mod$modelos$spline5_hazard$knots,
  #                                        timescale="log", offset=0, lower.tail=TRUE, log.p=FALSE))
  # 
  # mediana_spline6_hazard <- (qsurvspline(0.5, fco_mod$modelos$spline6_hazard$coefficients,
  #                                        beta=0, X=0, scale="hazard", fco_mod$modelos$spline6_hazard$knots,
  #                                        timescale="log", offset=0, lower.tail=TRUE, log.p=FALSE))
  # 
  # mediana_spline1_odds <- (qsurvspline(0.5, fco_mod$modelos$spline1_odds$coefficients,
  #                                      beta=0, X=0, scale="odds", fco_mod$modelos$spline1_odds$knots,
  #                                      timescale="log", offset=0, lower.tail=TRUE, log.p=FALSE))
  # 
  # mediana_spline2_odds <- (qsurvspline(0.5, fco_mod$modelos$spline2_odds$coefficients,
  #                                      beta=0, X=0, scale="odds", fco_mod$modelos$spline2_odds$knots,
  #                                      timescale="log", offset=0, lower.tail=TRUE, log.p=FALSE))
  # 
  # mediana_spline3_odds <- (qsurvspline(0.5, fco_mod$modelos$spline3_odds$coefficients,
  #                                      beta=0, X=0, scale="odds", fco_mod$modelos$spline3_odds$knots,
  #                                      timescale="log", offset=0, lower.tail=TRUE, log.p=FALSE))
  # 
  # mediana_spline4_odds <- (qsurvspline(0.5, fco_mod$modelos$spline4_odds$coefficients,
  #                                      beta=0, X=0, scale="odds", fco_mod$modelos$spline4_odds$knots,
  #                                      timescale="log", offset=0, lower.tail=TRUE, log.p=FALSE))
  # 
  # mediana_spline5_odds <- (qsurvspline(0.5, fco_mod$modelos$spline5_odds$coefficients,
  #                                      beta=0, X=0, scale="odds", fco_mod$modelos$spline5_odds$knots,
  #                                      timescale="log", offset=0, lower.tail=TRUE, log.p=FALSE))
  # 
  # mediana_spline6_odds <- (qsurvspline(0.5, fco_mod$modelos$spline6_odds$coefficients,
  #                                      beta=0, X=0, scale="odds", fco_mod$modelos$spline6_odds$knots,
  #                                      timescale="log", offset=0, lower.tail=TRUE, log.p=FALSE))
  # 
  # mediana_spline1_normal <- (qsurvspline(0.5, fco_mod$modelos$spline1_normal$coefficients,
  #                                        beta=0, X=0, scale="normal", fco_mod$modelos$spline1_normal$knots,
  #                                        timescale="log", offset=0, lower.tail=TRUE, log.p=FALSE))
  # 
  # mediana_spline2_normal <- (qsurvspline(0.5, fco_mod$modelos$spline2_normal$coefficients,
  #                                        beta=0, X=0, scale="normal", fco_mod$modelos$spline2_normal$knots,
  #                                        timescale="log", offset=0, lower.tail=TRUE, log.p=FALSE))
  # 
  # mediana_spline3_normal <- (qsurvspline(0.5, fco_mod$modelos$spline3_normal$coefficients,
  #                                        beta=0, X=0, scale="normal", fco_mod$modelos$spline3_normal$knots,
  #                                        timescale="log", offset=0, lower.tail=TRUE, log.p=FALSE))
  # 
  # mediana_spline4_normal <- (qsurvspline(0.5, fco_mod$modelos$spline4_normal$coefficients,
  #                                        beta=0, X=0, scale="normal", fco_mod$modelos$spline4_normal$knots,
  #                                        timescale="log", offset=0, lower.tail=TRUE, log.p=FALSE))
  # 
  # mediana_spline5_normal <- (qsurvspline(0.5, fco_mod$modelos$spline5_normal$coefficients,
  #                                        beta=0, X=0, scale="normal", fco_mod$modelos$spline5_normal$knots,
  #                                        timescale="log", offset=0, lower.tail=TRUE, log.p=FALSE))
  # 
  # mediana_spline6_normal <- (qsurvspline(0.5, fco_mod$modelos$spline6_normal$coefficients,
  #                                        beta=0, X=0, scale="normal", fco_mod$modelos$spline6_normal$knots,
  #                                        timescale="log", offset=0, lower.tail=TRUE, log.p=FALSE))

  output <- list("Exponential" = mediana_sExp,
                                 "Gamma" = mediana_sGamma,
                                 "Loglogistic"= mediana_sLlogis,
                                 "Weibull"=mediana_sLlogis,
                                 "Lognormal" = mediana_sLn,
                                 "Gompertz" = mediana_sGo,
                                 "Generalized_gamma" = mediana_sGen,
                                 "Generalized_F" = mediana_sGf)
                                 # "spline1_hazard" = mediana_spline1_hazard,
                                 # "spline2_hazard" = mediana_spline2_hazard,
                                 # "spline3_hazard" = mediana_spline3_hazard,
                                 # "spline4_hazard" = mediana_spline4_hazard,
                                 # "spline5_hazard" = mediana_spline5_hazard,
                                 # "spline6_hazard" = mediana_spline6_hazard,
                                 # "spline1_odds" = mediana_spline1_odds,
                                 # "spline2_odds" = mediana_spline2_odds,
                                 # "spline3_odds" = mediana_spline3_odds,
                                 # "spline4_odds" = mediana_spline4_odds,
                                 # "spline5_odds" = mediana_spline5_odds,
                                 # "spline6_odds" = mediana_spline6_odds,
                                 # "spline1_normal" = mediana_spline1_normal,
                                 # "spline2_normal" = mediana_spline2_normal,
                                 # "spline3_normal" = mediana_spline3_normal,
                                 # "spline4_normal" = mediana_spline4_normal,
                                 # "spline5_normal" = mediana_spline5_normal,
                                 # "spline6_normal" = mediana_spline6_normal)
  return(output)}


#### MARKOV ####

MARKOV <- function(fco_mod,
                   tmax = 100,
                   cortes = 1)
{
  tmax <- tmax
  cortes <- cortes
  wb_fco <- createWorkbook()
  lista <- c()
  for (i in 1:length(fco_mod$modelos)) {
    puntos_markov_fco <- data.frame(summary(fco_mod$modelos[[i]], type="survival", t=seq(0,tmax,by=cortes)))
    TABLA_fco <- cbind(seq(nrow(puntos_markov_fco)),puntos_markov_fco)
    colnames(TABLA_fco) <- c("cicle","time","est","lcl","ucl")
    TABLA.est_fco <- data.frame(c(1,TABLA_fco$est))
    TABLA.est_fco <- data.frame(TABLA.est_fco[-nrow(TABLA.est_fco),])
    TABLA.est_fco <- TABLA_fco$est/TABLA.est_fco
    colnames(TABLA.est_fco) <- "MARKOV"
    TABLA_fco <- cbind(TABLA_fco,TABLA.est_fco)
    TABLA_fco$cicle <- TABLA_fco$cicle-1
    lista <- c(lista, list(TABLA_fco))
    names(lista)[i] <- fco_mod$NOMBRE_MODELO[[i]]
    addWorksheet(wb_fco, fco_mod$NOMBRE_MODELO[[i]])
    writeData(wb_fco, sheet = fco_mod$NOMBRE_MODELO[[i]], x = TABLA_fco)
    saveWorkbook(wb_fco, file=paste0("./RESULTADOS/markov_",fco_mod$FARMACO,".xlsx"), overwrite = TRUE)
  }
  return(lista)}

#### AUC ####

AUC <- function(fco_mod,
                tmin = 0,
                tmax = 100,
                cortes = 1)
{
  tmin <- tmin
  tmax <- tmax
  cortes <- cortes
  puntos <- c()
  for (i in 1:length(fco_mod$modelos)) {
    add <- data.frame(summary(fco_mod$modelos[[i]], type="survival",t=seq(0,tmax,by=cortes)))
    puntos <- c(puntos, list(add))
    names(puntos)[i] <- fco_mod$NOMBRE_MODELO[[i]]
  }
  lista <- c()
  for (i in 1:length(fco_mod$modelos)){
    data <- auc(puntos[[i]]$time,
                puntos[[i]]$est,
                from=tmin,
                to=tmax)
    lista <- c(lista, list(data))
    names(lista)[i] <- fco_mod$NOMBRE_MODELO[[i]]
  }
  return(AUC = lista)
}

