R/diag-utils.R

In R4VN: Health Data Analysis and Publication-Ready Reporting

# Internal helpers for tabdiag() and tabdiagi()
# Not exported.

.r4vn_diag_resolve_data <- function(data = NULL) {
  if (!is.null(data)) {
    if (!is.data.frame(data)) stop("`data` must be a data frame.", call. = FALSE)
    return(data)
  }

  active <- NULL
  if (exists("usedf", mode = "function", inherits = TRUE)) {
    active <- tryCatch(usedf(), error = function(e) NULL)
  }

  if (is.null(active) || !is.data.frame(active)) {
    stop(
      "No data supplied and no active R4VN data set is available. ",
      "Supply `data = ...` or set active data with `usedf()`/`opendata(..., active = TRUE)`.",
      call. = FALSE
    )
  }

  active
}

.r4vn_diag_name <- function(expr, arg = "variable") {
  if (is.symbol(expr)) return(as.character(expr))
  if (is.character(expr) && length(expr) == 1L) return(expr)
  stop("`", arg, "` must be an unquoted variable name or a single character name.", call. = FALSE)
}

.r4vn_diag_names <- function(exprs) {
  if (!length(exprs)) return(character())
  vapply(
    exprs,
    function(z) {
      if (is.symbol(z)) return(as.character(z))
      if (is.character(z) && length(z) == 1L) return(z)
      stop("Diagnostic test variables in `...` must be unquoted names or character names.", call. = FALSE)
    },
    character(1)
  )
}

.r4vn_diag_safe_div <- function(a, b) {
  if (length(a) != 1L || length(b) != 1L || is.na(a) || is.na(b)) {
    return(NA_real_)
  }

  # A non-zero numerator divided by zero is an infinite ratio.
  # This is required, for example, for LR+ when specificity = 1
  # and LR- when specificity = 0.
  if (b == 0) {
    if (a == 0) return(NA_real_)
    return(sign(a) * Inf)
  }

  a / b
}

.r4vn_diag_positive_level <- function(x, positive = NULL, what = "variable") {
  z <- x[!is.na(x)]
  u <- unique(z)

  if (length(u) != 2L) {
    stop("`", what, "` must have exactly two non-missing levels.", call. = FALSE)
  }

  if (!is.null(positive)) {
    hit <- which(as.character(u) == as.character(positive))
    if (length(hit) != 1L) {
      stop(
        "`", positive, "` is not a unique observed level of `", what, "`.",
        call. = FALSE
      )
    }
    return(u[hit])
  }

  # Prefer common positive encodings.
  uc <- trimws(tolower(as.character(u)))
  positive_words <- c(
    "1", "true", "yes", "y", "positive", "pos", "case", "disease",
    "c\u00f3", "co", "d\u01b0\u01a1ng", "duong",
    "d\u01b0\u01a1ng t\u00ednh", "duong tinh"
  )
  hit <- match(positive_words, uc, nomatch = 0L)
  hit <- hit[hit > 0L]
  if (length(hit)) return(u[hit[1L]])

  # For factors, follow factor level order and use the second level.
  if (is.factor(x)) {
    lev <- levels(x)
    present <- lev[lev %in% as.character(u)]
    if (length(present) == 2L) {
      return(u[match(present[2L], as.character(u))])
    }
  }

  # For numeric/logical values use the larger/TRUE value.
  if (is.logical(z)) return(TRUE)
  if (is.numeric(z)) return(max(u))

  # Otherwise use the second observed level, but the printed result states it explicitly.
  u[2L]
}

.r4vn_diag_binary <- function(x, positive = NULL, what = "variable") {
  pos <- .r4vn_diag_positive_level(x, positive = positive, what = what)
  y <- rep(NA_integer_, length(x))
  ok <- !is.na(x)
  y[ok] <- as.integer(as.character(x[ok]) == as.character(pos))
  neg <- unique(x[ok & y == 0L])
  list(
    value = y,
    positive = as.character(pos),
    negative = if (length(neg)) as.character(neg[1L]) else NA_character_
  )
}

.r4vn_diag_is_binary <- function(x) {
  length(unique(x[!is.na(x)])) == 2L
}

.r4vn_diag_ci_requested <- function(ci, metric) {
  if (isTRUE(ci)) return(TRUE)
  if (isFALSE(ci) || is.null(ci)) return(FALSE)
  if (!is.character(ci)) stop("`ci` must be TRUE, FALSE, or a character vector of metric names.", call. = FALSE)

  aliases <- list(
    sens = c("sens", "sensitivity", "se", "recall", "tpr"),
    spec = c("spec", "specificity", "sp", "tnr"),
    ppv = c("ppv", "precision", "positive predictive value"),
    npv = c("npv", "negative predictive value"),
    accuracy = c("accuracy", "acc"),
    lr_pos = c("lr+", "lr_pos", "lrpos", "positive likelihood ratio"),
    lr_neg = c("lr-", "lr_neg", "lrneg", "negative likelihood ratio"),
    dor = c("dor", "diagnostic odds ratio"),
    auc = c("auc"),
    cutoff = c("cut", "cutoff", "threshold")
  )

  wanted <- tolower(trimws(ci))
  key <- aliases[[metric]]
  if (is.null(key)) return(FALSE)
  any(wanted %in% key)
}

