reliability_diagram = function(data_ls,u, stat_type,bins,

                               title_ls,color_ls,title) {

  data_ls_len = length(data_ls)

  ### create bin averages

  for (data_i in 1:data_ls_len){

    colnames(data_ls[[data_i]]) = c("status","survival","pred")

    temp_res = reliability_datapts(data_ls[[data_i]], u = u,bins = bins, stat_type = stat_type)

    temp_hl = hosmer_lemeshow(data_ls[[data_i]], u = u,bins = bins, stat_type = stat_type)

    assign(paste("recal_bins", data_i, sep=""),temp_res)

    assign(paste("HL_pvalue", data_i, sep=""),temp_hl[1])

    assign(paste("HL_statistic", data_i, sep=""),temp_hl[2])

  }

  hl_dt = data.frame(x=numeric(0),y=numeric(0),variable= character(0), label = character(0),HLpval = numeric(0),HLstatistic =numeric(0) )

  for (data_i in 1:data_ls_len){

    temp_res = get(paste('recal_bins', data_i, sep=''))

    temp_res$variable <- paste("Vol.x", data_i, sep='')

    colnames(temp_res)[colnames(temp_res)=='V1'] = 'value'

    assign(paste('melt', data_i, sep=''), temp_res)

    hl_p = get(paste('HL_pvalue', data_i, sep=''))

    hl_s = get(paste('HL_statistic', data_i, sep=''))

    hl_dt[data_i,] = cbind(0.05,1-0.1*data_i,paste("Vol.x", data_i, sep=''),title_ls[data_i],hl_p[[1]],hl_s[[1]])

  }

  hl_dt$x = as.numeric(hl_dt$x)

  hl_dt$y = as.numeric(hl_dt$y)

  hl_dt$HLpval = as.numeric(hl_dt$HLpval)

  hl_dt$HLstatistic = as.numeric(hl_dt$HLstatistic)

  data = melt1

  if (data_ls_len > 1){

    for (data_i in 2:data_ls_len){

      data = rbind(data, get(paste('melt', data_i, sep='')))

    }

  }

  line_plot = ggplot(data, aes(x=pred,  y=obs,color=variable,group= variable))   +  

    geom_point(aes(shape = variable),size = 3)+

    geom_line(linetype = "solid")+

    geom_text(data = hl_dt,aes(x=x,y=y,color = variable,hjust=0,

                               label = paste0(label," H-L test p=",round(HLpval,2),", C-statistic=",round(HLstatistic,2))))+

    scale_color_manual(labels = title_ls,

                       values = color_ls) +

    scale_shape_discrete(labels = title_ls)+

    scale_linetype_discrete(labels = title_ls)+

    guides(color=guide_legend(" "),shape= guide_legend(" "),linetype=guide_legend(" ")) + 

    xlab(paste("Predicted Survival Probability at",title)) + ylab(paste("Observed Proportion at",title)) +

    geom_abline(intercept = 0, slope = 1, color="black",

                linetype="dashed", size=1) +

    lims(x=c(0,1),y=c(0,1))+

    theme_bw()+

    theme(

      legend.position="bottom",

      axis.title = element_text(size=12),

      axis.text = element_text(size=12),

      legend.text = element_text(size=15),

      legend.title = element_text(size=15),

      text = element_text(size=12)

    )

  line_plot

  if (!file.exists("results/LOO_calibration")) dir.create('results/LOO_calibration')

  ggsave(paste0("results/LOO_calibration/",title,".pdf"),device = 'pdf',width = 6,height = 6)

}

reliability_datapts <- function(ndata, u,bins=10, stat_type ='C') {

  min.pred <- min(ndata$pred)

  max.pred <- max(ndata$pred)

  min.max.diff <- max.pred - min.pred

  if (stat_type == 'H'){

    ndata = ndata[order(ndata$pred),]

    res = data.frame(obs= numeric(0), pred = numeric(0),obs_lower=numeric(0),obs_upper=numeric(0))

    split_mtx = split(ndata, cut(ndata$pred, seq(0,1,1/bins), include.lowest=TRUE))

    for (i in 1:length(split_mtx)){

      sfit = summary(survfit(Surv(survival, status) ~ 1,data=split_mtx[[i]]),times=u,extend=T)

      obs = sfit[['surv']]

      obs_upper =  sfit[['upper']]

      obs_lower = sfit[['lower']]

      pred = mean(split_mtx[[i]]$pred)

      if (sum(is.na(col_mean)) > 0) {

        next

      }

      res[i,] = c(obs,pred,obs_lower,obs_upper)

    }

  }else{

    ## C statistics, same number of instances in each bin

    mtx = ndata[order(ndata$pred),]

    n <- length(ndata$pred)/bins

    nr <- nrow(mtx)

    split_mtx = split(mtx, rep(1:ceiling(nr/n), each=n, length.out=nr))

    res = data.frame(obs= numeric(0), pred = numeric(0),obs_lower=numeric(0),obs_upper=numeric(0))

    for (i in 1:length(split_mtx)){

      sfit = summary(survfit(Surv(survival, status) ~ 1,data=split_mtx[[i]]),times=u,extend=T)

      obs = sfit[['surv']]

      obs_upper =  sfit[['upper']]

      obs_lower = sfit[['lower']]

      pred = mean(split_mtx[[i]]$pred)

      res[i,] = c(obs,pred,obs_lower,obs_upper)

    }

  }

  res

}

hosmer_lemeshow <- function(ndata, u,bins=10, stat_type ='C') {

  min.pred <- min(ndata$pred)

  max.pred <- max(ndata$pred)

  min.max.diff <- max.pred - min.pred

  if (stat_type == 'H'){

    ndata = ndata[order(ndata$pred),]

    res = data.frame(obs= numeric(0), pred = numeric(0),obs_lower=numeric(0),obs_upper=numeric(0))

    split_mtx = split(ndata, cut(ndata$pred, seq(0,1,1/bins), include.lowest=TRUE))

  }else{

    mtx = ndata[order(ndata$pred),]

    n <- length(ndata$pred)/bins

    nr <- nrow(mtx)

    split_mtx = split(mtx, rep(1:ceiling(nr/n), each=n, length.out=nr))

  }

  H_stat = 0

  for (i in 1:length(split_mtx)){

    sfit = summary(survfit(Surv(survival, status) ~ 1,data=split_mtx[[i]]),times=u,extend=T)

    obs = sfit[['surv']]

    exp = mean(split_mtx[[i]]$pred)

    obs_not = 1-obs

    exp_not = 1-exp

    if (exp == 0 || exp_not == 0){

      next

    }

    bin_sum = ((obs - exp)**2)/exp + ((obs_not - exp_not)**2)/exp_not

    H_stat = H_stat + bin_sum

  }

  PVAL = 1 - pchisq(H_stat, bins - 2)

  cat('PVALUE', PVAL, '\n')

  cat('stat', H_stat, '\n')

  return(c(PVAL,H_stat))

}

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

# data_ls = list(subset(dt_risk_LOO,!is.na(prob1)) %>% select(DFS_status,DFS,prob1) ,

#      subset(dt_risk_LOO,!is.na(cox1_prob1))%>% select(DFS_status,DFS,cox1_prob1))

# title_ls= c('Joint model', 'Cox model')

# color_ls = c('blue', 'red')

# limits=c(0,1)

# u=15

# bins = 5

# title = "15 Months"