R/raptools.R

In raptools: Risk Assessment Plot and Reclassification Metrics

# Global variables to avoid NSE issues with ggplot2 and dplyr
globalVariables(c("Baseline", "New", "Event", "Model", "Probabilities", "ID", 
                  "n_event", "n_nonevent", "n_event_cumulative", "n_nonevent_cumulative",
                  "Sensitivity", "1-Specificity", "IDs", "ID_list", "Metric", "Metric name",
                  "Prediction", "n_ev", "TP", "n_nonev", "TN", "FN", "FP", "spec", "treated",
                  "untreated", "sens", "ppv", "1-spec", "x", "interval", "rn", "mid_rn",
                  "prediction", "calib_baseline", "calib_new", "pred_interval", "prop_event",
                  "mn_pred", "lci_prop_event", "uci_prop_event", "lci_pred", "uci_pred",
                  "Chi-Square", "d.f.", "Xsqmindf", "variables", "Xsqmindf.percent",
                  "prevalence", "Extreme models", "extremes", "V1", "V2", "V3", "metric",
                  "statistics", "event", "up", "down", "n_up", "n_down", "nri",
                  "risk.class.x2", "risk.class.x1", "difference", "IDI", "mse_x1", "mse_x2",
                  "brier_baseline", "brier_new", "lci", "uci"))

#' The Risk Assessment Plot
#'
#' The function ggrap() plots the Sensitivity and 1-Specificity curves against the calculated risk for the baseline (reference) and newmodels, thus graphically displaying the IDIs for those with and without the events.  These plots can aid interpretation of the NRI and IDI metrics.
#' 
#' @param x1 Either a logistic regression fitted using glm (base package) or lrm (rms package) or alculated probabilities (eg through a logistic regression model) of the baseline model.  Must be between 0 & 1
#' @param x2 Either a logistic regression fitted using glm (base package) or lrm (rms package) or calculated probabilities (eg through a logistic regression model) of the new (alternative) model.   Must be between 0 & 1
#' @param y Binary of outcome of interest. Must be 0 or 1 (if fitted models are provided this is extracted from the fit which for an rms fit must have x = TRUE, y = TRUE). 
#' @return a ggplot 
#' @references The Risk Assessment Plot in this form was described by Pickering, J. W., & Endre, Z. H. (2012). New Metrics for Assessing Diagnostic Potential of Candidate Biomarkers. Clinical Journal of the American Society of Nephrology, 7, 1355–1364. doi:10.2215/CJN.09590911
#' @import ggplot2
#' @import tidyr 
#' @import dplyr
#' @importFrom stats anova glm quantile
#' @export
ggrap <- function(x1, x2=NULL, y=NULL) {
  
  if (class(x1)[1] == "glm") {
    y = x1$y
    x1 = stats::predict(x1, type = "response")
    data_type = "glm"
  }
  if (class(x2)[1] == "glm") {
    x2 = stats::predict(x2, type = "response")
  }
  if (class(x1)[1] == "lrm") {
    if (length(x1$y) == 0 ) {
      stop("Fitted models with the rms package must be made using y = TRUE")
    }
    y = as.numeric(as.character(x1$y))
    x1 = stats::predict(x1, type = "fitted")
    data_type = "lrm"
  }
  if (class(x2)[1] == "lrm") {
    if (length(x2$y) == 0 ) {
      stop("Fitted models with the rms package must be made using y = TRUE")
    }
    x2 = stats::predict(x2, type = "fitted")
  }
  
  if (class(x1)[1] != "glm"  & class(x1)[1] != "lrm" ) {
    data_type = "User supplied"
  }
  if (!is.null(x2) & length(x1) != length(x2))
    stop("Reference (baseline) and New (Alt) model vectors must be the same length")
  if (is.null(y))
    stop("Oops - there must be event data (y)")
  
  
  
  if(!is.null(x2)){
    df <- data.frame(Baseline = x1, New = x2, Event = y)
    df <- df |> 
      filter(!is.na(Baseline)) |> 
      filter(!is.na(New)) |> 
      filter(!is.na(Event))
    df$ID = seq(1,nrow(df),1)
    
    df_long <- df |> 
      tidyr::pivot_longer(cols = c(Baseline,New), values_to = "Probabilities", names_to = "Model") |> 
      group_by(Model)  |> 
      arrange(Model, Probabilities) |> 
      group_by(Model,Probabilities) |> 
      reframe(n_event = sum(Event == 1),
              n_nonevent = sum(Event == 0),
              IDs = as.character(ID)) |> 
      ungroup() |> 
      group_by(Model) |> 
      mutate(n_event_cumulative = cumsum(n_event)) |> 
      mutate(n_nonevent_cumulative = cumsum(n_nonevent)) |> 
      mutate(Sensitivity = (max(n_event_cumulative) - n_event_cumulative)/max(n_event_cumulative)) |> 
      mutate(`1-Specificity` = (max(n_nonevent_cumulative) - n_nonevent_cumulative)/max(n_nonevent_cumulative)) |> 
      ungroup()
  }
  
  if(is.null(x2)){
    df <- data.frame(Baseline = x1,  Event = y)
    df <- df |> 
      filter(!is.na(Baseline)) |> 
      filter(!is.na(Event))
    df$ID = seq(1,nrow(df),1)
    
    df_long <- df |> 
      tidyr::pivot_longer(cols = c(Baseline), values_to = "Probabilities", names_to = "Model") |> 
      group_by(Model)  |> 
      arrange(Model, Probabilities) |> 
      group_by(Model,Probabilities) |> 
      reframe(n_event = sum(Event == 1),
              n_nonevent = sum(Event == 0),
              IDs = as.character(ID)) |> 
      ungroup() |> 
      group_by(Model) |> 
      mutate(n_event_cumulative = cumsum(n_event)) |> 
      mutate(n_nonevent_cumulative = cumsum(n_nonevent)) |> 
      mutate(Sensitivity = (max(n_event_cumulative) - n_event_cumulative)/max(n_event_cumulative)) |> 
      mutate(`1-Specificity` = (max(n_nonevent_cumulative) - n_nonevent_cumulative)/max(n_nonevent_cumulative)) |> 
      ungroup()
  }
  
  
  
  
  df_long_av <- df_long |> 
    group_by(Model, Probabilities) |> 
    summarise(Sensitivity = mean(Sensitivity),
              `1-Specificity` = mean(`1-Specificity`),
              ID_list = list(IDs)
    ) |> 
    ungroup()
  
  df_g <- df_long_av |> 
    select(Probabilities, Sensitivity, `1-Specificity`, Model, ID_list) |>
    tidyr::pivot_longer(cols = c(Sensitivity, `1-Specificity`), values_to = "Metric", names_to = "Metric name")
  
  
  if(!is.null(x2)){
    g <- ggplot(df_g, aes(x = Probabilities, y = Metric, colour = `Metric name`, linetype = Model )) + 
      scale_x_continuous(breaks = seq(0,1,0.1), expand = c(0.005,0.005)) + 
      scale_y_continuous(breaks = seq(0,1,0.1), expand = c(0.005,0.005)) +
      geom_line() +
      scale_linetype_manual(values = c("dotted", "solid")) 
  }
  if(is.null(x2)){
    g <- ggplot(df_g, aes(x = Probabilities, y = Metric, colour = `Metric name` )) + 
      scale_x_continuous(breaks = seq(0,1,0.1), expand = c(0.005,0.005)) + 
      scale_y_continuous(breaks = seq(0,1,0.1), expand = c(0.005,0.005)) +
      geom_line() 
  }
  
  g <- g +
    coord_cartesian(xlim = c(0,1), ylim = c(0,1))    
  
  return(g)
}

#' The Decision curve
#' 
#' ggdecision plots decision curves to assess the net benefit at different thresholds
#' 
#' @param x1 Either a logistic regression fitted using glm (base package) or lrm (rms package) or calculated probabilities (eg through a logistic regression model) of the baseline model.  Must be between 0 & 1
#' @param x2 Either a logistic regression fitted using glm (base package) or lrm (rms package) or calculated probabilities (eg through a logistic regression model) of the new (alternative) model.   Must be between 0 & 1
#' @param y Binary of outcome of interest. Must be 0 or 1 (if fitted models are provided this is extracted from the fit which for an rms fit must have x = TRUE, y = TRUE). 
#' @param show_smooth Logical, whether to display smoothed curves (default = TRUE)
#' @param smooth_method Smoothing method for geom_smooth. Options: "loess", "lm", "glm", "gam". Default is "loess"
#' @param smooth_span Span parameter for loess smoothing, controls the degree of smoothing (default = 0.75). Lower values = less smooth
#' @param smooth_se Logical, whether to display confidence interval around smooth (default = FALSE)
#' @return a ggplot 
#' @references Vickers AJ, van Calster B, Steyerberg EW. A simple, step-by-step guide to interpreting decision curve analysis. Diagn Progn Res 2019;3(1):18. 2. Zhang Z, Rousson V, Lee W-C, et al. Decision curve analysis: a technical note. Ann Transl Med 2018;6(15):308-308. 
#' @import ggplot2
#' @import tidyr
#' @import dplyr
#' @export
#'
ggdecision <- function(x1, x2=NULL, y=NULL, show_smooth = TRUE, smooth_method = "loess", smooth_span = 0.75, smooth_se = FALSE) {
  
  if (class(x1)[1] == "glm") {
    y = x1$y
    x1 = stats::predict(x1, type = "response")
    data_type = "glm"
  }
  if (class(x2)[1] == "glm") {
    x2 = stats::predict(x2, type = "response")
  }
  if (class(x1)[1] == "lrm") {
    if (length(x1$y) == 0 ) {
      stop("Fitted models with the rms package must be made using y = TRUE")
    }
    y = as.numeric(as.character(x1$y))
    x1 = stats::predict(x1, type = "fitted")
    data_type = "lrm"
  }
  if (class(x2)[1] == "lrm") {
    if (length(x2$y) == 0 ) {
      stop("Fitted models with the rms package must be made using y = TRUE")
    }
    x2 = stats::predict(x2, type = "fitted")
  }
  
  if (class(x1)[1] != "glm"  & class(x1)[1] != "lrm" ) {
    data_type = "User supplied"
  }
  if (!is.null(x2) & length(x1) != length(x2))
    stop("Reference (baseline) and New (Alt) model vectors must be the same length")
  if (is.null(y))
    stop("Oops - there must be event data (y)")
  
  
  
  # Two lines
  if(!is.null(x2)){
    df <- data.frame(Baseline = x1, New = x2, Event = y)
    df <- df |> 
      filter(!is.na(Baseline)) |> 
      filter(!is.na(New)) |> 
      filter(!is.na(Event))
    df$ID = seq(1,nrow(df),1)
    
    n <- nrow(df)
    n_event = sum(df$Event)
    n_nonevent = n - n_event
    incidence <- 100 * n_event/n
    
    
    # Net Benefit for Decision Curves
    benefit.1 <- df |> 
      mutate(Prediction = Baseline) |> 
      group_by(Prediction) |> 
      summarise(n_ev = sum(Event == 1),
                n_nonev = sum(Event == 0)) |> 
      ungroup() |> 
      mutate(TP = n_event - cumsum(n_ev),
             FN = n_event - TP, 
             TN = cumsum(n_nonev),
             FP = n_nonevent - TN,
             prevalence = (TP + FN)/n,
             sens = TP/(TP + FN),
             spec = TN/(TN + FP),
             `1-spec` = 1 - spec,
             npv = TN/(TN + FN),
             ppv = TP/(TP + FP),
             treated = TP/n - (FP/n) * Prediction/(1 - Prediction),
             untreated = TN/n - (FN/n) * (1 - Prediction)/Prediction,
             overall = treated + untreated,
             all = prevalence - (1 - prevalence)  * Prediction/(1 - Prediction),
             Model = "Baseline"
      )
    
    benefit.2 <- df |> 
      mutate(Prediction = New) |> 
      group_by(Prediction) |> 
      summarise(n_ev = sum(Event == 1),
                n_nonev = sum(Event == 0)) |> 
      ungroup() |> 
      mutate(TP = n_event - cumsum(n_ev),
             FN = n_event - TP, 
             TN = cumsum(n_nonev),
             FP = n_nonevent - TN,
             prevalence = (TP + FN)/n,
             sens = TP/(TP + FN),
             spec = TN/(TN + FP),
             `1-spec` = 1 - spec,
             npv = TN/(TN + FN),
             ppv = TP/(TP + FP),
             treated = TP/n - (FP/n) * Prediction/(1 - Prediction),
             untreated = TN/n - (FN/n) * (1 - Prediction)/Prediction,
             overall = treated + untreated,
             all = prevalence - (1 - prevalence)  * Prediction/(1 - Prediction),
             Model = "New"
      )
    
    benefit = bind_rows(benefit.1, benefit.2)
    
    benefit_all_none <- benefit |> 
      select(Prediction, all) |> 
      mutate(none = 0)  |> 
      tidyr::pivot_longer(cols = c("all", "none"), names_to = "Extreme models", values_to = "extremes") |> 
      arrange(`Extreme models`)
  }
  
  
  # Pne lines
  if(is.null(x2)){
    df <- data.frame(Baseline = x1, Event = y)
    df <- df |> 
      filter(!is.na(Baseline)) |> 
      filter(!is.na(Event))
    df$ID = seq(1,nrow(df),1)
    
    n <- nrow(df)
    n_event = sum(df$Event)
    n_nonevent = n - n_event
    incidence <- 100 * n_event/n
    
    
    # Net Benefit for Decision Curves
    benefit <- df |> 
      mutate(Prediction = Baseline) |> 
      group_by(Prediction) |> 
      summarise(n_ev = sum(Event == 1),
                n_nonev = sum(Event == 0)) |> 
      ungroup() |> 
      mutate(TP = n_event - cumsum(n_ev),
             FN = n_event - TP, 
             TN = cumsum(n_nonev),
             FP = n_nonevent - TN,
             prevalence = (TP + FN)/n,
             sens = TP/(TP + FN),
             spec = TN/(TN + FP),
             `1-spec` = 1 - spec,
             npv = TN/(TN + FN),
             ppv = TP/(TP + FP),
             treated = TP/n - (FP/n) * Prediction/(1 - Prediction),
             untreated = TN/n - (FN/n) * (1 - Prediction)/Prediction,
             overall = treated + untreated,
             all = prevalence - (1 - prevalence)  * Prediction/(1 - Prediction),
             Model = "Baseline"
      )
    
    benefit_all_none <- benefit |> 
      select(Prediction, all) |> 
      mutate(none = 0)  |> 
      tidyr::pivot_longer(cols = c("all", "none"), names_to = "Extreme models", values_to = "extremes") |> 
      arrange(`Extreme models`)
  }
  
  g <- ggplot() + 
    scale_x_continuous(breaks = seq(0,1,0.1), expand = c(0.005,0.005)) + 
    scale_y_continuous( expand = c(0.005,0.005)) + 
    geom_point(data = benefit, aes(x = Prediction, y = treated, colour = Model), alpha = 0.3, size = 0.5)
  
  if (show_smooth) {
    g <- g + geom_smooth(data = benefit, aes(x = Prediction, y = treated, colour = Model), 
                         method = smooth_method, span = smooth_span, se = smooth_se)
  }
  
  g <- g + 
    geom_line(data = benefit_all_none, aes(x = Prediction, y = extremes, linetype = `Extreme models`)) +
    ylab("Net benefit") + xlab("Prediction threshold") +
    coord_cartesian(ylim = c(1.5 * min(benefit$treated), 1.05 * incidence/100)) +
    NULL
  return(g)
}