.r4vn_diag_prop_ci <- function(x, n, level = 0.95, method = "wilson") {
  if (is.na(x) || is.na(n) || n <= 0) return(c(lower = NA_real_, upper = NA_real_))
  x <- as.numeric(x)
  n <- as.numeric(n)

  method <- match.arg(method, c("wilson", "exact"))

  if (method == "exact") {
    bt <- stats::binom.test(x, n, conf.level = level)
    return(c(lower = unname(bt$conf.int[1L]), upper = unname(bt$conf.int[2L])))
  }

  p <- x / n
  z <- stats::qnorm(1 - (1 - level) / 2)
  den <- 1 + z^2 / n
  center <- (p + z^2 / (2 * n)) / den
  half <- z * sqrt((p * (1 - p) / n) + z^2 / (4 * n^2)) / den
  c(lower = max(0, center - half), upper = min(1, center + half))
}

.r4vn_diag_log_ci <- function(est, se_log, level = 0.95) {
  if (!is.finite(est) || est <= 0 || !is.finite(se_log)) {
    return(c(lower = NA_real_, upper = NA_real_))
  }
  z <- stats::qnorm(1 - (1 - level) / 2)
  c(
    lower = exp(log(est) - z * se_log),
    upper = exp(log(est) + z * se_log)
  )
}

.r4vn_diag_2x2 <- function(
  tp, fp, fn, tn,
  prevalence = NULL,
  ci = FALSE,
  ci_level = 0.95,
  ci_method = c("auto", "wilson", "exact"),
  zero_correction = 0.5
) {
  ci_method <- match.arg(ci_method)
  prop_method <- if (ci_method == "exact") "exact" else "wilson"

  vals <- c(tp = tp, fp = fp, fn = fn, tn = tn)
  if (any(!is.finite(vals)) || any(vals < 0) || any(abs(vals - round(vals)) > .Machine$double.eps^0.5)) {
    stop("TP, FP, FN, and TN must be finite non-negative integer counts.", call. = FALSE)
  }

  tp <- as.numeric(tp); fp <- as.numeric(fp); fn <- as.numeric(fn); tn <- as.numeric(tn)
  n <- tp + fp + fn + tn
  if (n <= 0) stop("The 2 x 2 table contains no observations.", call. = FALSE)

  disease_n <- tp + fn
  nondisease_n <- fp + tn
  testpos_n <- tp + fp
  testneg_n <- fn + tn

  sens <- .r4vn_diag_safe_div(tp, disease_n)
  spec <- .r4vn_diag_safe_div(tn, nondisease_n)
  ppv_obs <- .r4vn_diag_safe_div(tp, testpos_n)
  npv_obs <- .r4vn_diag_safe_div(tn, testneg_n)
  accuracy <- .r4vn_diag_safe_div(tp + tn, n)
  bal_accuracy <- if (all(is.finite(c(sens, spec)))) (sens + spec) / 2 else NA_real_
  prevalence_obs <- .r4vn_diag_safe_div(disease_n, n)

  fpr <- if (is.finite(spec)) 1 - spec else NA_real_
  fnr <- if (is.finite(sens)) 1 - sens else NA_real_
  fdr <- if (is.finite(ppv_obs)) 1 - ppv_obs else NA_real_
  forate <- if (is.finite(npv_obs)) 1 - npv_obs else NA_real_

  lr_pos <- if (is.finite(sens) && is.finite(spec)) .r4vn_diag_safe_div(sens, 1 - spec) else NA_real_
  lr_neg <- if (is.finite(sens) && is.finite(spec)) .r4vn_diag_safe_div(1 - sens, spec) else NA_real_

  dor <- if (fp == 0 || fn == 0) {
    if (tp > 0 && tn > 0 && (fp == 0 || fn == 0)) Inf else NA_real_
  } else {
    (tp * tn) / (fp * fn)
  }

  youden <- if (all(is.finite(c(sens, spec)))) sens + spec - 1 else NA_real_
  f1 <- .r4vn_diag_safe_div(2 * tp, 2 * tp + fp + fn)

  mcc_den <- sqrt((tp + fp) * (tp + fn) * (tn + fp) * (tn + fn))
  mcc <- if (is.finite(mcc_den) && mcc_den > 0) (tp * tn - fp * fn) / mcc_den else NA_real_

  po <- accuracy
  pe <- (((tp + fp) * (tp + fn)) + ((fn + tn) * (fp + tn))) / (n^2)
  kappa <- if (is.finite(pe) && pe < 1) (po - pe) / (1 - pe) else NA_real_

  detection_rate <- tp / n
  detection_prevalence <- testpos_n / n

  ppv <- ppv_obs
  npv <- npv_obs
  prevalence_used <- prevalence_obs

  if (!is.null(prevalence)) {
    if (!is.numeric(prevalence) || length(prevalence) != 1L ||
        is.na(prevalence) || prevalence <= 0 || prevalence >= 1) {
      stop("`prevalence` must be a single number strictly between 0 and 1.", call. = FALSE)
    }
    prevalence_used <- prevalence
    if (all(is.finite(c(sens, spec)))) {
      ppv <- .r4vn_diag_safe_div(
        sens * prevalence,
        sens * prevalence + (1 - spec) * (1 - prevalence)
      )
      npv <- .r4vn_diag_safe_div(
        spec * (1 - prevalence),
        (1 - sens) * prevalence + spec * (1 - prevalence)
      )
    } else {
      ppv <- npv <- NA_real_
    }
  }

  estimates <- c(
    n = n,
    tp = tp, fp = fp, fn = fn, tn = tn,
    sensitivity = sens,
    specificity = spec,
    ppv = ppv,
    npv = npv,
    ppv_observed = ppv_obs,
    npv_observed = npv_obs,
    accuracy = accuracy,
    balanced_accuracy = bal_accuracy,
    prevalence = prevalence_used,
    prevalence_observed = prevalence_obs,
    fpr = fpr,
    fnr = fnr,
    fdr = fdr,
    `for` = forate,
    lr_pos = lr_pos,
    lr_neg = lr_neg,
    dor = dor,
    youden = youden,
    f1 = f1,
    mcc = mcc,
    kappa = kappa,
    detection_rate = detection_rate,
    detection_prevalence = detection_prevalence
  )

  ci_out <- list()

  if (.r4vn_diag_ci_requested(ci, "sens")) {
    ci_out$sensitivity <- .r4vn_diag_prop_ci(tp, disease_n, ci_level, prop_method)
  }
  if (.r4vn_diag_ci_requested(ci, "spec")) {
    ci_out$specificity <- .r4vn_diag_prop_ci(tn, nondisease_n, ci_level, prop_method)
  }

  # PPV/NPV binomial CIs apply to the observed predictive values.
  if (.r4vn_diag_ci_requested(ci, "ppv")) {
    ci_out$ppv <- if (is.null(prevalence)) {
      .r4vn_diag_prop_ci(tp, testpos_n, ci_level, prop_method)
    } else c(lower = NA_real_, upper = NA_real_)
  }
  if (.r4vn_diag_ci_requested(ci, "npv")) {
    ci_out$npv <- if (is.null(prevalence)) {
      .r4vn_diag_prop_ci(tn, testneg_n, ci_level, prop_method)
    } else c(lower = NA_real_, upper = NA_real_)
  }
  if (.r4vn_diag_ci_requested(ci, "accuracy")) {
    ci_out$accuracy <- .r4vn_diag_prop_ci(tp + tn, n, ci_level, prop_method)
  }

  if (.r4vn_diag_ci_requested(ci, "lr_pos")) {
    se <- if (tp > 0 && fp > 0 && disease_n > 0 && nondisease_n > 0) {
      sqrt(1 / tp - 1 / disease_n + 1 / fp - 1 / nondisease_n)
    } else NA_real_
    ci_out$lr_pos <- .r4vn_diag_log_ci(lr_pos, se, ci_level)
  }

  if (.r4vn_diag_ci_requested(ci, "lr_neg")) {
    se <- if (fn > 0 && tn > 0 && disease_n > 0 && nondisease_n > 0) {
      sqrt(1 / fn - 1 / disease_n + 1 / tn - 1 / nondisease_n)
    } else NA_real_
    ci_out$lr_neg <- .r4vn_diag_log_ci(lr_neg, se, ci_level)
  }

  if (.r4vn_diag_ci_requested(ci, "dor")) {
    cc <- c(tp, fp, fn, tn)
    corrected <- FALSE
    if (any(cc == 0)) {
      cc <- cc + zero_correction
      corrected <- TRUE
    }
    tp2 <- cc[1L]; fp2 <- cc[2L]; fn2 <- cc[3L]; tn2 <- cc[4L]
    est2 <- (tp2 * tn2) / (fp2 * fn2)
    se <- sqrt(1 / tp2 + 1 / fp2 + 1 / fn2 + 1 / tn2)
    ci_out$dor <- .r4vn_diag_log_ci(est2, se, ci_level)
    attr(ci_out$dor, "zero_corrected") <- corrected
  }

  list(
    estimates = estimates,
    ci = ci_out,
    prevalence_adjusted = !is.null(prevalence),
    zero_correction = zero_correction
  )
}

