## ----include = FALSE----------------------------------------------------------
#library(knitr)
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----include=FALSE------------------------------------------------------------
library(TwoCutoff)


## ----echo=FALSE, fig.width=14, fig.height=5-----------------------------------
# library(TwoCutoff)
if (requireNamespace("DiagrammeR", quietly = TRUE)) {
DiagrammeR::grViz("
digraph TwoCutoffPipeline {

  graph [layout = dot, rankdir = TB]

  node [
    shape = box,
    style = filled,
    fillcolor = lightblue,
    fontname = Helvetica,
    fontsize = 12,   # smaller font
    width = 2.2,     # narrower boxes
    height = 0.8     # shorter boxes
  ]

  adjust_score [label = 'adjust_score()\nMixture + XGBoost']

  cutoff [label = 'Cutoff Derivation\nPercentile-based or ROC-based']

  plotcut [label = 'plot_two_cutoff()\nVisualization + Confusion Matrix']

  eval [label = 'evaluate_performance()\nSensitivity, Specificity,\nPPV, NPV']

  compare [label = 'compare_performance()\nSide-by-side summary']

  dca [label = 'dca_analysis()\nDecision Curve Analysis']

  adjust_score -> cutoff
  cutoff -> eval
  eval -> compare
  cutoff -> plotcut
  compare -> dca
}
" , width = 800, height = 600)
}else{
    cat("Diagram omitted because the 'DiagrammeR' package is not installed.")

}

## -----------------------------------------------------------------------------

 n <- 1000

 # Set seed for reproducibility
 set.seed(123)

 # Outcome
 DISEASE <- rbinom(n, 1, 0.35)

 # Biomarker depends on disease
 PTAU <- rnorm(n, mean = 25 + 10 * DISEASE, sd = 8)

 # Confounders (optional)
 AGE <- round(rnorm(n, 35, 7))
 SEX <- sample(c(0, 1), n, replace = TRUE)
 BMI <- rnorm(n, 25, 4)
 HTN <- rbinom(n, 1, 0.3)
 DM <- rbinom(n, 1, 0.2)

 df <- data.frame(PTAU, AGE, SEX, BMI, HTN, DM, DISEASE)

## -----------------------------------------------------------------------------
res <- adjust_score(
  df,
  biomarker = "PTAU", outcome = "DISEASE",
  confounders = c("AGE","SEX","BMI","HTN","DM"))
head(res$data[,c("PTAU","p_disease","adjusted_risk")])

## -----------------------------------------------------------------------------
cutoffs_p <- derive_cutoffs_percentile(
  res,
  biomarker = "PTAU"
)

cutoffs_p$c_L
cutoffs_p$c_H

table(cutoffs_p$class)

## -----------------------------------------------------------------------------
cutoffs_s <- derive_cutoffs_sensspec(
  res,
  biomarker = "PTAU",
  outcome = "DISEASE",
  sens_target = 0.80,
  spec_target = 0.80
)

cutoffs_s$c_L
cutoffs_s$c_H

table(cutoffs_s$class)

## -----------------------------------------------------------------------------
perf_p <- evaluate_performance(
  cutoffs_p,
  outcome = "DISEASE"
)

## ----echo=F-------------------------------------------------------------------
knitr::kable(perf_p, align = "c")

## -----------------------------------------------------------------------------
perf_s <- evaluate_performance(
  cutoffs_s,
  outcome = "DISEASE"
)

## ----echo=F-------------------------------------------------------------------
knitr::kable(perf_s, align = "c")

## -----------------------------------------------------------------------------
comparison <- compare_performance(
  cutoffs_p,
  cutoffs_s,
  outcome = "DISEASE"
)

## ----echo=FALSE---------------------------------------------------------------
knitr::kable(comparison, align = "c")

## ----fig.width=8, fig.height=6------------------------------------------------
plot_two_cutoff(
  cutoffs_p,
  biomarker = "PTAU",
  outcome = "DISEASE"
)

## ----fig.width=8, fig.height=6------------------------------------------------
# Create combined plotting object
cut_plot <- cutoffs_p

cut_plot$c_L <- c(cutoffs_p$c_L, cutoffs_s$c_L)
cut_plot$c_H <- c(cutoffs_p$c_H, cutoffs_s$c_H)

plot_two_cutoff(
  cut_plot,
  biomarker = "PTAU",
  outcome = "DISEASE",
  panel_title = c("Percentile-based", "ROC-based")
)

## ----fig.width=8, fig.height=6------------------------------------------------
dca_out <- dca_analysis(
  res,
  outcome = "DISEASE"
)

head(dca_out$data)

print(dca_out$plot)