#' The ROC plot
#' 
#' ggroc plots Sensitivity v 1-Specificity
#' 
#' @param x1 Either a logistic regression fitted using glm (base package) or lrm (rms package) or alculated probabilities (eg through a logistic regression model) of the baseline model.  Must be between 0 & 1
#' @param x2 Either a logistic regression fitted using glm (base package) or lrm (rms package) or calculated probabilities (eg through a logistic regression model) of the new (alternative) model.   Must be between 0 & 1
#' @param y Binary of outcome of interest. Must be 0 or 1 (if fitted models are provided this is extracted from the fit which for an rms fit must have x = TRUE, y = TRUE). 
#' @param carrington_line The Useful Area is from the roc down to this line. It depends on prevalence and the costs of FP, FN, TP, TN.  Default is FALSE. See Carrington et al.  
#' @param costs Numeric vectors costs = c(cFP, cFN,cTP, cTN). The costs of FP, FN, TP, TN.  Default, c(0,0,1,1), is for there to be no costs for the FP & FN and identical costs for TN and TP.  See Carrington et al.  
#' @param label_number The number of points on the curve to label.The default has no labels.    
#' @return A ggplot object displaying the ROC curve(s) with sensitivity on the y-axis and 1-specificity on the x-axis. If two models are provided, both curves are shown for comparison.
#' @references  Carrington AM, Fieguth PW, Mayr F, James ND, Holzinger A, Pickering JW, et al. The ROC Diagonal is not Layperson's Chance: a New Baseline Shows the Useful Area. Machine Learning and Knowledge Extraction. Vienna, Austria: Springer; 2022. pp. 100-113. Available: 10.1007/978-3-031-14463-9_7.  
#' @import ggplot2
#' @import tidyr
#' @import dplyr
#' @import ggrepel
#' @export
#'
ggroc <- function(x1, x2=NULL, y=NULL,  carrington_line = FALSE, costs = c(0,0,1,1), label_number = NULL) {
  
  if (class(x1)[1] == "glm") {
    y = x1$y
    x1 = stats::predict(x1, type = "response")
    data_type = "glm"
  }
  if (class(x2)[1] == "glm") {
    x2 = stats::predict(x2, type = "response")
  }
  if (class(x1)[1] == "lrm") {
    if (length(x1$y) == 0 ) {
      stop("Fitted models with the rms package must be made using y = TRUE")
    }
    y = as.numeric(as.character(x1$y))
    x1 = stats::predict(x1, type = "fitted")
    data_type = "lrm"
  }
  if (class(x2)[1] == "lrm") {
    if (length(x2$y) == 0 ) {
      stop("Fitted models with the rms package must be made using y = TRUE")
    }
    x2 = stats::predict(x2, type = "fitted")
  }
  
  if (class(x1)[1] != "glm"  & class(x1)[1] != "lrm" ) {
    data_type = "User supplied"
  }
  if (!is.null(x2) & length(x1) != length(x2))
    stop("Reference (baseline) and New (Alt) model vectors must be the same length")
  if (is.null(y))
    stop("Oops - there must be event data (y)")
  
  
  if(!is.null(x2)) {
    df <- data.frame(Baseline = x1, New = x2, Event = y)
    df <- df |> 
      filter(!is.na(Baseline)) |> 
      filter(!is.na(New)) |> 
      filter(!is.na(Event))
  }
  if(is.null(x2)) {
    df <- data.frame(Baseline = x1, Event = y)
    df <- df |> 
      filter(!is.na(Baseline)) |> 
      filter(!is.na(Event))
  }
  
  
  df$ID = seq(1,nrow(df),1)
  
  n <- nrow(df)
  n_event = sum(df$Event)
  n_nonevent = n - n_event
  incidence <- 100 * n_event/n
  prevalence <- n_event/n
  incidence <- 100 * prevalence
  slope = ((1- prevalence)/prevalence)*(costs[1]-costs[4])/(costs[2]-costs[3])
  xintercept = 0.5*(1-1/slope)
  yintercept = 0.5*(1-slope)
  
  # roc curves
  roc.1 <- df |> 
    mutate(Prediction = Baseline) |> 
    group_by(Prediction) |> 
    summarise(n_ev = sum(Event == 1),
              n_nonev = sum(Event == 0)) |> 
    ungroup() |> 
    mutate(TP = n_event - cumsum(n_ev),
           FN = n_event - TP, 
           TN = cumsum(n_nonev),
           FP = n_nonevent - TN,
           prevalence = (TP + FN)/n,
           sens = TP/(TP + FN),
           spec = TN/(TN + FP),
           `1-spec` = 1 - spec,
           npv = TN/(TN + FN),
           ppv = TP/(TP + FP),
           treated = TP/n - (FP/n) * Prediction/(1 - Prediction),
           untreated = TN/n - (FN/n) * (1 - Prediction)/Prediction,
           overall = treated + untreated,
           all = prevalence - (1 - prevalence)  * Prediction/(1 - Prediction),
           Model = "Baseline"
    )
  
  roc_df <- roc.1
  
  if(!is.null(x2)){
    roc.2 <- df |> 
      mutate(Prediction = New) |> 
      group_by(Prediction) |> 
      summarise(n_ev = sum(Event == 1),
                n_nonev = sum(Event == 0)) |> 
      ungroup() |> 
      mutate(TP = n_event - cumsum(n_ev),
             FN = n_event - TP, 
             TN = cumsum(n_nonev),
             FP = n_nonevent - TN,
             prevalence = (TP + FN)/n,
             sens = TP/(TP + FN),
             spec = TN/(TN + FP),
             `1-spec` = 1 - spec,
             npv = TN/(TN + FN),
             ppv = TP/(TP + FP),
             treated = TP/n - (FP/n) * Prediction/(1 - Prediction),
             untreated = TN/n - (FN/n) * (1 - Prediction)/Prediction,
             overall = treated + untreated,
             all = prevalence - (1 - prevalence)  * Prediction/(1 - Prediction),
             Model = "New"
      )
    
    roc_df = bind_rows(roc.1, roc.2)
  }
  roc_df <- roc_df |> arrange(Model, sens, `1-spec` )
  
  g <- ggplot(data = roc_df) + 
    scale_x_continuous(breaks = seq(0,1,0.1), expand = c(0.005,0.005)) + 
    scale_y_continuous(breaks = seq(0,1,0.1), expand = c(0.005,0.005)) 
  
  df_polygon = data.frame(x = c(0,0 , (1-yintercept)/slope, -yintercept/slope), y = c(0,1, 1, 0))
  
  if(carrington_line == TRUE){
    g <- g + 
      geom_polygon(data = df_polygon, aes(x=x,y=y), fill = "grey60", alpha = 0.5) +
      geom_abline(aes(intercept = yintercept, slope = slope ),colour = "grey30", linetype = "dashed") +
      geom_point(data = data.frame(x=0.5, y=0.5), aes(x=x,y=y))
  } 
  
  g <- g +
    geom_step( aes(x = `1-spec`, y = sens, colour = Model)) + 
    geom_point(aes(x = `1-spec`, y = sens, colour = Model)) 
  
  if(!is.null(label_number)){
    labs <- roc_df |> 
      mutate(interval =cut_interval(Prediction, label_number, labels = FALSE)) |> 
      group_by(Model, interval) |> 
      mutate(rn = row_number()) |> 
      mutate(mid_rn = ceiling(max(rn)/2)) |> 
      filter(rn == mid_rn) |> 
      ungroup() |> 
      mutate(Prediction = as.character(round(Prediction,3)))
    
    g <- g + 
      geom_point(data = labs, aes(x=`1-spec`, y=sens, label=Prediction), shape = 22) +
      geom_text_repel(data = labs, aes(x=`1-spec`, y=sens, label=Prediction), min.segment.length = 0.1, colour = "black", nudge_x=0.05, nudge_y = 0.00) 
  } 
  
  g <- g + 
    ylab("Sensitivity") + xlab("1-Specificity") +
    coord_cartesian(ylim = c(0,1), xlim = c(0,1)) +
    NULL
  
  return(g)
}