.r4vn_diag_add_ci_columns <- function(df, res) {
  supported <- c("sensitivity", "specificity", "ppv", "npv", "accuracy", "lr_pos", "lr_neg", "dor")
  for (nm in supported) {
    lo <- hi <- NA_real_
    if (!is.null(res$ci[[nm]])) {
      lo <- unname(res$ci[[nm]][1L])
      hi <- unname(res$ci[[nm]][2L])
    }
    df[[paste0(nm, "_low")]] <- lo
    df[[paste0(nm, "_high")]] <- hi
  }
  df
}

.r4vn_diag_threshold_row <- function(
  marker, criterion, cutoff, direction, y, x,
  prevalence, ci, ci_level, ci_method, zero_correction
) {
  ok <- !is.na(y) & !is.na(x)
  yy <- y[ok]
  xx <- x[ok]

  pred <- if (direction == "<") {
    as.integer(xx >= cutoff)
  } else {
    as.integer(xx <= cutoff)
  }

  tp <- sum(pred == 1L & yy == 1L)
  fp <- sum(pred == 1L & yy == 0L)
  fn <- sum(pred == 0L & yy == 1L)
  tn <- sum(pred == 0L & yy == 0L)

  m <- .r4vn_diag_2x2(
    tp, fp, fn, tn,
    prevalence = prevalence,
    ci = ci,
    ci_level = ci_level,
    ci_method = ci_method,
    zero_correction = zero_correction
  )

  e <- m$estimates
  out <- data.frame(
    marker = marker,
    criterion = criterion,
    cutoff = as.numeric(cutoff),
    direction = if (direction == "<") ">=" else "<=",
    n = unname(e["n"]),
    tp = unname(e["tp"]),
    fp = unname(e["fp"]),
    fn = unname(e["fn"]),
    tn = unname(e["tn"]),
    sensitivity = unname(e["sensitivity"]),
    specificity = unname(e["specificity"]),
    ppv = unname(e["ppv"]),
    npv = unname(e["npv"]),
    accuracy = unname(e["accuracy"]),
    balanced_accuracy = unname(e["balanced_accuracy"]),
    lr_pos = unname(e["lr_pos"]),
    lr_neg = unname(e["lr_neg"]),
    dor = unname(e["dor"]),
    youden = unname(e["youden"]),
    f1 = unname(e["f1"]),
    mcc = unname(e["mcc"]),
    kappa = unname(e["kappa"]),
    fpr = unname(e["fpr"]),
    fnr = unname(e["fnr"]),
    fdr = unname(e["fdr"]),
    `for` = unname(e["for"]),
    prevalence = unname(e["prevalence"]),
    detection_rate = unname(e["detection_rate"]),
    detection_prevalence = unname(e["detection_prevalence"]),
    stringsAsFactors = FALSE
  )
  .r4vn_diag_add_ci_columns(out, m)
}

