R/tabdiagi.R

In R4VN: Health Data Analysis and Publication-Ready Reporting

#' Immediate diagnostic accuracy from a 2 x 2 table
#'
#' @description
#' `tabdiagi()` calculates diagnostic test accuracy directly from the four
#' cell counts of a 2 x 2 table. It is intended especially for teaching,
#' checking hand calculations, and analyses where only aggregated counts are
#' available.
#'
#' @param tp True positives. With positional input this is the first number.
#' @param fp False positives. With positional input this is the second number.
#' @param fn False negatives. With positional input this is the third number.
#' @param tn True negatives. With positional input this is the fourth number.
#' @param a,b,c,d Optional aliases for `tp`, `fp`, `fn`, and `tn`,
#'   respectively. Use either the TP/FP/FN/TN names or the a/b/c/d names,
#'   not both.
#' @param prevalence Optional population disease prevalence, strictly between
#'   0 and 1. When supplied, PPV and NPV are recalculated from sensitivity,
#'   specificity, and this prevalence. The observed predictive values remain
#'   available in the returned object.
#' @param ci `FALSE` (default), `TRUE`, or a character vector specifying which
#'   confidence intervals to calculate. Examples include
#'   `ci = c("sens", "spec")`,
#'   `ci = c("ppv", "npv", "lr+", "lr-", "dor")`, or `ci = TRUE`.
#' @param ci_level Confidence level, default 0.95.
#' @param ci_method Method for binomial proportion confidence intervals:
#'   `"auto"` or `"wilson"` uses Wilson intervals; `"exact"` uses exact
#'   binomial intervals. LR+/LR- and DOR use log-scale intervals.
#' @param zero_correction Continuity correction used only when a zero cell
#'   prevents a finite log-scale DOR confidence interval. Default 0.5.
#' @param digit Number of digits displayed for estimates.
#' @param show Logical; open the formatted result in the Viewer. Default `TRUE`.
#' @param console Logical; also print the teaching-oriented result in the Console. Default `FALSE`.
#'
#' @details
#' The cell layout is:
#'
#' ```
#'                         Reference standard
#'                       Positive     Negative
#' Test positive            TP           FP
#' Test negative            FN           TN
#' ```
#'
#' The function calculates the full set of commonly used 2 x 2 diagnostic
#' measures: sensitivity, specificity, PPV, NPV, accuracy, balanced accuracy,
#' LR+, LR-, diagnostic odds ratio (DOR), Youden index, F1 score, Matthews
#' correlation coefficient (MCC), Cohen's kappa, false-positive rate (FPR),
#' false-negative rate (FNR), false-discovery rate (FDR), false-omission rate
#' (FOR), disease prevalence, detection rate, and detection prevalence.
#'
#' Confidence intervals are optional because they can make console output much
#' wider. `ci = TRUE` requests all currently supported intervals.
#' Sensitivity, specificity, PPV, NPV, and accuracy use binomial intervals.
#' LR+ and LR- use conventional log-scale intervals. DOR uses a log-scale
#' interval; if any cell is zero, the requested `zero_correction` is applied
#' for the DOR interval calculation and this is recorded in the returned
#' object.
#'
#' When `prevalence` is supplied, PPV and NPV in the main results are
#' prevalence-adjusted. Because simple binomial confidence intervals no longer
#' apply to these adjusted predictive values, their CI entries are returned as
#' missing. The observed PPV and NPV remain in `result$estimates` as
#' `ppv_observed` and `npv_observed`.
#'
#' @return
#' Invisibly returns an object of class `r4vn_diagi` and `r4vn_stat` containing
#' `table`, `estimates`, `ci`, `settings`, and `notes`.
#'
#' @section Typical teaching use:
#' Use four counts directly:
#'
#' ```
#' tabdiagi(80, 20, 10, 90)
#' ```
#'
#' or use explicit names:
#'
#' ```
#' tabdiagi(tp = 80, fp = 20, fn = 10, tn = 90)
#' ```
#'
#' @section Confidence intervals:
#' Request only sensitivity and specificity confidence intervals:
#'
#' ```
#' tabdiagi(80, 20, 10, 90, ci = c("sens", "spec"))
#' ```
#'
#' Request all supported confidence intervals:
#'
#' ```
#' tabdiagi(80, 20, 10, 90, ci = TRUE)
#' ```
#'
#' @section Predictive values at a target prevalence:
#' To show how PPV and NPV change when disease prevalence is 10 percent:
#'
#' ```
#' tabdiagi(80, 20, 10, 90, prevalence = 0.10)
#' ```
#'
#' @examples
#' tabdiagi(80, 20, 10, 90)
#' tabdiagi(tp = 80, fp = 20, fn = 10, tn = 90)
#' tabdiagi(80, 20, 10, 90, ci = c("sens", "spec"))
#' tabdiagi(80, 20, 10, 90, ci = TRUE, show = FALSE)
#' tabdiagi(80, 20, 10, 90, prevalence = 0.10)
#'
#' @export
tabdiagi <- function(
  tp = NULL, fp = NULL, fn = NULL, tn = NULL,
  a = NULL, b = NULL, c = NULL, d = NULL,
  prevalence = NULL,
  ci = FALSE,
  ci_level = 0.95,
  ci_method = c("auto", "wilson", "exact"),
  zero_correction = 0.5,
  digit = 2,
  show = TRUE,
  console = FALSE
) {
  ci_method <- match.arg(ci_method)

  aliases_used <- !all(vapply(list(a, b, c, d), is.null, logical(1)))
  primary_used <- !all(vapply(list(tp, fp, fn, tn), is.null, logical(1)))

  if (aliases_used && primary_used) {
    stop("Use either `tp/fp/fn/tn` or `a/b/c/d`, not both.", call. = FALSE)
  }

  if (aliases_used) {
    tp <- a; fp <- b; fn <- c; tn <- d
  }

  if (any(vapply(list(tp, fp, fn, tn), is.null, logical(1)))) {
    stop("Four counts are required: TP, FP, FN, and TN.", call. = FALSE)
  }

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

  tab <- matrix(
    c(tp, fn, tp + fn,
      fp, tn, fp + tn,
      tp + fp, fn + tn, tp + fp + fn + tn),
    nrow = 3L,
    byrow = FALSE,
    dimnames = list(
      c("Test positive", "Test negative", "Total"),
      c("Disease positive", "Disease negative", "Total")
    )
  )

  notes <- character()
  if (!is.null(prevalence)) {
    notes <- c(
      notes,
      paste0(
        "PPV and NPV were adjusted to a disease prevalence of ",
        formatC(prevalence * 100, format = "f", digits = 1), "%."
      )
    )
  }
  if (!is.null(res$ci$dor) && isTRUE(attr(res$ci$dor, "zero_corrected"))) {
    notes <- c(
      notes,
      paste0(
        "A continuity correction of ", zero_correction,
        " was used for the DOR confidence interval because at least one cell was zero."
      )
    )
  }

  out <- list(
    table = tab,
    estimates = res$estimates,
    ci = res$ci,
    settings = list(
      prevalence = prevalence,
      ci = ci,
      ci_level = ci_level,
      ci_method = ci_method,
      zero_correction = zero_correction,
      digit = digit
    ),
    notes = notes,
    call = match.call()
  )
  class(out) <- c("r4vn_diagi", "r4vn_stat")

  .r4vn_show(out, show = show, console = console)
}