#' The Precision-Recall plot
#' 
#' ggprerec plots Precision (PPV) v Recall (Sensitivity)
#' 
#' @param x1 Either a logistic regression fitted using glm (base package) or lrm (rms package) or alculated probabilities (eg through a logistic regression model) of the baseline model.  Must be between 0 & 1
#' @param x2 Either a logistic regression fitted using glm (base package) or lrm (rms package) or calculated probabilities (eg through a logistic regression model) of the new (alternative) model.   Must be between 0 & 1
#' @param y Binary of outcome of interest. Must be 0 or 1 (if fitted models are provided this is extracted from the fit which for an rms fit must have x = TRUE, y = TRUE). 
#' @param show_smooth Logical, whether to display smoothed curves (default = TRUE)
#' @param smooth_method Smoothing method for geom_smooth. Options: "loess", "lm", "glm", "gam". Default is "loess"
#' @param smooth_span Span parameter for loess smoothing, controls the degree of smoothing (default = 0.75). Lower values = less smooth
#' @param smooth_se Logical, whether to display confidence interval around smooth (default = FALSE)
#' @return A ggplot object displaying the precision-recall curve(s) with recall (sensitivity) on the x-axis and precision (positive predictive value) on the y-axis. If two models are provided, both curves are shown for comparison.
#' @import ggplot2
#' @import tidyr
#' @import dplyr
#' @export
#'
ggprerec <- function(x1, x2=NULL, y=NULL, show_smooth = TRUE, smooth_method = "loess", smooth_span = 0.75, smooth_se = FALSE) {
  
  if (class(x1)[1] == "glm") {
    y = x1$y
    x1 = stats::predict(x1, type = "response")
    data_type = "glm"
  }
  if (class(x2)[1] == "glm") {
    x2 = stats::predict(x2, type = "response")
  }
  if (class(x1)[1] == "lrm") {
    if (length(x1$y) == 0 ) {
      stop("Fitted models with the rms package must be made using y = TRUE")
    }
    y = as.numeric(as.character(x1$y))
    x1 = stats::predict(x1, type = "fitted")
    data_type = "lrm"
  }
  if (class(x2)[1] == "lrm") {
    if (length(x2$y) == 0 ) {
      stop("Fitted models with the rms package must be made using y = TRUE")
    }
    x2 = stats::predict(x2, type = "fitted")
  }
  
  if (class(x1)[1] != "glm"  & class(x1)[1] != "lrm" ) {
    data_type = "User supplied"
  }
  if (!is.null(x2) & length(x1) != length(x2))
    stop("Reference (baseline) and New (Alt) model vectors must be the same length")
  if (is.null(y))
    stop("Oops - there must be event data (y)")
  
  
  
  if(!is.null(x2)) {
    df <- data.frame(Baseline = x1, New = x2, Event = y)
    df <- df |> 
      filter(!is.na(Baseline)) |> 
      filter(!is.na(New)) |> 
      filter(!is.na(Event))
  }
  if(is.null(x2)) {
    df <- data.frame(Baseline = x1, Event = y)
    df <- df |> 
      filter(!is.na(Baseline)) |> 
      filter(!is.na(Event))
  }
  
  
  df$ID = seq(1,nrow(df),1)
  
  n <- nrow(df)
  n_event = sum(df$Event)
  n_nonevent = n - n_event
  incidence <- 100 * n_event/n
  prevalence <- n_event/n
  incidence <- 100 * prevalence
  
  
  # pr curves
  pr.1 <- df |> 
    mutate(Prediction = Baseline) |> 
    group_by(Prediction) |> 
    summarise(n_ev = sum(Event == 1),
              n_nonev = sum(Event == 0)) |> 
    ungroup() |> 
    mutate(TP = n_event - cumsum(n_ev),
           FN = n_event - TP, 
           TN = cumsum(n_nonev),
           FP = n_nonevent - TN,
           prevalence = (TP + FN)/n,
           sens = TP/(TP + FN),
           spec = TN/(TN + FP),
           `1-spec` = 1 - spec,
           npv = TN/(TN + FN),
           ppv = TP/(TP + FP),
           treated = TP/n - (FP/n) * Prediction/(1 - Prediction),
           untreated = TN/n - (FN/n) * (1 - Prediction)/Prediction,
           overall = treated + untreated,
           all = prevalence - (1 - prevalence)  * Prediction/(1 - Prediction),
           Model = "Baseline"
    )
  
  pr_df <- pr.1
  
  if(!is.null(x2)){
    pr.2 <- df |> 
      mutate(Prediction = New) |> 
      group_by(Prediction) |> 
      summarise(n_ev = sum(Event == 1),
                n_nonev = sum(Event == 0)) |> 
      ungroup() |> 
      mutate(TP = n_event - cumsum(n_ev),
             FN = n_event - TP, 
             TN = cumsum(n_nonev),
             FP = n_nonevent - TN,
             prevalence = (TP + FN)/n,
             sens = TP/(TP + FN),
             spec = TN/(TN + FP),
             `1-spec` = 1 - spec,
             npv = TN/(TN + FN),
             ppv = TP/(TP + FP),
             treated = TP/n - (FP/n) * Prediction/(1 - Prediction),
             untreated = TN/n - (FN/n) * (1 - Prediction)/Prediction,
             overall = treated + untreated,
             all = prevalence - (1 - prevalence)  * Prediction/(1 - Prediction),
             Model = "New"
      )
    
    pr_df = bind_rows(pr.1, pr.2)
  }
  
  g <- ggplot(data = pr_df) + 
    scale_x_continuous(breaks = seq(0,1,0.1), expand = c(0.005,0.005)) + 
    scale_y_continuous(breaks = seq(0,1,0.1), expand = c(0.005,0.005))  +
    geom_point(aes(x = sens, y = ppv, colour = Model), alpha = 0.3, size = 0.5)
  
  if (show_smooth) {
    g <- g + geom_smooth(aes(x = sens, y = ppv, colour = Model), 
                         method = smooth_method, span = smooth_span, se = smooth_se)
  }
  
  g <- g + 
    ylab("Precision (PPV)") + xlab("Recall (Sensitivity)") +
    coord_cartesian(ylim = c(0,1), xlim = c(0,1)) +
    NULL
  
  return(g)
}

#' Prepare calibration inputs
#'
#' Extracts predictions and outcomes from fitted glm/lrm models (or uses the
#' supplied vectors), checks them and drops incomplete cases. Used by
#' ggcalibrate and ggcalibrate_BA.
#' @noRd
calibration_inputs <- function(x1, x2 = NULL, y = NULL) {

  if (class(x1)[1] == "glm") {
    y = x1$y
    x1 = stats::predict(x1, type = "response")
  }
  if (class(x2)[1] == "glm") {
    x2 = stats::predict(x2, type = "response")
  }
  if (class(x1)[1] == "lrm") {
    if (length(x1$y) == 0 ) {
      stop("Fitted models with the rms package must be made using y = TRUE")
    }
    y = as.numeric(as.character(x1$y))
    x1 = stats::predict(x1, type = "fitted")
  }
  if (class(x2)[1] == "lrm") {
    if (length(x2$y) == 0 ) {
      stop("Fitted models with the rms package must be made using y = TRUE")
    }
    x2 = stats::predict(x2, type = "fitted")
  }

  if (!is.null(x2) & length(x1) != length(x2))
    stop("Reference (baseline) and New (alternative) model vectors must be the same length")
  if (is.null(y))
    stop("Oops - there must be event data (y)")

  # Keep complete cases only, so both curves are fitted to the same patients
  if (is.null(x2)) {
    df <- data.frame(Baseline = x1, Event = y)
  } else {
    df <- data.frame(Baseline = x1, New = x2, Event = y)
  }
  df[stats::complete.cases(df), , drop = FALSE]
}

#' Calibration curve with confidence interval
#'
#' Fits a logistic regression of the outcome on a restricted cubic spline of
#' the prediction. The confidence interval is calculated on the log-odds scale
#' and back-transformed, so it always lies between 0 and 1.
#' @noRd
calibration_curve <- function(prediction, event, model, n_knots, ci_level) {

  fit <- glm(event ~ rms::rcs(prediction, n_knots), family = "binomial")
  lp <- stats::predict(fit, type = "link", se.fit = TRUE)
  z <- stats::qnorm(ci_level + (1 - ci_level)/2)

  data.frame(Model = model,
             x = prediction,
             y = stats::plogis(lp$fit),
             lci = stats::plogis(lp$fit - z * lp$se.fit),
             uci = stats::plogis(lp$fit + z * lp$se.fit),
             Event = event)
}

#' Calibration curves for the baseline and (optionally) new model
#' @noRd
calibration_curves <- function(df, n_knots, ci_level) {

  curves <- calibration_curve(df$Baseline, df$Event, "Baseline", n_knots, ci_level)
  if ("New" %in% names(df)) {
    curves <- bind_rows(curves,
                        calibration_curve(df$New, df$Event, "New", n_knots, ci_level))
  }
  curves
}

#' The Calibration plot
#'
#' ggcalibrate plots the stats::predicted events against the actual event rate
#'
#' The calibration curve for each model is a logistic regression of the outcome on a restricted cubic spline of the predicted probability. The confidence interval is calculated from the standard error of that fit on the log-odds scale and back-transformed, so it is bounded by 0 and 1.
#'
#' @param x1 Either a logistic regression fitted using glm (base package) or lrm (rms package) or calculated probabilities (eg through a logistic regression model) of the baseline model.  Must be between 0 & 1
#' @param x2 Either a logistic regression fitted using glm (base package) or lrm (rms package) or calculated probabilities (eg through a logistic regression model) of the new (alternative) model.   Must be between 0 & 1
#' @param y Binary of outcome of interest. Must be 0 or 1 (if fitted models are provided this is extracted from the fit which for an rms fit must have x = TRUE, y = TRUE).
#' @param n_knots The curves are made by fitting a restricted cubic spline (rms package). The default 5-knots is usually enough.
#' @param ci_level Confidence interval of the curve (default = 0.95).
#' @param alpha_level Transparency (alpha) of the shaded confidence interval (default = 0.25).
#' @param actuals Logical, whether to also plot the actual events (0 or 1) against the predictions as short vertical marks, like a rug plot, at 0 and 1 (default = FALSE).
#' @param smooth_method Deprecated and ignored. The curves are no longer drawn with geom_smooth().
#' @param smooth_span Deprecated and ignored. The curves are no longer drawn with geom_smooth().
#' @return a ggplot
#' @seealso [ggcalibrate_BA()] for the same curves plotted as deviations from perfect calibration.
#' @examples
#' # Quick example with subset of data
#' data(data_risk)
#' data_subset <- data_risk[1:100, ]  # Use first 100 rows for speed
#' complete_cases <- complete.cases(data_subset)
#' data_clean <- data_subset[complete_cases, ]
#' y <- data_clean$outcome
#' x1 <- data_clean$baseline
#' x2 <- data_clean$new
#' output <- ggcalibrate(x1, x2, y, n_knots = 3, ci_level = 0.95)
#'
#' \donttest{
#' # Full dataset example
#' data(data_risk)
#' complete_cases <- complete.cases(data_risk)
#' data_clean <- data_risk[complete_cases, ]
#' y <- data_clean$outcome
#' x1 <- data_clean$baseline
#' x2 <- data_clean$new
#' output <- ggcalibrate(x1, x2, y, n_knots = 5, ci_level = 0.95)
#'
#' # Show the actual events and a darker confidence interval
#' output <- ggcalibrate(x1, x2, y, alpha_level = 0.5, actuals = TRUE)
#' }
#' @import forcats
#' @import ggplot2
#' @import tidyr
#' @import dplyr
#' @importFrom pracma trapz
#' @export
ggcalibrate <- function(x1, x2 = NULL, y = NULL,  n_knots = 5, ci_level = 0.95, alpha_level = 0.25, actuals = FALSE, smooth_method = NULL, smooth_span = NULL) {

  if (!is.null(smooth_method) | !is.null(smooth_span))
    warning("`smooth_method` and `smooth_span` are deprecated and ignored: ggcalibrate() no longer uses geom_smooth().", call. = FALSE)

  df <- calibration_inputs(x1, x2, y)
  df_calib <- calibration_curves(df, n_knots, ci_level)

  # Plot
  g <- ggplot(data = df_calib, aes(x=x,y=y, colour = Model )) +
    scale_x_continuous(breaks = seq(0,1,0.1), expand = c(0.005,0.005)) +
    scale_y_continuous(breaks = seq(0,1,0.1), expand = c(0.005,0.005)) +
    geom_abline(slope = 1, intercept = 0, colour = "grey50", linetype = "dashed")  +
    geom_ribbon(aes(ymin = lci, ymax = uci, fill = Model), colour = NA, alpha = alpha_level) +
    geom_line() +
    xlab("Predicted probability") +
    ylab("Actual probability") +
    coord_cartesian(xlim = c(0,1), ylim = c(0,1)) +
    NULL

  if (actuals) {
    g <- g +
      geom_point(aes(y = Event), shape = "I", alpha = 0.5)
  }

  return(g)
}