.r4vn_diag_binary_row <- function(
  marker, y, test, test_positive = NULL,
  prevalence, ci, ci_level, ci_method, zero_correction
) {
  tb <- .r4vn_diag_binary(test, positive = test_positive, what = marker)
  ok <- !is.na(y) & !is.na(tb$value)
  yy <- y[ok]
  pred <- tb$value[ok]

  tp <- sum(pred == 1L & yy == 1L)
  fp <- sum(pred == 1L & yy == 0L)
  fn <- sum(pred == 0L & yy == 1L)
  tn <- sum(pred == 0L & yy == 0L)

  m <- .r4vn_diag_2x2(
    tp, fp, fn, tn,
    prevalence = prevalence,
    ci = ci,
    ci_level = ci_level,
    ci_method = ci_method,
    zero_correction = zero_correction
  )

  e <- m$estimates
  out <- data.frame(
    marker = marker,
    criterion = paste0("Positive: ", tb$positive),
    cutoff = NA_real_,
    direction = "=",
    n = unname(e["n"]),
    tp = unname(e["tp"]),
    fp = unname(e["fp"]),
    fn = unname(e["fn"]),
    tn = unname(e["tn"]),
    sensitivity = unname(e["sensitivity"]),
    specificity = unname(e["specificity"]),
    ppv = unname(e["ppv"]),
    npv = unname(e["npv"]),
    accuracy = unname(e["accuracy"]),
    balanced_accuracy = unname(e["balanced_accuracy"]),
    lr_pos = unname(e["lr_pos"]),
    lr_neg = unname(e["lr_neg"]),
    dor = unname(e["dor"]),
    youden = unname(e["youden"]),
    f1 = unname(e["f1"]),
    mcc = unname(e["mcc"]),
    kappa = unname(e["kappa"]),
    fpr = unname(e["fpr"]),
    fnr = unname(e["fnr"]),
    fdr = unname(e["fdr"]),
    `for` = unname(e["for"]),
    prevalence = unname(e["prevalence"]),
    detection_rate = unname(e["detection_rate"]),
    detection_prevalence = unname(e["detection_prevalence"]),
    stringsAsFactors = FALSE
  )
  .r4vn_diag_add_ci_columns(out, m)
}

.r4vn_diag_test_positive_for <- function(positive, marker, index) {
  if (is.null(positive)) return(NULL)
  if (length(positive) == 1L && is.null(names(positive))) return(positive[[1L]])
  if (!is.null(names(positive)) && marker %in% names(positive)) return(positive[[marker]])
  if (is.null(names(positive)) && length(positive) >= index) return(positive[[index]])
  NULL
}

.r4vn_diag_prepare_predictor <- function(x, marker) {
  if (is.ordered(x)) return(as.numeric(x))
  if (is.numeric(x)) return(as.numeric(x))
  if (inherits(x, "Date") || inherits(x, "POSIXt")) return(as.numeric(x))
  if (is.factor(x) && nlevels(x) > 2L) {
    stop(
      "`", marker, "` is an unordered factor with more than two levels. ",
      "Use an ordered factor or a numeric diagnostic score.",
      call. = FALSE
    )
  }
  suppressWarnings(z <- as.numeric(as.character(x)))
  if (all(is.na(z) == is.na(x))) return(z)
  stop(
    "`", marker, "` must be binary, numeric, date/time, or an ordered factor for ROC analysis.",
    call. = FALSE
  )
}

.r4vn_diag_auc_score <- function(y, score) {
  ok <- !is.na(y) & !is.na(score)
  y <- y[ok]
  score <- score[ok]
  m <- sum(y == 1L)
  n <- sum(y == 0L)
  if (m < 1L || n < 1L) return(NA_real_)
  r <- rank(score, ties.method = "average")
  (sum(r[y == 1L]) - m * (m + 1) / 2) / (m * n)
}

.r4vn_diag_auc_placements <- function(y, score) {
  ok <- !is.na(y) & !is.na(score)
  y <- y[ok]
  score <- score[ok]
  cases <- score[y == 1L]
  controls <- score[y == 0L]
  m <- length(cases)
  n <- length(controls)
  if (!m || !n) {
    return(list(auc = NA_real_, variance = NA_real_, case = numeric(), control = numeric()))
  }

  case_place <- vapply(
    cases,
    function(v) mean((v > controls) + 0.5 * (v == controls)),
    numeric(1)
  )
  control_place <- vapply(
    controls,
    function(v) mean((cases > v) + 0.5 * (cases == v)),
    numeric(1)
  )
  auc <- mean(case_place)
  v_case <- if (m > 1L) stats::var(case_place) / m else NA_real_
  v_control <- if (n > 1L) stats::var(control_place) / n else NA_real_
  variance <- if (is.finite(v_case) && is.finite(v_control)) v_case + v_control else NA_real_

  list(auc = auc, variance = variance, case = case_place, control = control_place)
}