#' @export
print.r4vn_diagi <- function(x, ...) {
  digit <- x$settings$digit
  ci <- x$settings$ci
  e <- x$estimates

  cat("\nDiagnostic accuracy from a 2 x 2 table\n\n")
  print(x$table, quote = FALSE)

  metrics <- c(
    sensitivity = "Sensitivity",
    specificity = "Specificity",
    ppv = "PPV",
    npv = "NPV",
    accuracy = "Accuracy",
    balanced_accuracy = "Balanced accuracy",
    lr_pos = "LR+",
    lr_neg = "LR-",
    dor = "Diagnostic odds ratio",
    youden = "Youden index",
    f1 = "F1 score",
    mcc = "Matthews correlation coefficient",
    kappa = "Cohen's kappa",
    fpr = "False-positive rate",
    fnr = "False-negative rate",
    fdr = "False-discovery rate",
    `for` = "False-omission rate",
    prevalence = "Disease prevalence",
    detection_rate = "Detection rate",
    detection_prevalence = "Detection prevalence"
  )

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

  display <- data.frame(
    Measure = unname(metrics),
    Result = character(length(metrics)),
    stringsAsFactors = FALSE,
    check.names = FALSE
  )

  for (i in seq_along(metrics)) {
    nm <- names(metrics)[i]
    lo <- hi <- NA_real_
    if (!is.null(x$ci[[nm]])) {
      lo <- unname(x$ci[[nm]][1L])
      hi <- unname(x$ci[[nm]][2L])
    }

    display$Result[i] <- if (nm %in% c("sensitivity", "specificity", "ppv", "npv", "accuracy")) {
      .r4vn_diag_format_est_ci(e[[nm]], lo, hi, digit, percent = TRUE)
    } else if (nm %in% c("lr_pos", "lr_neg", "dor")) {
      .r4vn_diag_format_est_ci(e[[nm]], lo, hi, digit, percent = FALSE)
    } else {
      val <- e[[nm]]
      if (nm %in% percent) val <- val * 100
      .r4vn_diag_format_num(val, digit)
    }
  }

  cat("\nDiagnostic measures\n")
  .r4vn_diag_print_df(display)

  if (length(x$notes)) {
    cat("\nNotes:\n")
    for (z in x$notes) cat("- ", z, "\n", sep = "")
  }

  invisible(x)
}