#' The Bland-Altman style Calibration plot
#'
#' ggcalibrate_BA plots the deviation of the actual event rate from the stats::predicted events (Actual - Prediction) against the prediction, similar to a Bland-Altman plot.
#'
#' Perfect calibration is the horizontal line at zero. A curve above the line means the model under-predicts the outcome, and a curve below it means the model over-predicts. The curves and confidence intervals are the same as in [ggcalibrate()], with the prediction subtracted. To focus on the region of interest (eg around a decision threshold), add `coord_cartesian(xlim = ...)` to the plot.
#'
#' @param x1 Either a logistic regression fitted using glm (base package) or lrm (rms package) or calculated probabilities (eg through a logistic regression model) of the baseline model.  Must be between 0 & 1
#' @param x2 Either a logistic regression fitted using glm (base package) or lrm (rms package) or calculated probabilities (eg through a logistic regression model) of the new (alternative) model.   Must be between 0 & 1
#' @param y Binary of outcome of interest. Must be 0 or 1 (if fitted models are provided this is extracted from the fit which for an rms fit must have x = TRUE, y = TRUE).
#' @param n_knots The curves are made by fitting a restricted cubic spline (rms package). The default 5-knots is usually enough.
#' @param ci_level Confidence interval of the curve (default = 0.95).
#' @param alpha_level Transparency (alpha) of the shaded confidence interval (default = 0.25).
#' @return a ggplot
#' @seealso [ggcalibrate()]
#' @examples
#' # Quick example with subset of data
#' data(data_risk)
#' data_subset <- data_risk[1:100, ]  # Use first 100 rows for speed
#' complete_cases <- complete.cases(data_subset)
#' data_clean <- data_subset[complete_cases, ]
#' y <- data_clean$outcome
#' x1 <- data_clean$baseline
#' x2 <- data_clean$new
#' output <- ggcalibrate_BA(x1, x2, y, n_knots = 3, ci_level = 0.95, alpha_level = 0.25)
#'
#' \donttest{
#' # Full dataset example
#' data(data_risk)
#' complete_cases <- complete.cases(data_risk)
#' data_clean <- data_risk[complete_cases, ]
#' y <- data_clean$outcome
#' x1 <- data_clean$baseline
#' x2 <- data_clean$new
#' output <- ggcalibrate_BA(x1, x2, y, n_knots = 5, ci_level = 0.95, alpha_level = 0.5)
#'
#' # Zoom in on predictions below 30%
#' output + ggplot2::coord_cartesian(xlim = c(0, 0.3), ylim = c(-0.3, 0.3))
#' }
#' @import ggplot2
#' @import dplyr
#' @export
ggcalibrate_BA <- function(x1, x2 = NULL, y = NULL, n_knots = 5, ci_level = 0.95, alpha_level = 0.25) {

  df <- calibration_inputs(x1, x2, y)
  df_deviation <- calibration_curves(df, n_knots, ci_level) |>
    mutate(y = y - x, lci = lci - x, uci = uci - x)

  # Plot
  g <- ggplot(data = df_deviation, aes(x=x,y=y, colour = Model )) +
    geom_hline(yintercept = 0, colour = "grey50", linetype = "dashed") +
    geom_ribbon(aes(ymin = lci, ymax = uci, fill = Model), colour = NA, alpha = alpha_level) +
    geom_line() +
    scale_x_continuous(breaks = seq(0,1,0.1), expand = c(0.005,0.005)) +
    scale_y_continuous(breaks = seq(-1,1,0.1)) +
    coord_cartesian(xlim = c(0,1), ylim = c(-1,1)) +
    ylab("Deviation from the outcome \n(Actual - Prediction)") +
    xlab("Prediction") +
    NULL

  return(g)
}



#' The Original Calibration plot
#' 
#' ggcalibrate_original plots the stats::predicted events against the actual event rate using the "old" form.  
#' 
#' @param x1 Either a logistic regression fitted using glm (base package) or lrm (rms package) or calculated probabilities (eg through a logistic regression model) of the baseline model.  Must be between 0 & 1
#' @param x2 Either a logistic regression fitted using glm (base package) or lrm (rms package) or calculated probabilities (eg through a logistic regression model) of the new (alternative) model.   Must be between 0 & 1
#' @param y Binary of outcome of interest. Must be 0 or 1 (if fitted models are provided this is extracted from the fit which for an rms fit must have x = TRUE, y = TRUE). 
#' @param n_cut An integer indicating either the number of intervals of the same width, the number of intervals of the same number of subjects, or the width (as a percentage) of the intervals.  
#' @param cut_type One of three strings: "interval", "number", or "width". 
#'   - "interval": uses cut_interval() to get n_cut intervals of approximately equal width.
#'   - "number": uses cut_number() to get n_cut intervals with approximately equal counts.
#'   - "width": uses cut_width() to get intervals of a fixed width (approximately 100/n_cut).
#' @param include_margin TRUE for including producing a bar plot of the counts of in each of the intervals. Default is FALSE.  Note if the output is saved to my_graphs then using the library gridExtra the function grid.arrange(graphs$g, graphs$g_marg , nrow = 2, heights = c(2,1)) will produce a plot with both the calibration plot and the marginal plot.  
#' @return a list of one or two ggplots
#' @examples
#' # Quick example with subset of data
#' data(data_risk)
#' data_subset <- data_risk[1:100, ]  # Use first 100 rows for speed
#' complete_cases <- complete.cases(data_subset)
#' data_clean <- data_subset[complete_cases, ]
#' y <- data_clean$outcome 
#' x1 <- data_clean$baseline
#' x2 <- data_clean$new
#' output <- ggcalibrate_original(
#'   x1, x2, y,
#'   n_cut = 3, cut_type = "interval",
#'   include_margin = FALSE
#' )
#' 
#' \donttest{
#' # Full dataset example
#' data(data_risk)
#' complete_cases <- complete.cases(data_risk)
#' data_clean <- data_risk[complete_cases, ]
#' y <- data_clean$outcome 
#' x1 <- data_clean$baseline
#' x2 <- data_clean$new
#' output <- ggcalibrate_original(
#'   x1, x2, y,
#'   n_cut = 5, cut_type = "interval",
#'   include_margin = FALSE
#' )
#' }
#' @import forcats
#' @import ggplot2
#' @import tidyr
#' @import dplyr
#' @importFrom pracma trapz
#' @export
ggcalibrate_original <- function(x1, x2 = NULL, y = NULL, n_cut = 5, 
                                 cut_type = c("interval","number","width"), include_margin = FALSE) {
  
  # Default to "interval" if cut_type not specified
  if (length(cut_type) > 1) {
    cut_type <- "interval"
  }
  
  if (class(x1)[1] == "glm") {
    y = x1$y
    x1 = stats::predict(x1, type = "response")
    data_type = "glm"
  }
  if (class(x2)[1] == "glm") {
    x2 = stats::predict(x2, type = "response")
  }
  if (class(x1)[1] == "lrm") {
    if (length(x1$y) == 0 ) {
      stop("Fitted models with the rms package must be made using y = TRUE")
    }
    y = as.numeric(as.character(x1$y))
    x1 = stats::predict(x1, type = "fitted")
    data_type = "lrm"
  }
  if (class(x2)[1] == "lrm") {
    if (length(x2$y) == 0 ) {
      stop("Fitted models with the rms package must be made using y = TRUE")
    }
    x2 = stats::predict(x2, type = "fitted")
  }
  
  if (class(x1)[1] != "glm"  & class(x1)[1] != "lrm" ) {
    data_type = "User supplied"
  }
  if (!is.null(x2) & length(x1) != length(x2))
    stop("Reference (baseline) and New (Alt) model vectors must be the same length")
  if (is.null(y))
    stop("Oops - there must be event data (y)")
  
  
  if (!is.null(x2)){
    df <- data.frame(Baseline = x1, New = x2, Event = y)
  }
  if (is.null(x2)){
    df <- data.frame(Baseline = x1, Event = y)
  }
  
  df$ID = seq(1,nrow(df),1)
  
  n <- nrow(df)
  n_event = sum(df$Event)
  n_nonevent = n - n_event
  incidence <- 100 * n_event/n
  
  
  # Calibration
  if (!is.null(x2)){
    df_long <- df %>% 
      tidyr::pivot_longer(cols = c(Baseline,New), values_to = "prediction", names_to = "Model") %>% 
      select(Model,Event, prediction) %>% # or whatever your outcome and predictions are
      mutate(prediction = 100 * as.numeric(prediction)) %>%  # If prediciton is in the 0-1 range I prefer it in the 0-100 range
      filter(!is.na(prediction)) 
  }
  
  if (is.null(x2)){
    df_long <- df %>% 
      tidyr::pivot_longer(cols = c(Baseline), values_to = "prediction", names_to = "Model") %>% 
      select(Model,Event, prediction) %>% # or whatever your outcome and predictions are
      mutate(prediction = 100 * as.numeric(prediction)) %>%  # If prediciton is in the 0-1 range I prefer it in the 0-100 range
      filter(!is.na(prediction)) 
  }
  
  
  if(cut_type == "interval"){df_long <- df_long %>% group_by(Model) %>% mutate(pred_interval = cut_interval(prediction, n = n_cut))}
  if(cut_type == "number"){df_long <- df_long %>% group_by(Model) %>% mutate(pred_interval = cut_number(prediction, n = n_cut))}
  if(cut_type == "width"){df_long <- df_long %>% group_by(Model) %>% mutate(pred_interval = cut_width(prediction, width = n_cut))}
  
  # prepare file to use in plot - sometimes only the CIs are made for the actual events and not the prediction, but that is just laziness
  gdf <- df_long %>%
    group_by(Model, pred_interval) %>%
    summarise(n = n(),
              mn_pred = mean(prediction, na.rm = TRUE),
              lci_pred = quantile(prediction, c(0.025), na.rm = TRUE),
              uci_pred = quantile(prediction, c(0.975), na.rm = TRUE),
              n_event = sum(Event == "1")) %>%
    mutate(prop_event = 100 * n_event/n) %>%
    mutate(lci_prop_event = prop_event - 1.96 * sqrt(prop_event * (100 - prop_event)/n)) %>%
    mutate(uci_prop_event = prop_event + 1.96 * sqrt(prop_event * (100 - prop_event)/n))
  
  g <- ggplot() +
    geom_abline(slope = 1, intercept = 0, colour = "grey50", linetype = "dashed") 
  
  if(!is.null(x2)){
    g <- g +
      geom_segment(data = gdf, aes(x = mn_pred, xend = mn_pred, y = lci_prop_event, yend = uci_prop_event, colour = Model))  +
      geom_segment(data = gdf, aes(x = lci_pred, xend = uci_pred, y = prop_event, yend = prop_event, colour = Model)) +
      geom_point(data = gdf, aes(x = mn_pred, y = prop_event, colour = Model))
    
  }
  if(is.null(x2)){
    g <- g +
      geom_segment(data = gdf, aes(x = mn_pred, xend = mn_pred, y = lci_prop_event, yend = uci_prop_event))  +
      geom_segment(data = gdf, aes(x = lci_pred, xend = uci_pred, y = prop_event, yend = prop_event)) +
      geom_point(data = gdf, aes(x = mn_pred, y = prop_event))
  }
  
  g <- g +
    scale_x_continuous(breaks = seq(0,100,10)) +
    scale_y_continuous(breaks = seq(0,100,10)) +
    xlab("Predicted percentage") +
    ylab("Actual percentage") +
    coord_cartesian(xlim = c(0,100), ylim = c(0,100)) +
    NULL
  
  ifelse(include_margin == TRUE,
         {
           gdf <- gdf %>% 
             group_by(Model) %>% 
             arrange(mn_pred) %>% 
             mutate(width_ci = uci_pred - lci_pred)
           width = mean(gdf$width_ci, na.rm = TRUE)
           gdf$pred_interval <- fct_inorder(gdf$pred_interval)
           
           
           if(!is.null(x2)){ 
             g_marg <- ggplot(gdf) +
               geom_col(aes(x = mn_pred, y = n, fill = Model), width = width * 2/3, alpha = 0.6)  
           }
           
           if(is.null(x2)){
             g_marg <- ggplot(gdf) +
               geom_col(aes(x = mn_pred, y = n), width = width * 2/3, alpha = 0.6)  
           }
           
           g_marg <- g_marg +   
             scale_x_continuous(breaks = seq(0,100,10)) +
             xlab("Predicted percentage") +
             ylab("Count") +
             coord_cartesian(xlim = c(0,100)) +
             NULL
           
           output <- list(g = g, g_marg = g_marg)
         },
         {
           output <- list(g = g)
         }
  )  
  
  return(output)
}