.r4vn_diag_roc_coords <- function(y, x, direction) {
  u <- sort(unique(x[is.finite(x)]))
  if (length(u) < 2L) return(data.frame())
  mids <- u[-length(u)] + (u[-1L] - u[-length(u)]) / 2
  thresholds <- c(-Inf, mids, Inf)

  rows <- lapply(thresholds, function(cut) {
    pred <- if (direction == "<") x >= cut else x <= cut
    tp <- sum(pred & y == 1L)
    fn <- sum(!pred & y == 1L)
    fp <- sum(pred & y == 0L)
    tn <- sum(!pred & y == 0L)
    data.frame(
      threshold = cut,
      sensitivity = .r4vn_diag_safe_div(tp, tp + fn),
      specificity = .r4vn_diag_safe_div(tn, tn + fp),
      stringsAsFactors = FALSE
    )
  })
  do.call(rbind, rows)
}

.r4vn_diag_roc <- function(y, x, direction = c("auto", "<", ">")) {
  direction <- match.arg(direction)
  ok <- !is.na(y) & !is.na(x) & is.finite(x)
  y <- as.integer(y[ok])
  x <- as.numeric(x[ok])
  if (length(unique(y)) != 2L) stop("ROC requires both outcome levels.", call. = FALSE)
  if (length(unique(x)) < 2L) stop("ROC requires a non-constant marker.", call. = FALSE)

  if (direction == "auto") {
    med_case <- stats::median(x[y == 1L])
    med_control <- stats::median(x[y == 0L])
    direction <- if (med_case >= med_control) "<" else ">"
  }
  score <- if (direction == "<") x else -x
  placements <- .r4vn_diag_auc_placements(y, score)
  coords <- .r4vn_diag_roc_coords(y, x, direction)

  structure(
    list(
      response = y,
      predictor = x,
      score = score,
      direction = direction,
      auc = placements$auc,
      variance = placements$variance,
      case_placements = placements$case,
      control_placements = placements$control,
      coordinates = coords,
      cases = sum(y == 1L),
      controls = sum(y == 0L)
    ),
    class = "r4vn_diag_roc"
  )
}

.r4vn_diag_auc_ci <- function(roc, level = 0.95) {
  auc <- roc$auc
  v <- roc$variance
  if (!is.finite(auc) || !is.finite(v) || v < 0) {
    return(c(lower = NA_real_, upper = NA_real_))
  }
  z <- stats::qnorm(1 - (1 - level) / 2)
  se <- sqrt(v)
  c(lower = max(0, auc - z * se), upper = min(1, auc + z * se))
}

.r4vn_diag_curve_area <- function(x, y, from, to) {
  ok <- is.finite(x) & is.finite(y)
  x <- x[ok]
  y <- y[ok]
  if (length(x) < 2L || !is.finite(from) || !is.finite(to)) return(NA_real_)
  lim <- sort(c(from, to))
  from <- lim[1L]
  to <- lim[2L]
  if (to <= from) return(0)

  ord <- order(x, y)
  x <- x[ord]
  y <- y[ord]
  ux <- unique(x)
  uy <- vapply(ux, function(v) max(y[x == v], na.rm = TRUE), numeric(1))
  if (length(ux) < 2L) return(NA_real_)

  boundary <- stats::approx(ux, uy, xout = c(from, to), rule = 2, ties = "ordered")$y
  inside <- ux > from & ux < to
  xx <- c(from, ux[inside], to)
  yy <- c(boundary[1L], uy[inside], boundary[2L])
  sum(diff(xx) * (yy[-length(yy)] + yy[-1L]) / 2)
}

.r4vn_diag_partial_auc <- function(roc, range, focus = c("specificity", "sensitivity"), correct = FALSE) {
  focus <- match.arg(focus)
  if (is.null(range) || length(range) != 2L) return(NA_real_)
  coords <- roc$coordinates
  if (!nrow(coords)) return(NA_real_)

  if (focus == "specificity") {
    x <- 1 - coords$specificity
    y <- coords$sensitivity
    xr <- 1 - range
    lo <- min(xr)
    hi <- max(xr)
    random_area <- (hi^2 - lo^2) / 2
  } else {
    x <- coords$sensitivity
    y <- coords$specificity
    lo <- min(range)
    hi <- max(range)
    random_area <- (hi - lo) - (hi^2 - lo^2) / 2
  }
  raw <- .r4vn_diag_curve_area(x, y, lo, hi)
  if (!isTRUE(correct) || !is.finite(raw)) return(raw)

  max_area <- hi - lo
  den <- max_area - random_area
  if (!is.finite(den) || den <= 0) return(NA_real_)
  corrected <- 0.5 * (1 + (raw - random_area) / den)
  max(0, min(1, corrected))
}

.r4vn_diag_partial_auc_ci <- function(
  roc, range, focus = "specificity", correct = FALSE,
  level = 0.95, boot = 2000
) {
  y <- roc$response
  x <- roc$predictor
  case_idx <- which(y == 1L)
  control_idx <- which(y == 0L)
  if (!length(case_idx) || !length(control_idx)) return(c(lower = NA_real_, upper = NA_real_))

  vals <- rep(NA_real_, boot)
  for (b in seq_len(boot)) {
    idx <- c(
      sample(case_idx, length(case_idx), replace = TRUE),
      sample(control_idx, length(control_idx), replace = TRUE)
    )
    rr <- tryCatch(.r4vn_diag_roc(y[idx], x[idx], direction = roc$direction), error = function(e) NULL)
    if (!is.null(rr)) vals[b] <- .r4vn_diag_partial_auc(rr, range, focus, correct)
  }
  vals <- vals[is.finite(vals)]
  if (length(vals) < max(20L, floor(boot * 0.10))) return(c(lower = NA_real_, upper = NA_real_))
  alpha <- (1 - level) / 2
  q <- stats::quantile(vals, probs = c(alpha, 1 - alpha), na.rm = TRUE, names = FALSE, type = 7)
  c(lower = q[1L], upper = q[2L])
}

