##### library #####
if (!require('DT')) install.packages('DT'); library('DT') #### functions: datatable
if (!require('flexsurv')) install.packages('flexsurv'); library('flexsurv') #### functions: ggsurvplot
if (!require('openxlsx')) install.packages('openxlsx'); library('openxlsx') #### functions: write.xlsx
if (!require('readr')) install.packages('readr'); library('readr') #### functions: read_delim 
if (!require('readxl')) install.packages('readxl'); library('readxl')
if (!require('survminer')) install.packages('survminer'); library('survminer') #### functions: ggsurvplot
if (!require('MESS')) install.packages('MESS'); library('MESS') #### functions: auc
if (!require('heemod')) install.packages('heemod'); library('heemod') #### functions: auc

#### GITHUB #####
if (!require('remotes')) install.packages('remotes'); library('remotes')
if (!require('IPDfromKM')) remotes::install_github("NaLiuStat/IPDfromKM"); library('IPDfromKM')

# setwd - comando funciona desde rstudio
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
#### 
source("FUNCION_GUYOT_1.R")

##### input file ####
input <- c()
  input$trial <- "PEMBRO"
input$fco1 <- "PEMBRO+QT"
input$fco2 <- "QT" 
input$cycle_length <- 1  ##### length of the cycle
input$tmax <- 100

###### PFS fco1 GUYOT #####
brazo <- 1
input$PFS_triski_fco1 <- seq(from=0, to=48, by=3)
input$PFS_Nriski_fco1 <- c(784,781,769,751,728,718,702,692,681,671,652,551,433,303,165,28,0)
input$PFS_total.events_fco1 <- NULL   # IF YOU DO NOT KNOW ENTER NULL
input$PFS_file_points_fco1 <- "PFS_PEMBRO_INTV.csv"
PFS_fco1 <- GUYOT(
  ensayo = input$trial,
  farmaco = input$fco1,
  brazo = brazo,
  triski_fco = input$PFS_triski_fco1,
  Nriski_fco = input$PFS_Nriski_fco1,
  totalpts =  input$PFS_Nriski_fco1[1],
  total.events_fco = input$PFS_total.events_fco1,
  datos_digitalizados = input$PFS_file_points_fco1,
  cortes_FP = input$cycle_length)

###### PFS fco2 GUYOT #####
brazo <- 2
input$PFS_triski_fco2 <- seq(from=0, to=48, by=3)
input$PFS_Nriski_fco2 <- c(390,386,382,358,342,328,319,310,304,297,250,195,140,83,17,0)
input$PFS_total.events_fco2 <- NULL  # IF YOU DO NOT KNOW ENTER NULL
input$PFS_file_points_fco2 <- "PFS_PEMBRO_CONT.csv"
PFS_fco2 <- GUYOT(
  ensayo = input$trial,
  farmaco = input$fco2,
  brazo = brazo,
  triski_fco = input$PFS_triski_fco2,
  Nriski_fco = input$PFS_Nriski_fco2,
  totalpts =  input$PFS_Nriski_fco2[1],
  total.events_fco = input$PFS_total.events_fco2,
  datos_digitalizados = input$PFS_file_points_fco2,
  cortes_FP = input$cycle_length)




#### PFS fco1 modeling ####
PFS_fco1_mod <- MODELIZACION(PFS_fco1)
PFS_fco1_mod$modelos
MEDIANAS(PFS_fco1_mod)

#### PFS fco2 modeling ####
PFS_fco2_mod <- MODELIZACION(PFS_fco2)
PFS_fco2_mod$modelos
MEDIANAS(PFS_fco2_mod)


save(PFS_fco1_mod,PFS_fco2_mod,
     file="modelos_PEMBRO.RData")

library(survival)
PFS_fco2_mod$IPD$treat<-0 #recode control group
cox_data<-data.frame(rbind(PFS_fco1_mod$IPD,PFS_fco2_mod$IPD))
cox_model<-coxph(Surv(time,status)~treat,cox_data, x=TRUE)

#HR and 95% confidence limit
exp(cox_model$coefficients)
exp(cox_model$coefficients-1.96*sqrt(cox_model$var))
exp(cox_model$coefficients+1.96*sqrt(cox_model$var))

#test of proportional hazards assumption
temp<-cox.zph(cox_model)
temp
round(temp$table["treat","p"],2) # test of proportional hazards assumption 