#'
#' The function ggrap() plots the Sensitivity and 1-Specificity curves against the calculated risk for the baseline (reference) and newmodels, thus graphically displaying the IDIs for those with and without the events.  These plots can aid interpretation of the NRI and IDI metrics.
#' 
#' @param x1 Either a logistic regression fitted using glm (base package) or lrm (rms package) or alculated probabilities (eg through a logistic regression model) of the baseline model.  Must be between 0 & 1
#' @param x2 Either a logistic regression fitted using glm (base package) or lrm (rms package) or calculated probabilities (eg through a logistic regression model) of the new (alternative) model.   Must be between 0 & 1
#' @param y Binary of outcome of interest. Must be 0 or 1 (if fitted models are provided this is extracted from the fit which for an rms fit must have x = TRUE, y = TRUE). 
#' @return a ggplot 
#' @references The Risk Assessment Plot in this form was described by Pickering, J. W., & Endre, Z. H. (2012). New Metrics for Assessing Diagnostic Potential of Candidate Biomarkers. Clinical Journal of the American Society of Nephrology, 7, 1355–1364. doi:10.2215/CJN.09590911
#' @import ggplot2
#' @import tidyr 
#' @import dplyr
#' @export
#'
ggrap <- function(x1, x2=NULL, y=NULL) {
  
  if (class(x1)[1] == "glm") {
    y = x1$y
    x1 = stats::predict(x1, type = "response")
    data_type = "glm"
  }
  if (class(x2)[1] == "glm") {
    x2 = stats::predict(x2, type = "response")
  }
  if (class(x1)[1] == "lrm") {
    if (length(x1$y) == 0 ) {
      stop("Fitted models with the rms package must be made using y = TRUE")
    }
    y = as.numeric(as.character(x1$y))
    x1 = stats::predict(x1, type = "fitted")
    data_type = "lrm"
  }
  if (class(x2)[1] == "lrm") {
    if (length(x2$y) == 0 ) {
      stop("Fitted models with the rms package must be made using y = TRUE")
    }
    x2 = stats::predict(x2, type = "fitted")
  }
  
  if (class(x1)[1] != "glm"  & class(x1)[1] != "lrm" ) {
    data_type = "User supplied"
  }
  if (!is.null(x2) & length(x1) != length(x2))
    stop("Reference (baseline) and New (Alt) model vectors must be the same length")
  
  
  if(!is.null(x2)){
    df <- data.frame(Baseline = x1, New = x2, Event = y)
    df <- df |> 
      filter(!is.na(Baseline)) |> 
      filter(!is.na(New)) |> 
      filter(!is.na(Event))
    df$ID = seq(1,nrow(df),1)
    
    df_long <- df |> 
      tidyr::pivot_longer(cols = c(Baseline,New), values_to = "Probabilities", names_to = "Model") |> 
      group_by(Model)  |> 
      arrange(Model, Probabilities) |> 
      group_by(Model,Probabilities) |> 
      reframe(n_event = sum(Event == 1),
              n_nonevent = sum(Event == 0),
              IDs = as.character(ID)) |> 
      ungroup() |> 
      group_by(Model) |> 
      mutate(n_event_cumulative = cumsum(n_event)) |> 
      mutate(n_nonevent_cumulative = cumsum(n_nonevent)) |> 
      mutate(Sensitivity = (max(n_event_cumulative) - n_event_cumulative)/max(n_event_cumulative)) |> 
      mutate(`1-Specificity` = (max(n_nonevent_cumulative) - n_nonevent_cumulative)/max(n_nonevent_cumulative)) |> 
      ungroup()
  }
  
  if(is.null(x2)){
    df <- data.frame(Baseline = x1,  Event = y)
    df <- df |> 
      filter(!is.na(Baseline)) |> 
      filter(!is.na(Event))
    df$ID = seq(1,nrow(df),1)
    
    df_long <- df |> 
      tidyr::pivot_longer(cols = c(Baseline), values_to = "Probabilities", names_to = "Model") |> 
      group_by(Model)  |> 
      arrange(Model, Probabilities) |> 
      group_by(Model,Probabilities) |> 
      reframe(n_event = sum(Event == 1),
              n_nonevent = sum(Event == 0),
              IDs = as.character(ID)) |> 
      ungroup() |> 
      group_by(Model) |> 
      mutate(n_event_cumulative = cumsum(n_event)) |> 
      mutate(n_nonevent_cumulative = cumsum(n_nonevent)) |> 
      mutate(Sensitivity = (max(n_event_cumulative) - n_event_cumulative)/max(n_event_cumulative)) |> 
      mutate(`1-Specificity` = (max(n_nonevent_cumulative) - n_nonevent_cumulative)/max(n_nonevent_cumulative)) |> 
      ungroup()
  }
  
  df_long_av <- df_long |> 
    group_by(Model, Probabilities) |> 
    summarise(Sensitivity = mean(Sensitivity),
              `1-Specificity` = mean(`1-Specificity`),
              ID_list = list(IDs)
    ) |> 
    ungroup()
  
  df_g <- df_long_av |> 
    select(Probabilities, Sensitivity, `1-Specificity`, Model, ID_list) |>
    tidyr::pivot_longer(cols = c(Sensitivity, `1-Specificity`), values_to = "Metric", names_to = "Metric name")
  
  
  if(!is.null(x2)){
    g <- ggplot(df_g, aes(x = Probabilities, y = Metric, colour = `Metric name`, linetype = Model )) + 
      scale_x_continuous(breaks = seq(0,1,0.1), expand = c(0.005,0.005)) + 
      scale_y_continuous(breaks = seq(0,1,0.1), expand = c(0.005,0.005)) +
      geom_line() +
      scale_linetype_manual(values = c("dotted", "solid")) 
  }
  if(is.null(x2)){
    g <- ggplot(df_g, aes(x = Probabilities, y = Metric, colour = `Metric name` )) + 
      scale_x_continuous(breaks = seq(0,1,0.1), expand = c(0.005,0.005)) + 
      scale_y_continuous(breaks = seq(0,1,0.1), expand = c(0.005,0.005)) +
      geom_line() 
  }
  
  g <- g +
    coord_cartesian(xlim = c(0,1), ylim = c(0,1))    
  
  return(g)
}

#' The Decision curve
#' 
#' ggdecision plots decision curves to assess the net benefit at different thresholds
#' 
#' @param x1 Either a logistic regression fitted using glm (base package) or lrm (rms package) or calculated probabilities (eg through a logistic regression model) of the baseline model.  Must be between 0 & 1
#' @param x2 Either a logistic regression fitted using glm (base package) or lrm (rms package) or calculated probabilities (eg through a logistic regression model) of the new (alternative) model.   Must be between 0 & 1
#' @param y Binary of outcome of interest. Must be 0 or 1 (if fitted models are provided this is extracted from the fit which for an rms fit must have x = TRUE, y = TRUE). 
#' @param show_smooth Logical, whether to display smoothed curves (default = TRUE)
#' @param smooth_method Smoothing method for geom_smooth. Options: "loess", "lm", "glm", "gam". Default is "loess"
#' @param smooth_span Span parameter for loess smoothing, controls the degree of smoothing (default = 0.75). Lower values = less smooth
#' @param smooth_se Logical, whether to display confidence interval around smooth (default = FALSE)
#' @return a ggplot 
#' @references Vickers AJ, van Calster B, Steyerberg EW. A simple, step-by-step guide to interpreting decision curve analysis. Diagn Progn Res 2019;3(1):18. 2. Zhang Z, Rousson V, Lee W-C, et al. Decision curve analysis: a technical note. Ann Transl Med 2018;6(15):308–308. 
#' @import ggplot2
#' @import tidyr
#' @import dplyr
#' @export
#'
ggdecision <- function(x1, x2=NULL, y=NULL, show_smooth = TRUE, smooth_method = "loess", smooth_span = 0.75, smooth_se = FALSE) {
  
  if (class(x1)[1] == "glm") {
    y = x1$y
    x1 = stats::predict(x1, type = "response")
    data_type = "glm"
  }
  if (class(x2)[1] == "glm") {
    x2 = stats::predict(x2, type = "response")
  }
  if (class(x1)[1] == "lrm") {
    if (length(x1$y) == 0 ) {
      stop("Fitted models with the rms package must be made using y = TRUE")
    }
    y = as.numeric(as.character(x1$y))
    x1 = stats::predict(x1, type = "fitted")
    data_type = "lrm"
  }
  if (class(x2)[1] == "lrm") {
    if (length(x2$y) == 0 ) {
      stop("Fitted models with the rms package must be made using y = TRUE")
    }
    x2 = stats::predict(x2, type = "fitted")
  }
  
  if (class(x1)[1] != "glm"  & class(x1)[1] != "lrm" ) {
    data_type = "User supplied"
  }
  if (!is.null(x2) & length(x1) != length(x2))
    stop("Reference (baseline) and New (Alt) model vectors must be the same length")
  
  # Two lines
  if(!is.null(x2)){
    df <- data.frame(Baseline = x1, New = x2, Event = y)
    df <- df |> 
      filter(!is.na(Baseline)) |> 
      filter(!is.na(New)) |> 
      filter(!is.na(Event))
    df$ID = seq(1,nrow(df),1)
    
    n <- nrow(df)
    n_event = sum(df$Event)
    n_nonevent = n - n_event
    incidence <- 100 * n_event/n
    
    
    # Net Benefit for Decision Curves
    benefit.1 <- df |> 
      mutate(Prediction = Baseline) |> 
      group_by(Prediction) |> 
      summarise(n_ev = sum(Event == 1),
                n_nonev = sum(Event == 0)) |> 
      ungroup() |> 
      mutate(TP = n_event - cumsum(n_ev),
             FN = n_event - TP, 
             TN = cumsum(n_nonev),
             FP = n_nonevent - TN,
             prevalence = (TP + FN)/n,
             sens = TP/(TP + FN),
             spec = TN/(TN + FP),
             `1-spec` = 1 - spec,
             npv = TN/(TN + FN),
             ppv = TP/(TP + FP),
             treated = TP/n - (FP/n) * Prediction/(1 - Prediction),
             untreated = TN/n - (FN/n) * (1 - Prediction)/Prediction,
             overall = treated + untreated,
             all = prevalence - (1 - prevalence)  * Prediction/(1 - Prediction),
             Model = "Baseline"
      )
    
    benefit.2 <- df |> 
      mutate(Prediction = New) |> 
      group_by(Prediction) |> 
      summarise(n_ev = sum(Event == 1),
                n_nonev = sum(Event == 0)) |> 
      ungroup() |> 
      mutate(TP = n_event - cumsum(n_ev),
             FN = n_event - TP, 
             TN = cumsum(n_nonev),
             FP = n_nonevent - TN,
             prevalence = (TP + FN)/n,
             sens = TP/(TP + FN),
             spec = TN/(TN + FP),
             `1-spec` = 1 - spec,
             npv = TN/(TN + FN),
             ppv = TP/(TP + FP),
             treated = TP/n - (FP/n) * Prediction/(1 - Prediction),
             untreated = TN/n - (FN/n) * (1 - Prediction)/Prediction,
             overall = treated + untreated,
             all = prevalence - (1 - prevalence)  * Prediction/(1 - Prediction),
             Model = "New"
      )
    
    benefit = bind_rows(benefit.1, benefit.2)
    
    benefit_all_none <- benefit |> 
      select(Prediction, all) |> 
      mutate(none = 0)  |> 
      tidyr::pivot_longer(cols = c("all", "none"), names_to = "Extreme models", values_to = "extremes") |> 
      arrange(`Extreme models`)
  }
  
  
  # Pne lines
  if(is.null(x2)){
    df <- data.frame(Baseline = x1, Event = y)
    df <- df |> 
      filter(!is.na(Baseline)) |> 
      filter(!is.na(Event))
    df$ID = seq(1,nrow(df),1)
    
    n <- nrow(df)
    n_event = sum(df$Event)
    n_nonevent = n - n_event
    incidence <- 100 * n_event/n
    
    
    # Net Benefit for Decision Curves
    benefit <- df |> 
      mutate(Prediction = Baseline) |> 
      group_by(Prediction) |> 
      summarise(n_ev = sum(Event == 1),
                n_nonev = sum(Event == 0)) |> 
      ungroup() |> 
      mutate(TP = n_event - cumsum(n_ev),
             FN = n_event - TP, 
             TN = cumsum(n_nonev),
             FP = n_nonevent - TN,
             prevalence = (TP + FN)/n,
             sens = TP/(TP + FN),
             spec = TN/(TN + FP),
             `1-spec` = 1 - spec,
             npv = TN/(TN + FN),
             ppv = TP/(TP + FP),
             treated = TP/n - (FP/n) * Prediction/(1 - Prediction),
             untreated = TN/n - (FN/n) * (1 - Prediction)/Prediction,
             overall = treated + untreated,
             all = prevalence - (1 - prevalence)  * Prediction/(1 - Prediction),
             Model = "Baseline"
      )
    
    benefit_all_none <- benefit |> 
      select(Prediction, all) |> 
      mutate(none = 0)  |> 
      tidyr::pivot_longer(cols = c("all", "none"), names_to = "Extreme models", values_to = "extremes") |> 
      arrange(`Extreme models`)
  }
  
  g <- ggplot() + 
    scale_x_continuous(breaks = seq(0,1,0.1), expand = c(0.005,0.005)) + 
    scale_y_continuous( expand = c(0.005,0.005)) + 
    geom_point(data = benefit, aes(x = Prediction, y = treated, colour = Model), alpha = 0.3, size = 0.5)
  
  if (show_smooth) {
    g <- g + geom_smooth(data = benefit, aes(x = Prediction, y = treated, colour = Model), 
                         method = smooth_method, span = smooth_span, se = smooth_se)
  }
  
  g <- g + 
    geom_line(data = benefit_all_none, aes(x = Prediction, y = extremes, linetype = `Extreme models`)) +
    ylab("Net benefit") + xlab("Prediction threshold") +
    coord_cartesian(ylim = c(1.5 * min(benefit$treated), 1.05 * incidence/100)) +
    NULL
  return(g)
}