.r4vn_diag_best_thresholds <- function(roc, best, target = 0.95) {
  if (is.null(best) || !length(best)) return(data.frame())
  best <- unique(tolower(best))
  allc <- roc$coordinates
  allc <- allc[is.finite(allc$threshold), , drop = FALSE]
  if (!nrow(allc)) return(data.frame())

  out <- list()
  k <- 0L
  add <- function(criteria, thresholds) {
    thresholds <- unique(as.numeric(thresholds[is.finite(thresholds)]))
    if (!length(thresholds)) return(invisible(NULL))
    k <<- k + 1L
    out[[k]] <<- data.frame(
      criterion = rep(criteria, length(thresholds)),
      cutoff = thresholds,
      stringsAsFactors = FALSE
    )
    invisible(NULL)
  }

  tol <- sqrt(.Machine$double.eps)
  for (b in best) {
    if (b %in% c("youden", "y")) {
      score <- allc$sensitivity + allc$specificity - 1
      mx <- max(score, na.rm = TRUE)
      add("Youden", allc$threshold[abs(score - mx) <= tol])
    } else if (b %in% c("closest", "closest.topleft", "topleft")) {
      distance <- (1 - allc$sensitivity)^2 + (1 - allc$specificity)^2
      mn <- min(distance, na.rm = TRUE)
      add("Closest top-left", allc$threshold[abs(distance - mn) <= tol])
    } else if (b %in% c("ruleout", "rule-out", "sensitivity", "sens")) {
      cand <- allc[allc$sensitivity >= target, , drop = FALSE]
      if (nrow(cand)) {
        mx <- max(cand$specificity, na.rm = TRUE)
        cand <- cand[abs(cand$specificity - mx) <= tol, , drop = FALSE]
        mx2 <- max(cand$sensitivity, na.rm = TRUE)
        cand <- cand[abs(cand$sensitivity - mx2) <= tol, , drop = FALSE]
        add(paste0("Rule-out: sensitivity >= ", round(target * 100, 1), "%"), cand$threshold)
      }
    } else if (b %in% c("rulein", "rule-in", "specificity", "spec")) {
      cand <- allc[allc$specificity >= target, , drop = FALSE]
      if (nrow(cand)) {
        mx <- max(cand$sensitivity, na.rm = TRUE)
        cand <- cand[abs(cand$sensitivity - mx) <= tol, , drop = FALSE]
        mx2 <- max(cand$specificity, na.rm = TRUE)
        cand <- cand[abs(cand$specificity - mx2) <= tol, , drop = FALSE]
        add(paste0("Rule-in: specificity >= ", round(target * 100, 1), "%"), cand$threshold)
      }
    } else {
      stop(
        "Unknown `best` method: ", b,
        ". Use 'youden', 'closest', 'ruleout', or 'rulein'.",
        call. = FALSE
      )
    }
  }

  if (!length(out)) return(data.frame())
  unique(do.call(rbind, out))
}

.r4vn_diag_cut_ci <- function(roc, method, level = 0.95, boot = 2000) {
  method <- tolower(method)
  if (!method %in% c("youden", "closest", "closest.topleft", "topleft")) {
    return(c(lower = NA_real_, upper = NA_real_))
  }
  method <- if (method == "youden") "youden" else "closest"
  y <- roc$response
  x <- roc$predictor
  case_idx <- which(y == 1L)
  control_idx <- which(y == 0L)
  vals <- rep(NA_real_, boot)

  for (b in seq_len(boot)) {
    idx <- c(
      sample(case_idx, length(case_idx), replace = TRUE),
      sample(control_idx, length(control_idx), replace = TRUE)
    )
    rr <- tryCatch(.r4vn_diag_roc(y[idx], x[idx], direction = roc$direction), error = function(e) NULL)
    if (is.null(rr)) next
    bt <- .r4vn_diag_best_thresholds(rr, best = method, target = 0.95)
    if (nrow(bt)) vals[b] <- stats::median(bt$cutoff, na.rm = TRUE)
  }
  vals <- vals[is.finite(vals)]
  if (length(vals) < max(20L, floor(boot * 0.10))) return(c(lower = NA_real_, upper = NA_real_))
  alpha <- (1 - level) / 2
  q <- stats::quantile(vals, probs = c(alpha, 1 - alpha), na.rm = TRUE, names = FALSE, type = 7)
  c(lower = q[1L], upper = q[2L])
}