# ============================================================================
# Viewer renderer for tabdiagi()
# ============================================================================
.r4vn_diagi_display_table <- function(x) {
  digit <- x$settings$digit
  e <- x$estimates
  metrics <- c(
    sensitivity = "Sensitivity", specificity = "Specificity", ppv = "PPV", npv = "NPV",
    accuracy = "Accuracy", balanced_accuracy = "Balanced accuracy", lr_pos = "LR+", lr_neg = "LR-",
    dor = "Diagnostic odds ratio", youden = "Youden index", f1 = "F1 score",
    mcc = "Matthews correlation coefficient", kappa = "Cohen's kappa",
    fpr = "False-positive rate", fnr = "False-negative rate", fdr = "False-discovery rate",
    `for` = "False-omission rate", prevalence = "Disease prevalence",
    detection_rate = "Detection rate", detection_prevalence = "Detection prevalence"
  )
  percent <- c("sensitivity", "specificity", "ppv", "npv", "accuracy", "balanced_accuracy",
               "fpr", "fnr", "fdr", "for", "prevalence", "detection_rate", "detection_prevalence")
  display <- data.frame(Measure = unname(metrics), Result = character(length(metrics)), stringsAsFactors = FALSE, check.names = FALSE)
  for (i in seq_along(metrics)) {
    nm <- names(metrics)[i]; lo <- hi <- NA_real_
    if (!is.null(x$ci[[nm]])) { lo <- unname(x$ci[[nm]][1L]); hi <- unname(x$ci[[nm]][2L]) }
    display$Result[i] <- if (nm %in% c("sensitivity", "specificity", "ppv", "npv", "accuracy")) {
      .r4vn_diag_format_est_ci(e[[nm]], lo, hi, digit, percent = TRUE)
    } else if (nm %in% c("lr_pos", "lr_neg", "dor")) {
      .r4vn_diag_format_est_ci(e[[nm]], lo, hi, digit, percent = FALSE)
    } else {
      val <- e[[nm]]; if (nm %in% percent) val <- val * 100; .r4vn_diag_format_num(val, digit)
    }
  }
  display
}

.r4vn_viewer_tabdiagi <- function(x) {
  body <- paste0(
    .r4vn_view_section("2 x 2 table", x$table),
    .r4vn_view_section("Diagnostic measures", .r4vn_diagi_display_table(x))
  )
  subtitle <- if (is.null(x$settings$prevalence)) "Diagnostic accuracy from entered counts" else
    paste0("Diagnostic accuracy \u00b7 population prevalence = ", formatC(100 * x$settings$prevalence, format = "f", digits = 1), "%")
  .r4vn_view_document("Diagnostic accuracy from a 2 x 2 table", body, notes = x$notes,
                      subtitle = subtitle, prefix = "r4vn-tabdiagi-")
}