#' The Contribution plot
#' 
#' ggcontribute plots the contribution of each variable to the model
#' 
#' @param x1 Either a logistic regression fitted using glm (base package) or lrm (rms package)  of the baseline model. 
#' @param x2 Either a logistic regression fitted using glm (base package) or lrm (rms package)  of the new (alternative) model.  
#' @param option_flag A flag to choose if the relative percentage of the Chi2-degrees of freedom are plotted.  
#' @return A ggplot object displaying the contribution of each variable to the model(s) using either Chi-square minus degrees of freedom or relative percentage contribution. If two models are provided, arrows show the change in contribution between models.
#' @import ggplot2
#' @import tidyr
#' @import dplyr
#' @export
#'
ggcontribute <- function(x1, x2 = NULL, option_flag = c("chi2","percent")){
  
  if (class(x1)[1] != "glm"  & class(x1)[1] != "lrm" ) { 
    stop("Oops - there must be a glm or lrm fit for the first variable")
  }
  if (!is.null(x2) & class(x2)[1] != "glm"  & class(x2)[1] != "lrm" ) {
    stop("Oops - there must be a glm or lrm fit for the second variable")
  } 
  if (class(x1)[1] == "glm") {
    an1 <- anova_glm(x1)
    data_type1 <- "glm"
  }
  if (class(x2)[1] == "glm") {
    an2 <- anova_glm(x2)
    data_type2 <- "glm"
  }
  if (class(x1)[1] == "lrm") {
    an1 <- anova(x1)
    data_type1  <- "lrm"
  }
  if (class(x2)[1] == "lrm") {
    an2 <- anova(x2)
    data_type1 <- "lrm"
  }
  if (length(option_flag) > 1){
    option_flag <- "chi2"
  }
  
  
  base_anova <- as.data.frame(an1[1:nrow(an1)-1,]) |> 
    mutate(variables= rownames(an1[1:nrow(an1)-1,])) |> 
    mutate(Xsqmindf = `Chi-Square` - `d.f.`) |> 
    arrange(-Xsqmindf) |> 
    filter(variables != "TOTAL")  
  
  vars = row.names(base_anova) 
  base_anova <- base_anova |> mutate(variables = vars)
  
  total.df <- sum(base_anova$Df)
  total.Xsqmindf <- sum(base_anova$Xsqmindf)
  base_anova <-  base_anova |> 
    mutate(Xsqmindf.percent = 100*Xsqmindf/total.Xsqmindf) |> 
    mutate(variables=factor(variables, levels=variables[order(Xsqmindf.percent)]))  |> 
    mutate(Model = "Baseline")
  
  if(!is.null(x2)){
    
    new_anova <- as.data.frame(an2[1:nrow(an2)-1,]) |> 
      mutate(variables= rownames(an2[1:nrow(an2)-1,])) |> 
      mutate(Xsqmindf = `Chi-Square` - `d.f.`) |> 
      arrange(-Xsqmindf) |> 
      filter(variables != "TOTAL")  
    
    vars = row.names(new_anova) 
    new_anova <- new_anova |> mutate(variables = vars)
    
    total.df <- sum(new_anova$Df)
    total.Xsqmindf <- sum(new_anova$Xsqmindf)
    new_anova <-  new_anova |> 
      mutate(Xsqmindf.percent = 100*Xsqmindf/total.Xsqmindf) |> 
      mutate(variables=factor(variables, levels=variables[order(Xsqmindf.percent)]))  |> 
      mutate(Model = "New")
    
    combined_anova = bind_rows(base_anova, new_anova)
    
    if(option_flag != "percent"){
      combined_anova_wide <- combined_anova |> 
        pivot_wider(id_cols = variables, names_from = Model, values_from = Xsqmindf) |> 
        filter(!is.na(Baseline))
    }
    
    
    if(option_flag == "percent"){
      combined_anova_wide <- combined_anova |> 
        pivot_wider(id_cols = variables, names_from = Model, values_from = Xsqmindf.percent) |> 
        filter(!is.na(Baseline))
    }
  }
  
  
  if(is.null(x2)){
    
    if(option_flag !="percent"){
      g <- ggplot()+
        geom_point(data = base_anova , aes(x=Xsqmindf, y=variables),size=4) +
        ylab(NULL) +
        xlab(expression(chi^{2}~{"-"}~df)) +
        theme(panel.background = element_rect(fill=NA),
              panel.grid.major.x = element_line(colour="grey70"),
              axis.ticks = element_blank(),
              axis.title.y = element_blank(),
              legend.title = element_blank(),
              legend.position = "bottom",
              text = element_text(size=12))
    }
    if(option_flag == "percent"){
      g <- ggplot() + 
        geom_point(data = base_anova , aes(x=Xsqmindf.percent, y=variables),size=4) +
        ylab(NULL) + 
        xlab("Relative variable contribution to the model (%)") + 
        scale_x_continuous(breaks=seq(0,100,5)) +
        theme(panel.background = element_rect(fill=NA),
              panel.grid.major.x = element_line(colour="grey70"),
              axis.ticks = element_blank(),
              axis.title.y = element_blank(),
              legend.title = element_blank(),
              legend.position = "bottom",
              text = element_text(size=12))
    }
  }
  
  
  if(!is.null(x2)){
    
    if(option_flag != "percent"){
      
      g <- ggplot()+
        geom_point(data = combined_anova , aes(x=Xsqmindf, y=variables, colour = Model),size=4) +
        geom_segment(data = combined_anova_wide, aes(x = Baseline, xend = New, y=variables, yend=variables), arrow = arrow(angle = 30, length = unit(0.15, "inches"),ends = "last", type = "open") ) +
        ylab(NULL) +
        xlab(expression(chi^{2}~{"-"}~df)) +
        theme(panel.background = element_rect(fill=NA),
              panel.grid.major.x = element_line(colour="grey70"),
              axis.ticks = element_blank(),
              axis.title.y = element_blank(),
              legend.title = element_blank(),
              legend.position = "bottom",
              text = element_text(size=12))
    } 
    
    if(option_flag == "percent"){
      g <- ggplot() + 
        geom_point(data = combined_anova , aes(x=Xsqmindf.percent, y=variables, colour = Model),size=4) +
        geom_segment(data = combined_anova_wide, aes(x = Baseline, xend = New, y=variables, yend=variables), arrow = arrow(angle = 30, length = unit(0.15, "inches"),ends = "last", type = "open") ) +
        ylab(NULL) +
        xlab("Relative variable contribution to the model (%)") + 
        scale_x_continuous(breaks=seq(0,100,5)) +
        theme(panel.background = element_rect(fill=NA),
              panel.grid.major.x = element_line(colour="grey70"),
              axis.ticks = element_blank(),
              axis.title.y = element_blank(),
              legend.title = element_blank(),
              legend.position = "bottom",
              text = element_text(size=12)) 
    }
  }
  
  return(g)
}