.r4vn_diag_compare_delong <- function(y, score1, score2, level = 0.95) {
  p1 <- .r4vn_diag_auc_placements(y, score1)
  p2 <- .r4vn_diag_auc_placements(y, score2)
  m <- length(p1$case)
  n <- length(p1$control)
  auc1 <- p1$auc
  auc2 <- p2$auc
  diff <- auc1 - auc2

  cov_case <- if (m > 1L) stats::cov(p1$case, p2$case) / m else NA_real_
  cov_control <- if (n > 1L) stats::cov(p1$control, p2$control) / n else NA_real_
  cov12 <- if (is.finite(cov_case) && is.finite(cov_control)) cov_case + cov_control else NA_real_
  var_diff <- p1$variance + p2$variance - 2 * cov12
  if (is.finite(var_diff) && var_diff < 0 && abs(var_diff) < 1e-12) var_diff <- 0
  se <- if (is.finite(var_diff) && var_diff >= 0) sqrt(var_diff) else NA_real_
  z <- if (is.finite(se) && se > 0) diff / se else NA_real_
  p <- if (is.finite(z)) 2 * stats::pnorm(-abs(z)) else if (is.finite(se) && se == 0 && diff == 0) 1 else NA_real_
  zz <- stats::qnorm(1 - (1 - level) / 2)
  lower <- if (is.finite(se)) diff - zz * se else NA_real_
  upper <- if (is.finite(se)) diff + zz * se else NA_real_

  list(auc1 = auc1, auc2 = auc2, difference = diff, lower = lower, upper = upper, p = p, se = se)
}

.r4vn_diag_compare_bootstrap <- function(y, score1, score2, level = 0.95, boot = 2000) {
  auc1 <- .r4vn_diag_auc_score(y, score1)
  auc2 <- .r4vn_diag_auc_score(y, score2)
  diff <- auc1 - auc2
  case_idx <- which(y == 1L)
  control_idx <- which(y == 0L)
  vals <- rep(NA_real_, boot)

  for (b in seq_len(boot)) {
    idx <- c(
      sample(case_idx, length(case_idx), replace = TRUE),
      sample(control_idx, length(control_idx), replace = TRUE)
    )
    vals[b] <- .r4vn_diag_auc_score(y[idx], score1[idx]) -
      .r4vn_diag_auc_score(y[idx], score2[idx])
  }
  vals <- vals[is.finite(vals)]
  if (length(vals) < max(20L, floor(boot * 0.10))) {
    return(list(auc1 = auc1, auc2 = auc2, difference = diff,
                lower = NA_real_, upper = NA_real_, p = NA_real_, se = NA_real_))
  }
  alpha <- (1 - level) / 2
  q <- stats::quantile(vals, probs = c(alpha, 1 - alpha), na.rm = TRUE, names = FALSE, type = 7)
  se <- stats::sd(vals)
  z <- if (is.finite(se) && se > 0) diff / se else NA_real_
  p <- if (is.finite(z)) 2 * stats::pnorm(-abs(z)) else NA_real_
  list(auc1 = auc1, auc2 = auc2, difference = diff,
       lower = q[1L], upper = q[2L], p = p, se = se)
}

.r4vn_diag_format_num <- function(x, digit = 2) {
  ifelse(
    is.na(x), "",
    ifelse(is.infinite(x), ifelse(x > 0, "Inf", "-Inf"), formatC(x, format = "f", digits = digit))
  )
}

.r4vn_diag_format_p <- function(p, digit = 3) {
  ifelse(
    is.na(p), "",
    ifelse(
      p < 10^(-digit),
      paste0("<", formatC(10^(-digit), format = "f", digits = digit)),
      formatC(p, format = "f", digits = digit)
    )
  )
}

.r4vn_diag_format_est_ci <- function(est, lo, hi, digit = 2, percent = FALSE) {
  mult <- if (percent) 100 else 1
  est2 <- est * mult
  lo2 <- lo * mult
  hi2 <- hi * mult
  base <- .r4vn_diag_format_num(est2, digit)

  has_ci <- !is.na(lo2) & !is.na(hi2)
  base[has_ci] <- paste0(
    .r4vn_diag_format_num(est2[has_ci], digit),
    " (",
    .r4vn_diag_format_num(lo2[has_ci], digit),
    "\u2013",
    .r4vn_diag_format_num(hi2[has_ci], digit),
    ")"
  )
  base
}

.r4vn_diag_metric_alias <- function(x) {
  x <- tolower(trimws(x))
  map <- c(
    sens = "sensitivity", sensitivity = "sensitivity", se = "sensitivity",
    spec = "specificity", specificity = "specificity", sp = "specificity",
    ppv = "ppv", npv = "npv",
    acc = "accuracy", accuracy = "accuracy",
    balanced = "balanced_accuracy", balanced_accuracy = "balanced_accuracy",
    `lr+` = "lr_pos", lr_pos = "lr_pos", lrpos = "lr_pos",
    `lr-` = "lr_neg", lr_neg = "lr_neg", lrneg = "lr_neg",
    dor = "dor", youden = "youden", f1 = "f1", mcc = "mcc", kappa = "kappa",
    tp = "tp", fp = "fp", fn = "fn", tn = "tn",
    fpr = "fpr", fnr = "fnr", fdr = "fdr", `for` = "for",
    prevalence = "prevalence",
    detection_rate = "detection_rate",
    detection_prevalence = "detection_prevalence",
    auc = "auc", cutoff = "cutoff", cut = "cutoff", threshold = "cutoff",
    n = "n"
  )
  out <- unname(map[x])
  out[is.na(out)] <- x[is.na(out)]
  out
}

