## cross validation

runCV = function(i,Tstart,predict_times,data,id){

  library(dplyr)

  trainingData <- data[which(data$PatientID %in% id[i,1]),]

  testingData <- data[which(data$PatientID %in% id[-i,1]),]

  trainingData.id <- trainingData[!duplicated(trainingData$PatientID), ]

  testingData.id <- testingData[!duplicated(testingData$PatientID), ]

  CV_results <- function(trainingData, testingData.id, predict_times){

    test <- try({

      ##################

      # Joint Modeling #

      ##################

      library("JMbayes")

      library("splines")

      library("survminer")

      library(survival)

      library(dplyr)

      lmeFit1 <- lme( logAF ~ ns(TestDate,2),data=trainingData,random = ~  ns(TestDate,2) | PatientID,

                      control = lmeControl(opt = "optim",msMaxIter =10000))

      coxFit1 <- coxph(Surv(DFS,DFS_status)~T_stage+TP53,data=trainingData.id,x = TRUE)

      # origin

      jointFit1 = try(jointModelBayes(lmeFit1, coxFit1, timeVar = "TestDate",n.iter = 20000L),TRUE)

      # value + slope

      dForm <- list(fixed = ~ 0 + dns(TestDate, 2), random = ~ 0 + dns(TestDate, 2),

                    indFixed = 2:3, indRandom = 2:3)

      jointFit2 <- try(jointModelBayes(lmeFit1, coxFit1, timeVar = "TestDate",

                                       param = "td-both", extraForm = dForm),T)

      # value + cumulative

      iForm <- list(fixed = ~ 0 + TestDate + ins(TestDate, 2), random = ~ 0 + TestDate + ins(TestDate, 2),

                    indFixed = 1:3, indRandom = 1:3)

      jointFit3 <-try(jointModelBayes(lmeFit1, coxFit1, timeVar = "TestDate",

                                      param = "td-extra", extraForm = iForm) ,T)

      ######################################

      JM_Models <- list("JM_value" = jointFit1, 

                        "JM_value-slope" = jointFit2,

                        "JM_cumulative" = jointFit3

      )

      combos <- expand.grid("model_name" = names(JM_Models), "time" = predict_times)

      ## test 

      # AUC

      # time = 12

      # jm = jointFit3

      auc_fun <- function (time) {

        auc_objs <- mapply(function(jm,newdata=testingData,Tstart=Tstart,Thoriz =time){

          res = try(aucJM(jm,newdata=testingData,Tstart=Tstart,Thoriz =time,idVar ="PatientID" ,simulate=T,M=100),T)

          if (class(res)=="try-error") res = list('auc'=NA)

          return(res)

        }, JM_Models, 

        MoreArgs = list(newdata=testingData,Tstart=Tstart,Thoriz =time), 

        SIMPLIFY = FALSE

        )

        sapply(auc_objs, "[[", "auc")

      }

      JM_AUCs <- sapply(predict_times, auc_fun)

      # PE

      pe_fun <- function (time) {

        pe_objs <- mapply(function(jm,newdata=testingData,Tstart=Tstart,Thoriz =time){

          res = try(prederrJM(jm,newdata=testingData,Tstart=Tstart,Thoriz =time,idVar ="PatientID", lossFun = "square" ,simulate=T,M=100),T)

          if (class(res)=="try-error") res = list('prederr'=NA)

          return(res)

        }, JM_Models, MoreArgs = list(newdata=testingData,Tstart=Tstart,Thoriz =time), SIMPLIFY = FALSE)

        sapply(pe_objs, "[[", "prederr")

      }

      JM_PEs <- sapply(predict_times, pe_fun)

      ## training 

      auc_fun_tr <- function (time) {

        auc_objs <- mapply(function(jm,newdata=trainingData,Tstart=Tstart,Thoriz =time){

          res = try(aucJM(jm,newdata=trainingData,Tstart=Tstart,Thoriz =time,idVar ="PatientID" ,simulate=T,M=100),T)

          if (class(res)=="try-error") res = list('auc'=NA)

          return(res)

        }, JM_Models, 

        MoreArgs = list(newdata=trainingData,Tstart=Tstart,Thoriz =time), SIMPLIFY = FALSE)

        sapply(auc_objs, "[[", "auc")

      }

      JM_AUCs_tr <- sapply(predict_times, auc_fun_tr)

      # PE

      pe_fun_tr <- function (time) {

        pe_objs <- mapply(function(jm,newdata=trainingData,Tstart=Tstart,Thoriz =time){

          res = try(prederrJM(jm,newdata=trainingData,Tstart=Tstart,Thoriz =time,idVar ="PatientID", lossFun = "square" ,simulate=T,M=100),T)

          if (class(res)=="try-error") res = list('prederr'=NA)

          return(res)

        }, JM_Models, MoreArgs = list(newdata=trainingData,Tstart=Tstart,Thoriz =time), SIMPLIFY = FALSE)

        sapply(pe_objs, "[[", "prederr")

      }

      JM_PEs_tr <- sapply(predict_times, pe_fun_tr)

      ###############

      # cox #

      ###############

      dataLM <- JMbayes:::dataLM

      LM_models_fun <- function () {

        # Landmark data sets

        D1 = dataLM(trainingData, Tstart, idVar ="PatientID" ,respVar = "Test_status", timeVar = "TestDate", 

                    evTimeVar = "DFS")

        CoxLM1 <- coxph(Surv(DFS, DFS_status) ~ Test_status+T_stage+TP53, 

                        data = D1)

        CoxLM2 <- coxph(Surv(DFS, DFS_status) ~ P2+T_stage+TP53, 

                        data = D1)

        list("LM_bi" = CoxLM1,

             LM_P2 = CoxLM2 

        )

      }  

      ######################################################################################

      # Calculate Performance Measures

      ## test

      # AUC

      LM_models <- LM_models_fun()

      auc_fun <- function (time) {

        auc_objs <- mapply(function(jm,newdata0,Tstart0,Thoriz0){

          if ("P2" %in% names(jm$assign)){

              newdata0 = subset(newdata0,!is.na(P2))

          }

          res = try(aucJM(jm,newdata=newdata0,Tstart=Tstart0,Thoriz =Thoriz0,

                          idVar ="PatientID" ,respVar = "Test_status", 

                          timeVar = "TestDate",evTimeVar = "DFS"),TRUE)

          if (class(res)=="try-error"){res=list('auc'=NA)}

          return(res)

        }, 

        LM_models, 

        MoreArgs = list(newdata0=testingData,Tstart0=Tstart,Thoriz0 =time), SIMPLIFY = FALSE)

        sapply(auc_objs, "[[", "auc")

      }

      LM_AUCs <- sapply(predict_times, auc_fun)

      # PE

      pe_fun <- function (time) {

        pe_objs <- mapply(function(jm,newdata0,Tstart0,Thoriz0){

          if ("P2" %in% names(jm$assign)){

            newdata0 = subset(newdata0,!is.na(P2))

          }

          res = try(prederrJM(jm,newdata=newdata0,Tstart=Tstart0,Thoriz =Thoriz0,

                              idVar ="PatientID" ,respVar = "Test_status", 

                              timeVar = "TestDate",evTimeVar = "DFS", 

                              lossFun = "square"),TRUE)

          if (class(res)=="try-error"){res=list('prederr'=NA)}

          return(res)

        },

        LM_models, 

        MoreArgs = list(newdata0=testingData,Tstart0=Tstart,Thoriz0 =time), SIMPLIFY = FALSE)

        sapply(pe_objs, "[[", "prederr")

      }

      LM_PEs <- sapply(predict_times, pe_fun)

      ## training

      # AUC

      auc_fun_tr <- function (time) {

        auc_objs <- mapply(function(jm,newdata0,Tstart0,Thoriz0){

          if ("P2" %in% names(jm$assign)){

            newdata0 = subset(newdata0,!is.na(P2))

          }

          res = try(aucJM(jm,newdata=newdata0,Tstart=Tstart0,Thoriz =Thoriz0,

                          idVar ="PatientID" ,respVar = "Test_status", 

                          timeVar = "TestDate",evTimeVar = "DFS"),TRUE)

          if (class(res)=="try-error"){res=list('auc'=NA)}

          return(res)

        }, 

        LM_models, 

        MoreArgs = list(newdata0=trainingData,Tstart0=Tstart,Thoriz0 =time), SIMPLIFY = FALSE)

        sapply(auc_objs, "[[", "auc")

      }

      LM_AUCs_tr <- sapply(predict_times, auc_fun_tr)

      # PE

      pe_fun_tr <- function (time) {

        pe_objs <- mapply(function(jm,newdata0,Tstart0,Thoriz0){

          if ("P2" %in% names(jm$assign)){

            newdata0 = subset(newdata0,!is.na(P2))

          }

          res = try(prederrJM(jm,newdata=newdata0,Tstart=Tstart0,Thoriz =Thoriz0,

                              idVar ="PatientID" ,respVar = "Test_status", 

                              timeVar = "TestDate",evTimeVar = "DFS", 

                              lossFun = "square"),TRUE)

          if (class(res)=="try-error"){res=list('prederr'=NA)}

          return(res)

        },

        LM_models, 

        MoreArgs = list(newdata0=trainingData,Tstart0=Tstart,Thoriz0 =time), SIMPLIFY = FALSE)

        sapply(pe_objs, "[[", "prederr")

      }

      LM_PEs_tr <- sapply(predict_times, pe_fun_tr)

      ######################################################################################

      # Collect results

      colnames(JM_AUCs) = colnames(JM_PEs) = colnames(LM_AUCs) = colnames(LM_PEs)  = predict_times

      colnames(JM_AUCs_tr) = colnames(JM_PEs_tr) =  colnames(LM_AUCs_tr) = colnames(LM_PEs_tr) = predict_times

      JM = rbind(as.data.frame(JM_AUCs)%>%tibble::rownames_to_column(var="model2")%>%mutate(variable="AUC",type="testing"),

                 as.data.frame(JM_PEs)%>%tibble::rownames_to_column(var="model2")%>%mutate(variable="PE",type="testing"),

                 as.data.frame(JM_AUCs_tr)%>%tibble::rownames_to_column(var="model2")%>%mutate(variable="AUC",type="training"),

                 as.data.frame(JM_PEs_tr)%>%tibble::rownames_to_column(var="model2")%>%mutate(variable="PE",type="training")

      )

      cox = rbind(as.data.frame(LM_AUCs)%>%tibble::rownames_to_column(var="model2")%>%mutate(variable="AUC",type="testing"),

                  as.data.frame(LM_PEs)%>%tibble::rownames_to_column(var="model2")%>%mutate(variable="PE",type="testing"),

                  as.data.frame(LM_AUCs_tr)%>%tibble::rownames_to_column(var="model2")%>%mutate(variable="AUC",type="training"),

                  as.data.frame(LM_PEs_tr)%>%tibble::rownames_to_column(var="model2")%>%mutate(variable="PE",type="training")

      )

      results_accuracy = rbind(JM%>%mutate(model1="JM"),cox%>%mutate(model1="cox"))

      results_accuracy

    },TRUE)

    if (!inherits(test, "try-error")) {

      return(test)

    } else {

      JM_AUCs = matrix(nrow = 3,ncol=length(predict_times),

                       dimnames = list(c("JM_value","JM_value-slope","JM_cumulative"),predict_times))

      JM_PEs = matrix(nrow = 3,ncol=length(predict_times),

                      dimnames = list(c("JM_value","JM_value-slope","JM_cumulative"),predict_times))

      LM_AUCs = matrix(nrow = 2,ncol=length(predict_times),

                       dimnames = list(c("LM_bi","LM_P2"),predict_times))

      LM_PEs = matrix(nrow = 2,ncol=length(predict_times),

                      dimnames = list(c("LM_bi","LM_P2"),predict_times))

      JM_AUCs_tr = matrix(nrow = 3,ncol=length(predict_times),

                          dimnames = list(c("JM_value","JM_value-slope","JM_cumulative"),predict_times))

      JM_PEs_tr = matrix(nrow = 3,ncol=length(predict_times),

                         dimnames = list(c("JM_value","JM_value-slope","JM_cumulative"),predict_times))

      LM_AUCs_tr = matrix(nrow = 2,ncol=length(predict_times),

                          dimnames = list(c("LM_bi","LM_P2"),predict_times))

      LM_PEs_tr = matrix(nrow = 2,ncol=length(predict_times),

                         dimnames = list(c("LM_bi","LM_P2"),predict_times))

      JM = rbind(as.data.frame(JM_AUCs)%>%tibble::rownames_to_column(var="model2")%>%mutate(variable="AUC",type="testing"),

                 as.data.frame(JM_PEs)%>%tibble::rownames_to_column(var="model2")%>%mutate(variable="PE",type="testing"),

                 as.data.frame(JM_AUCs_tr)%>%tibble::rownames_to_column(var="model2")%>%mutate(variable="AUC",type="training"),

                 as.data.frame(JM_PEs_tr)%>%tibble::rownames_to_column(var="model2")%>%mutate(variable="PE",type="training")

      )

      cox = rbind(as.data.frame(LM_AUCs)%>%tibble::rownames_to_column(var="model2")%>%mutate(variable="AUC",type="testing"),

                  as.data.frame(LM_PEs)%>%tibble::rownames_to_column(var="model2")%>%mutate(variable="PE",type="testing"),

                  as.data.frame(LM_AUCs_tr)%>%tibble::rownames_to_column(var="model2")%>%mutate(variable="AUC",type="training"),

                  as.data.frame(LM_PEs_tr)%>%tibble::rownames_to_column(var="model2")%>%mutate(variable="PE",type="training")

      )

      results_accuracy = rbind(JM%>%mutate(model1="JM"),cox%>%mutate(model1="cox"))

      results_accuracy

    }

  }

  CV_results(trainingData,testingData,predict_times)

}

##### test #########

# i= splits$Fold1.Rep02

# data = mydata

# id = mydata.id

# runCV(i,Tstart,predict_times,mydata,mydata.id)