#' Statistical metrics
#' 
#' The function statistics.raplot calculates the reclassification metrics. Used by CI.raplot.
#' 
#' @param x1 Either a logistic regression fitted using glm (base package) or lrm (rms package) or calculated probabilities (eg through a logistic regression model) of the baseline model.  Must be between 0 & 1
#' @param x2 Either a logistic regression fitted using glm (base package) or lrm (rms package) or calculated probabilities (eg through a logistic regression model) of the new (alternative) model.   Must be between 0 & 1
#' @param y Binary of outcome of interest. Must be 0 or 1 (if fitted models are provided this is extracted from the fit which for an rms fit must have x = TRUE, y = TRUE). 
#' @param  t The risk threshold(s) for groups. eg t<-c(0,0.1,1) is a two group scenario with a threshold of 0.1 & t<-c(0,0.1,0.3,1)  is a three group scenario with thresholds at 0.1 and 0.3. Nb. If no t is provided it defaults to a single threshold at the prevalence of the cohort.  
#' @param  NRI_return Flag to return NRI metrics, default is FALSE.  
#' @return A matrix of metrics for use within CI.raplot
#' @import pROC
#' @import dplyr
#' @import tidyr
#' @export
#' 
statistics.raplot <- function(x1, x2, y,  t = NULL, NRI_return = FALSE) {   
  if (is.null(y))
    stop("Oops - there must be event data (y)")
  
  df <- data.frame(x1 = x1, x2 = x2, y = y) # y = 1 for the event, 0 for not the event.  
  
  #Remove rows with missing data: Unecessary - need to make sure no missing data in original input.
  #df <- df |>  
  #  filter(!is.na(x1)) |> 
  #  filter(!is.na(x2)) |> 
  #  filter(!is.na(y))
  
  n <- nrow(df)
  
  df <- df |>
    mutate(y == as.numeric(y)) # just in case
  
  u <- sort(unique(df$y))
  if (length(u) != 2 || u[1] != 0 || u[2] != 1)
    stop("Event (y) must have two values: 0 and 1")
  
  r <- range(df$x1, df$x2)
  if (r[1] < 0 || r[2] > 1)
    stop("x1 must be in [0,1]")
  
  n_event = sum(df$y)
  incidence <- 100*n_event/n
  
  if (is.null(t)) {t <- c(0, incidence/100,1) }
  
  # ROC
  roc.1 <- pROC::roc(df$y, df$x1)
  auc.x1 <- as.numeric(auc(roc.1))
  df_model.1 = data.frame(sens = roc.1$sensitivities, spec = roc.1$specificities, Prediction = roc.1$thresholds)
  df_model.1 <- df_model.1 |> mutate(Model = "Baseline")
  
  
  roc.2 <- pROC::roc(df$y, df$x2)
  auc.x2 <- as.numeric(auc(roc.2))
  auc.difference = auc.x2 - auc.x1
  df_model.2 = data.frame(sens = roc.2$sensitivities, spec = roc.2$specificities, Prediction = roc.2$thresholds)
  df_model.2 <- df_model.2 |> mutate(Model = "New")
  
  
  df <- df |> 
    mutate(risk.class.x1 = Hmisc::cut2(x1,t)) |> #risk groups based on thresholds
    mutate(mse_x1 = (y - x1)^2) |>  # mean squared errors of the baseline model
    mutate(risk.class.x2 = Hmisc::cut2(x2,t)) |> 
    mutate(difference = x2 - x1)  |>   # difference in probabilities
    mutate(mse_x2 = (y - x2)^2)  
  
  
  n_event = sum(df$y)
  n_non_event = nrow(df) - n_event
  n = n_event + n_non_event
  prevalence = n_event/n
  
  #NRI
  
  df_NRI <- df |> 
    group_by(y) |> 
    mutate(up = ifelse(as.numeric(risk.class.x2) > as.numeric(risk.class.x1),1,0))  |> 
    mutate(down = ifelse(as.numeric(risk.class.x2) < as.numeric(risk.class.x1),1,0)) |> 
    dplyr::summarise(
      n = n(),
      n_up = sum(up),
      n_down = sum(down),
      nri = (n_up - n_down)/n # for the non-event take the negative of this (below)
    ) |> 
    ungroup() |> 
    mutate(nri = ifelse(y == 0, -nri, nri))
  
  #IDI
  df_IDI <- df |> 
    group_by(y) |>  # only report for those with and without the disease seperately
    dplyr::summarise(IDI = mean(difference, na.rm = TRUE)) |> # the IDI is the mean change in Prediction
    ungroup() |> 
    mutate(IDI = ifelse(y == 0, -IDI, IDI)) # for the non-event convert this so that a positive mean change means the move in the right direction
  
  #Brier
  df_Brier <- df |> 
    dplyr::summarise(
      brier_baseline = mean(mse_x1),
      brier_new = mean(mse_x2),
      brier_skill = 100 * (brier_baseline - brier_new)/brier_baseline
    )
  df_model = bind_rows(df_model.1, df_model.2) 
  
  df_model <- df_model |> 
    mutate(Prediction = ifelse(sens == 1 & spec == 0, 0, ifelse(sens == 0 & spec == 1, 1, Prediction) ))
  
  df_model_Istats <- df_model |> 
    group_by(Model) |> 
    arrange(Prediction) |> 
    dplyr::summarise(
      IS = pracma::trapz(x = c(Prediction), y = c(sens)),  
      IP = pracma::trapz(x = c(Prediction), y = 1 - spec) )
  
  # Output
  if(NRI_return == TRUE  ){
    output <- list(n = n, 
                   n_event = n_event, 
                   n_non_event = n_non_event, 
                   Prevalence = prevalence,
                   NRI_up_event = as.integer(df_NRI[df_NRI$y == 1,]$n_up),
                   NRI_up_nonevent = as.integer(df_NRI[df_NRI$y == 0,]$n_up),
                   NRI_down_event = as.integer(df_NRI[df_NRI$y == 1,]$n_down),
                   NRI_down_nonevent = as.integer(df_NRI[df_NRI$y == 0,]$n_down),
                   NRI_event = df_NRI[df_NRI$y == 1,]$nri,
                   NRI_nonevent = df_NRI[df_NRI$y == 0,]$nri,
                   IDI_event = df_IDI[df_IDI$y == 1,]$IDI,
                   IDI_nonevent = df_IDI[df_IDI$y == 0,]$IDI,
                   IP_baseline = df_model_Istats[df_model_Istats$Model == "Baseline",]$IP,
                   IS_baseline = df_model_Istats[df_model_Istats$Model == "Baseline",]$IS,
                   IP_new = df_model_Istats[df_model_Istats$Model == "New",]$IP, 
                   IS_new = df_model_Istats[df_model_Istats$Model == "New",]$IS,
                   Brier_baseline = df_Brier$brier_baseline,
                   Brier_new = df_Brier$brier_new,
                   Brier_skill = df_Brier$brier_skill,
                   AUC_baseline = auc.x1,
                   AUC_new = auc.x2,
                   AUC_difference = auc.difference
    )
  }
  if(NRI_return == FALSE ){
    output <- list(n = n, 
                   n_event = n_event, 
                   n_non_event = n_non_event, 
                   Prevalence = prevalence,
                   IDI_event = df_IDI[df_IDI$y == 1,]$IDI,
                   IDI_nonevent = df_IDI[df_IDI$y == 0,]$IDI,
                   IP_baseline = df_model_Istats[df_model_Istats$Model == "Baseline",]$IP,
                   IS_baseline = df_model_Istats[df_model_Istats$Model == "Baseline",]$IS,
                   IP_new = df_model_Istats[df_model_Istats$Model == "New",]$IP, 
                   IS_new = df_model_Istats[df_model_Istats$Model == "New",]$IS,
                   Brier_baseline = df_Brier$brier_baseline,
                   Brier_new = df_Brier$brier_new,
                   Brier_skill = df_Brier$brier_skill,
                   AUC_baseline = auc.x1,
                   AUC_new = auc.x2,
                   AUC_difference = auc.difference
    )
  }
  
  return(output)
}


#' Extract confidence interval
#' 
#' Extract a confidence in interval from the bootstrapped results. Used by CI.raplot
#' 
#' @param results.boot The matrix of n.boot metrics from within CI.raplot
#' @param conf.level The confidence interval expressed between 0 & 1 (eg 95\%CI is conf.level = 0.95)
#' @param n.boot The number of bootstrapped samples
#' @param dp the number of decimal places to report the point estimate and confidence interval
#' @return A two column matrix with the metric name and statistic with a confidence interval
#' @import tidyr
#' @import dplyr
#' @export
extractCI <- function(results.boot, conf.level, n.boot , dp){
  
  n_vars = length(results.boot[[1]])
  results.df <- data.frame(matrix(nrow = n.boot, ncol = n_vars))   
  
  for (i in 1:n.boot) {
    temp_df <- as.data.frame(results.boot[[i]])
    results.df[i,] <- temp_df[1,]
  } 
  
  results.matrix_est <- results.df |> 
    summarise(across(where(is.numeric), ~round(quantile(.x,c(0.5)),dp)))
  results.matrix_lower_CI <- results.df |> 
    summarise(across(where(is.numeric), ~round(quantile(.x,c((1 - conf.level)/2)),dp)))
  results.matrix_upper_CI <- results.df |> 
    summarise(across(where(is.numeric), ~round(quantile(.x,c(1 - (1-conf.level)/2)),dp)))
  
  temp <- bind_rows(results.matrix_est, results.matrix_lower_CI , results.matrix_upper_CI)
  results.matrix <- t(temp)
  colnames(results.matrix) <- c("V1", "V2", "V3")
  
  results.matrix <- as_tibble(results.matrix)
  
  results.matrix$metric <- names(results.boot[[1]]) 
  
  results.matrix <- results.matrix |> 
    mutate(statistics = paste0( V1, " (CI: ",V2, " to ", V3,")" )) |> 
    select(metric, statistics)
  
  return(results.matrix)
}                

#' Extract NRI confidence intervals
#' 
#' Extract a confidence in interval from the bootstrapped results. Used by CI.NRI
#' 
#' @param results.boot The matrix of n.boot metrics from within CI.NRI
#' @param conf.level The confidence interval expressed between 0 & 1 (eg 95\%CI is conf.level = 0.95)
#' @param n.boot The number of bootstrapped samples
#' @param dp the number of decimal places to report the point estimate and confidence interval
#' @return A two column matrix with the metric name and statistic with a confidence interval
#' @import tidyr
#' @import dplyr
#' @export
extract_NRI_CI <- function(results.boot, conf.level, n.boot, dp){
  
  n_vars = sum(unlist(lapply( results.boot[[1]], class)) == "numeric" | unlist(lapply( results.boot[[1]], class)) == "integer" ) #because we don't calculate CIs for the confusion matrices
  results.df <- data.frame(matrix(nrow = n.boot, ncol = n_vars))   
  
  for (i in 1:n.boot) {
    temp <- results.boot[[i]]
    temp <- unlist(temp[1:n_vars])
    results.df[i,] <- temp 
  } 
  
  results.matrix_est <- results.df |> 
    summarise(across(where(is.numeric), ~round(quantile(.x,c(0.5)),dp)))
  results.matrix_lower_CI <- results.df |> 
    summarise(across(where(is.numeric), ~round(quantile(.x,c((1-conf.level)/2)),dp)))
  results.matrix_upper_CI <- results.df |> 
    summarise(across(where(is.numeric), ~round(quantile(.x,c(1 - (1-conf.level)/2)),dp)))
  
  temp <- bind_rows(results.matrix_est, results.matrix_lower_CI , results.matrix_upper_CI)
  results.matrix <- t(temp)
  colnames(results.matrix) <- c("V1", "V2", "V3")
  
  results.matrix <- as_tibble(results.matrix)
  results.matrix$metric <- names(results.boot[[1]])[1:10]
  
  results.matrix <- results.matrix |> 
    mutate(statistics = paste0( V1, " (CI: ",V2, " to ", V3,")" )) |> 
    select(metric, statistics)
  
  return(results.matrix)
} 


#' Statistical metrics with confidence intervals
#' 
#' The CI.raplot function produces summary metrics for risk assessment. Outputs the NRI, IDI, weighted NRI and category Free NRI all for those with events and those without events.  Also the AUCs of the two models and the comparison (DeLong) between AUCs. Output includes confidence intervals. Uses statistics.raplot.  Displayed graphically by raplot.
#' 
#' @param x1 Either a logistic regression fitted using glm (base package) or lrm (rms package) or calculated probabilities (eg through a logistic regression model) of the baseline model.  Must be between 0 & 1
#' @param x2 Either a logistic regression fitted using glm (base package) or lrm (rms package) or calculated probabilities (eg through a logistic regression model) of the new (alternative) model.   Must be between 0 & 1
#' @param y Binary of outcome of interest. Must be 0 or 1 (if fitted models are provided this is extracted from the fit which for an rms fit must have x = TRUE, y = TRUE). 
#' @param t The risk threshold(s) for groups. eg t<-c(0,0.1,1) is a two group model with a threshold of 0.1 & t<-c(0,0.1,0.3,1)  is a three group model with thresholds at 0.1 and 0.3.
#' @param NRI_return If NRI statistics are required (default = FALSE).
#' @param conf.level The confidence interval expressed as a fraction of 1 (ie 0.95 is the 95\% confidence interval )
#' @param n.boot  The number of "bootstraps" to use. Performance slows down with more bootstraps. For trialling result, use a low number (eg 5), for accuracy use a large number (eg 2000)
#' @param dp The number of decimal places to display
#' @return A list with the following elements:
#' \describe{
#'   \item{meta_data}{A data.frame with thresholds, confidence interval, number of bootstraps, input data type and decimal places.}
#'   \item{Metrics}{Point estimates of the statistical metrics (see function docs).}
#'   \item{Each_bootstrap_metrics}{List of per-bootstrap metric results.}
#'   \item{Summary_metrics}{A table of summary metrics with confidence intervals (e.g. Total, Events, Non-events, NRI, IDI, AUCs, Brier scores, etc.).}
#' }
#' @export 
#' @examples
#' # Quick example with subset of data and fewer bootstraps
#' data(data_risk)
#' data_subset <- data_risk[1:100, ]  # Use first 100 rows for speed
#' complete_cases <- complete.cases(data_subset)
#' data_clean <- data_subset[complete_cases, ]
#' y <- data_clean$outcome 
#' x1 <- data_clean$baseline
#' x2 <- data_clean$new
#' t <- c(0, 0.19, 1) 
#' output <- CI.raplot(x1, x2, y, t, conf.level = 0.95, n.boot = 10, dp = 2)
#' 
#' \donttest{
#' # Full dataset example with more bootstraps
#' data(data_risk)
#' complete_cases <- complete.cases(data_risk)
#' data_clean <- data_risk[complete_cases, ]
#' y <- data_clean$outcome 
#' x1 <- data_clean$baseline
#' x2 <- data_clean$new
#' t <- c(0, 0.19, 1) 
#' output <- CI.raplot(x1, x2, y, t, conf.level = 0.95, n.boot = 1000, dp = 2) 
#' }
#' @references  Pencina, M. J., D'Agostino, R. B., & Vasan, R. S. (2008). Evaluating the added stats::predictive ability of a new marker: From area under the ROC curve to reclassification and beyond. Statistics in Medicine, 27(2), 157-172. doi:10.1002/sim.2929
CI.raplot <- function(x1, x2 = NULL, y = NULL,  t = NULL, NRI_return = FALSE,  conf.level = 0.95, n.boot = 1000, dp = 3) {
  
  if (class(x1)[1] == "glm") {
    y = x1$y # must come first
    x1 = stats::predict(x1, type = "response")
    data_type = "glm"
  }
  if (!is.null(x2) & class(x2)[1] == "glm") {
    x2 = stats::predict(x2, type = "response")
  }
  if (class(x1)[1] == "lrm") {
    y = x1$y # must come first
    if (length(y) == 0 ) {
      stop("Fitted models with the rms package must be made using y = TRUE")
    }
    x1 = stats::predict(x1, type = "fitted")
    data_type = "lrm"
  }
  if (!is.null(x2) & class(x2)[1] == "lrm") {
    if (length(x2$y) == 0 ) {
      stop("Fitted models with the rms package must be made using y = TRUE")
    }
    x2 = stats::predict(x2, type = "fitted")
  }
  
  if (length(x2) ==1) {
    x2 = rep(x2,length(x1)) # If no x2 entered, then make the probability the value entered for all.
  }
  if (is.null(x2) ) {
    x2 = rep(0.5,length(x1)) # If no x2 entered, then make the probability 0.5 for all.
  }
  if (class(x1)[1] != "glm"  & class(x1)[1] != "lrm" ) {
    data_type = "User supplied"
  }
  if (length(x1) != length(x2) )
    stop("Reference (baseline) and New (Alt) model vectors must be the same length")
  if (is.null(y))
    stop("Oops - there must be event data (y)")
  
  ifelse(is.null(t),
         results <- statistics.raplot(x1, x2, y, t = NULL, NRI_return),  
         results <- statistics.raplot(x1, x2, y,  t, NRI_return)
  )
  
  # Get results for each bootstrapped sample
  results.boot <- list()
  
  for (i in 1:n.boot) {
    boot.index <- sample(length(y), replace = TRUE)
    risk.model1.boot <- x1[boot.index]
    cc.status.boot <- y[boot.index]    
    risk.model2.boot <- x2[boot.index]
    
    results.boot[[i]] <- statistics.raplot(x1 = risk.model1.boot, 
                                           x2 = risk.model2.boot, 
                                           y = cc.status.boot,
                                           t,
                                           NRI_return)
    
  }
  
  results.matrix <- extractCI(results.boot = results.boot, conf.level = conf.level, n.boot = n.boot, dp = dp)
  
  thresh = ifelse(is.null(t), "baseline",t)
  
  meta_data <- data.frame(Thresholds = thresh, `Confidence interval` = 100 * conf.level, `Number of bootstraps` = n.boot, `Input data type` = data_type, `# decimal places` = dp)
  return(list(meta_data = meta_data, Metrics = results, Each_bootstrap_metrics = results.boot, Summary_metrics = results.matrix))
}