.r4vn_diag_threshold_display <- function(df, ci, show = "all", hide = NULL, digit = 2) {
  if (!nrow(df)) return(df)

  # Cell counts are intentionally not presented as ordinary performance
  # measures.  They remain fully available in the returned object and are
  # presented as a conventional 2 x 2 classification table in
  # `$tables$Confusion_matrix` / the Viewer.
  keep_metrics <- c(
    "n",
    "sensitivity", "specificity", "ppv", "npv",
    "lr_pos", "lr_neg", "dor",
    "accuracy", "balanced_accuracy", "youden",
    "f1", "mcc", "kappa",
    "fpr", "fnr", "fdr", "for",
    "prevalence", "detection_rate", "detection_prevalence"
  )

  if (length(show) == 1L && is.character(show)) {
    profile <- tolower(trimws(show))
    if (profile %in% c("core", "default", "publication")) {
      keep_metrics <- c(
        "n", "sensitivity", "specificity", "ppv", "npv",
        "lr_pos", "lr_neg", "accuracy"
      )
      if (identical(profile, "publication")) keep_metrics <- c(keep_metrics, "dor")
    } else if (profile %in% c("minimal", "brief")) {
      keep_metrics <- c("n", "sensitivity", "specificity", "ppv", "npv")
    } else if (profile %in% c("counts", "cells")) {
      # `counts` requests the 2 x 2 classification table rather than a
      # vertical TP/FP/FN/TN metric list.  N is retained only as a compact
      # descriptor in the performance object; the Viewer suppresses this
      # otherwise redundant table for the counts-only profile.
      keep_metrics <- "n"
    } else if (!identical(profile, "all")) {
      keep_metrics <- intersect(keep_metrics, .r4vn_diag_metric_alias(show))
    }
  } else if (!(length(show) == 1L && identical(tolower(show), "all"))) {
    keep_metrics <- intersect(keep_metrics, .r4vn_diag_metric_alias(show))
  }
  if (!is.null(hide)) {
    keep_metrics <- setdiff(keep_metrics, .r4vn_diag_metric_alias(hide))
  }

  out <- data.frame(
    Marker = if ("marker_label" %in% names(df)) df$marker_label else df$marker,
    Criterion = df$criterion,
    Cutoff = ifelse(
      is.na(df$cutoff), "",
      paste0(df$direction, .r4vn_diag_format_num(df$cutoff, digit))
    ),
    stringsAsFactors = FALSE,
    check.names = FALSE
  )

  labels <- c(
    n = "N", tp = "TP", fp = "FP", fn = "FN", tn = "TN",
    sensitivity = "Sensitivity", specificity = "Specificity",
    ppv = "PPV", npv = "NPV",
    lr_pos = "LR+", lr_neg = "LR-", dor = "DOR",
    accuracy = "Accuracy", balanced_accuracy = "Balanced accuracy",
    youden = "Youden", f1 = "F1", mcc = "MCC", kappa = "Kappa",
    fpr = "FPR", fnr = "FNR", fdr = "FDR", `for` = "FOR",
    prevalence = "Prevalence",
    detection_rate = "Detection rate",
    detection_prevalence = "Detection prevalence"
  )

  percent_metrics <- c(
    "sensitivity", "specificity", "ppv", "npv",
    "accuracy", "balanced_accuracy",
    "fpr", "fnr", "fdr", "for",
    "prevalence", "detection_rate", "detection_prevalence"
  )

  ci_names <- c("sensitivity", "specificity", "ppv", "npv", "accuracy", "lr_pos", "lr_neg", "dor")

  # Some threshold rows/analysis modes do not carry every possible diagnostic
  # metric. Keep the display layer defensive: only render metrics that actually
  # exist in the result object and have a display label.
  keep_metrics <- keep_metrics[keep_metrics %in% names(df) & keep_metrics %in% names(labels)]

  for (nm in keep_metrics) {
    if (nm %in% c("n", "tp", "fp", "fn", "tn")) {
      out[[labels[[nm]]]] <- .r4vn_diag_format_num(df[[nm]], 0)
      next
    }

    ci_cols <- c(paste0(nm, "_low"), paste0(nm, "_high"))
    if (nm %in% ci_names && all(ci_cols %in% names(df))) {
      lo <- df[[ci_cols[1L]]]
      hi <- df[[ci_cols[2L]]]
      out[[labels[[nm]]]] <- .r4vn_diag_format_est_ci(
        df[[nm]], lo, hi, digit = digit, percent = nm %in% percent_metrics
      )
    } else {
      val <- df[[nm]]
      if (nm %in% percent_metrics) val <- val * 100
      out[[labels[[nm]]]] <- .r4vn_diag_format_num(val, digit)
    }
  }

  out
}

.r4vn_diag_summary_display <- function(df, ci, digit = 2, p_digit = 3) {
  if (!nrow(df)) return(df)
  auc_text <- .r4vn_diag_format_est_ci(
    df$auc, df$auc_low, df$auc_high,
    digit = digit, percent = FALSE
  )
  out <- data.frame(
    Marker = if ("marker_label" %in% names(df)) df$marker_label else df$marker,
    N = df$n,
    Cases = df$cases,
    Controls = df$controls,
    AUC = auc_text,
    `AUC SE` = .r4vn_diag_format_num(df$auc_se, digit + 1L),
    `AUC p-value` = .r4vn_diag_format_p(df$auc_p, p_digit),
    Direction = df$direction,
    stringsAsFactors = FALSE,
    check.names = FALSE
  )

  if ("partial_auc" %in% names(df) && any(!is.na(df$partial_auc))) {
    out$`Partial AUC` <- .r4vn_diag_format_est_ci(
      df$partial_auc,
      df$partial_auc_low,
      df$partial_auc_high,
      digit = digit,
      percent = FALSE
    )
  }
  out
}

.r4vn_diag_print_df <- function(x) {
  if (!is.data.frame(x) || !nrow(x)) return(invisible(NULL))
  print(x, row.names = FALSE, right = FALSE)
  invisible(NULL)
}