#' List risk assessment metrics
#' 
#' Display the summary metrics
#' 
#' @param l List returned from CI.raplot
#' @return A tibble
summary.rap = function(l) {
  summary_metrics <- l[["Summary_metrics"]]
  return(summary_metrics)
}

#' List meta data
#' 
#' Display the meta data
#' 
#' @param l List returned from CI.raplot
#' @return A tibble
meta.rap = function(l) {
  meta_data <- l[["meta_data"]]
  return(meta_data)
}

#' Reclassification metrics with classes (ordinals) as inputs
#' 
#' The function statistics.classNRI calculates the NRI metrics for reclassification of data already in classes. For use by CI.classNRI.
#' 
#' @param c1 Risk class of Reference model (ordinal factor).
#' @param c2 Risk class of New model (ordinal factor)
#' @param y Binary of outcome of interest. Must be 0 or 1. 
#' @param s1 The savings or benefit when an event is reclassified to a higher group by the new model. i.e instead of counting as 1 an event classified to a higher group, it is counted as s1.
#' @param s2 The benefit when a non-event is reclassified to a lower group.  i.e instead of counting as 1 an event classified to a lower group, it is counted as s2.  
#' @return A matrix of metrics for use within CI.classNRI
#' @examples
#' # Quick example
#' data(data_class)
#' data_subset <- data_class[1:100, ]  # Use first 100 rows for speed
#' y <- data_subset$outcome 
#' c1 <- data_subset$base_class
#' c2 <- data_subset$new_class
#' output <- statistics.classNRI(c1, c2, y)
#' 
#' \donttest{
#' # Full dataset example
#' data(data_class)
#' y <- data_class$outcome 
#' c1 <- data_class$base_class
#' c2 <- data_class$new_class
#' output <- statistics.classNRI(c1, c2, y) 
#' }
#' @export
statistics.classNRI <- function(c1, c2, y,s1 = NULL, s2 = NULL) {    
  
  c_both <- c(c1, c2)
  if (!is.factor(c1)) {c1 <- factor(c1, levels = levels(factor(c_both))) }
  if (!is.factor(c2)) {c2 <- factor(c2, levels = levels(factor(c_both))) }
  
  df <- data.frame(c1 = c1, c2 = c2, event = y)
  
  #Remove rows with missing data
  df <- df |>  
    mutate(event = as.numeric(event)) |> 
    filter(!is.na(c1)) |> 
    filter(!is.na(c2)) |> 
    filter(!is.na(event))
  
  u <- sort(unique(df$event))
  if (length(u) != 2 || u[1] != 0 || u[2] != 1)
    stop("Outcome/Event must have two values: 0 and 1")
  
  
  n <- nrow(df)
  n_event = sum(df$event)
  incidence <- 100 * n_event/n
  
  #NRI
  df_NRI <- df |> 
    group_by(event) |> 
    mutate(up = ifelse(as.numeric(c2) > as.numeric(c1),1,0))  |> 
    mutate(down = ifelse(as.numeric(c2) < as.numeric(c1),1,0)) |> 
    summarise(
      n = n(),
      n_up = sum(up),
      n_down = sum(down),
      nri = (n_up - n_down)/n # for the non-event take the negative of this (below)
    ) |> 
    ungroup() |> 
    mutate(nri = ifelse(event == 0, -nri, nri))
  
  
  df_wNRI = NULL
  if (!is.null(s1) & !is.null(s2)) {
    df_wNRI <- df |> 
      group_by(event) |> 
      mutate(up = ifelse(as.numeric(c2) > as.numeric(c1),s1,0))  |> 
      mutate(down = ifelse(as.numeric(c2) < as.numeric(c1),s2,0)) |> 
      summarise(
        n = n(),
        n_up = sum(up),
        n_down = sum(down),
        nri = (n_up - n_down)/n # for the non-event take the negative of this (below)
      ) |> 
      ungroup() |> 
      mutate(nri = ifelse(event == 0, -nri, nri))
  }
  
  df_event <- df |> filter(event == 1)
  confusion.matrix_event <- table(df_event$c1, df_event$c2, dnn = c("Baseline", "New"))
  
  df_nonevent <- df |> filter(event == 0)
  confusion.matrix_nonevent <- table(df_nonevent$c1, df_nonevent$c2, dnn = c("Baseline", "New"))
  
  
  # Output
  output <- list(n = n, 
                 n_event = n_event, 
                 n_non_event = n - n_event, 
                 Prevalence = n_event/n,
                 NRI_up_event = as.integer(df_NRI[df_NRI$event == 1,]$n_up),
                 NRI_up_nonevent = as.integer(df_NRI[df_NRI$event == 0,]$n_up),
                 NRI_down_event = as.integer(df_NRI[df_NRI$event == 1,]$n_down),
                 NRI_down_nonevent = as.integer(df_NRI[df_NRI$event == 0,]$n_down),
                 NRI_event = df_NRI[df_NRI$event == 1,]$nri,
                 NRI_nonevent = df_NRI[df_NRI$event == 0,]$nri,
                 wNRI_event = df_wNRI[df_wNRI$event == 1,]$nri,
                 wNRI_nonevent = df_wNRI[df_wNRI$event == 0,]$nri,
                 confusion.matrix_event = confusion.matrix_event,
                 confusion.matrix_nonevent = confusion.matrix_nonevent
  )
  return(output)
}

#' Statistical metrics and confidence intervals for classes
#' 
#' The function CI.classNRI calculates the NRI statistics for reclassification of data already in classes with confidence intervals.  Uses statistics.classNRI.
#' 
#' @param c1 Risk classes of the baseline model (ordinal)
#' @param c2 Risk classes of new model
#' @param y Binary of outcome of interest. Must be 0 or 1. 
#' @param s1 The savings or benefit when am event is reclassified to a higher group by the new model (positive numeric)
#' @param s2 The benefit when a non-event is reclassified to a lower group (positive numeric)
#' @param conf.level The confidence interval expressed as a fraction of 1 (ie 0.95 is the 95\% confidence interval )
#' @param n.boot  The number of "bootstraps" to use. Performance slows down with more bootstraps. For trialling result, use a low number (eg 2), for accuracy use a large number (eg 2000)
#' @param dp The number of decimal places to display
#' @return A list with the following elements:
#' \describe{
#'   \item{meta_data}{Some overall meta data - Confidence Interval, number of bootstraps, s1, s2}
#'   \item{Metrics}{Point estimates of the statistical metrics.}
#'   \item{Each_bootstrap_metrics}{Point estimates of the statistical metrics for each bootstrapped sample.}
#'   \item{Summary_metrics}{Point estimates with confidence intervals of the statistical metrics (e.g. Total, Events, Non-events, Prevalence, NRI, IDI, confusion matrices).}
#' }
#' @return A matrix of metrics
#' @export
CI.classNRI <- function(c1, c2, y, s1 = NULL, s2 = NULL,  conf.level = 0.95, n.boot = 1000, dp = 3) {
  
  c_both <- c(c1, c2) 
  
  if (!is.factor(c1)) {c1 <- factor(c1, levels = levels(factor(c_both))) }
  if (!is.factor(c2)) {c2 <- factor(c2, levels = levels(factor(c_both))) }
  
  
  results <- statistics.classNRI(c1, c2,y,s1,s2)  
  
  # Get results for each bootstrapped sample
  results.boot <- list()
  
  for (i in 1:n.boot) {
    boot.index <- sample(length(y), replace = TRUE)
    risk.model1.boot <- c1[boot.index]
    risk.model2.boot <- c2[boot.index]
    cc.status.boot <- y[boot.index]           
    results.boot[[i]] <- statistics.classNRI(c1 = risk.model1.boot, 
                                             c2 = risk.model2.boot, 
                                             y = cc.status.boot,
                                             s1,
                                             s2)
  }
  
  results.matrix <- extract_NRI_CI(results.boot = results.boot, conf.level = conf.level, n.boot = n.boot, dp = dp) #Note, no CIs on the confusion matrices
  ifelse(is.null(s1),
         meta_data <- data.frame(`Confidence interval` = 100 * conf.level, `Number of bootstraps` = n.boot, `# decimal places` = dp),
         meta_data <- data.frame(`Confidence interval` = 100 * conf.level, `Number of bootstraps` = n.boot, `Weight s1` = s1, `Weight s2` = s2,`# decimal places` = dp))
  
  
  return(list(meta_data = meta_data, Metrics = results, Each_bootstrap_metrics = results.boot, Summary_metrics = results.matrix))
}


#' The function anova_glm() returns the Chi^2 and degrees of freedom for each variable & the same was anova.rms() does from lrm() in the rms package.
#' @param f A logistic regression fit created using glm (base package) 
#' @return A data frame with Chi-Square values and degrees of freedom for each variable in the model, plus a TOTAL row summarizing the overall model statistics.
#' @import tidyr 
#' @import dplyr
#' @export
anova_glm <- function(f){
  if (class(f)[1] != "glm") {
    stop("The fit must be a binary logistic model created with glm()")
  }
  df <- as.data.frame(matrix(nrow=length(f[["model"]]) ,ncol=2))
  rownames(df) = c(names(f[["model"]])[2:length(f[["model"]])],"TOTAL")
  colnames(df) <- c("Chi-Square","d.f.")
  total_deg_free = 0
  for (i in 2:length(f[["model"]])){
    deg_free <- length(levels(f[["model"]][[i]]))
    deg_free <- ifelse(deg_free == 0, 1, deg_free -1)
    var =stats::vcov(f)[(total_deg_free+2):(total_deg_free+deg_free+1),(total_deg_free+2):(total_deg_free+deg_free+1)]
    coef=f[["coefficients"]][(total_deg_free+2):(total_deg_free+deg_free+1)]
    df[i-1,1] <- coef %*% solve(var,coef)
    df[i-1,2] <- deg_free
    total_deg_free = total_deg_free + deg_free
    
  }
  
  total_deg_free <- sum(df[1:(nrow(df)-1),2])
  df[nrow(df),2] <- total_deg_free
  var = stats::vcov(f)[2:(total_deg_free+1),2:(total_deg_free+1)]
  coef=f[["coefficients"]][2:(total_deg_free+1)]
  df[nrow(df),1] <- coef %*% solve(var,coef)
  
  return(df)
}