R/statistics-utils.R

In R4VN: Health Data Analysis and Publication-Ready Reporting

# R4VN internal statistical utilities
#
# This file contains non-exported helpers shared by the immediate commands,
# data-based statistical commands, epidemiological commands, and models.
# Public functions are defined in immediate.R, statistics.R, epi.R, and models.R.

# ============================================================================
# Consolidated from: 000-utils.R
# ============================================================================
# Internal R4VN helpers. Not exported.

.r4vn_stat_data <- function(data = NULL) {
  if (!is.null(data)) {
    if (!is.data.frame(data)) stop("`data` must be a data frame.", call. = FALSE)
    return(data)
  }
  if (exists(".r4vn_resolve_analysis_data", mode = "function", inherits = TRUE)) {
    out <- tryCatch(get(".r4vn_resolve_analysis_data", mode = "function", inherits = TRUE)(NULL), error = function(e) NULL)
    if (is.data.frame(out)) return(out)
  }
  if (exists(".r4vn_get_active", mode = "function", inherits = TRUE)) {
    out <- tryCatch(get(".r4vn_get_active", mode = "function", inherits = TRUE)(), error = function(e) NULL)
    if (is.data.frame(out)) return(out)
  }
  stop("No data supplied and no active data is available. Use `data = ...` or `usedf()` first.", call. = FALSE)
}

.r4vn_eval_var <- function(expr, data, env, arg = "variable") {
  if (is.null(expr) || identical(expr, quote(NULL))) return(NULL)
  if (is.character(expr) && length(expr) == 1L && expr %in% names(data)) return(data[[expr]])
  out <- tryCatch(eval(expr, envir = data, enclos = env), error = function(e) NULL)
  if (is.null(out)) stop(sprintf("Could not evaluate `%s` in `data`.", arg), call. = FALSE)
  if (length(out) != nrow(data)) stop(sprintf("`%s` must have one value per row of `data`.", arg), call. = FALSE)
  out
}

.r4vn_name <- function(expr, fallback = "Variable") {
  if (is.symbol(expr)) return(as.character(expr))
  paste(deparse(expr), collapse = "") %||% fallback
}

`%||%` <- function(x, y) if (is.null(x) || !length(x)) y else x

.r4vn_num <- function(x, digits = 3L) {
  ifelse(is.na(x), "", ifelse(is.infinite(x), ifelse(x > 0, "Inf", "-Inf"), formatC(x, format = "f", digits = digits)))
}

.r4vn_p <- function(x, digits = 3L) {
  lim <- 10^(-digits)
  ifelse(is.na(x), "", ifelse(x < lim, paste0("<", formatC(lim, format = "f", digits = digits)), formatC(x, format = "f", digits = digits)))
}

.r4vn_ci <- function(est, low, high, digits = 3L) {
  paste0(.r4vn_num(est, digits), " (", .r4vn_num(low, digits), ", ", .r4vn_num(high, digits), ")")
}

.r4vn_result <- function(title, sections = list(), notes = NULL, raw = list(), call = NULL) {
  structure(list(title = title, sections = sections, notes = notes, raw = raw, call = call), class = "r4vn_stat")
}

#' Print an R4VN statistical result
#' @param x An object of class `r4vn_stat`.
#' @param ... Unused.
#' @return The input object, invisibly.
#' @export
print.r4vn_stat <- function(x, ...) {
  cat("\n", x$title, "\n", paste(rep("-", nchar(x$title)), collapse = ""), "\n", sep = "")
  for (nm in names(x$sections)) {
    obj <- x$sections[[nm]]
    if (!is.null(nm) && nzchar(nm)) cat("\n", nm, "\n", sep = "")
    if (is.matrix(obj)) print(obj, quote = FALSE, right = TRUE)
    else print(obj, row.names = FALSE, right = TRUE)
  }
  if (length(x$notes)) cat("\n", paste0("Note: ", x$notes, collapse = "\n"), "\n", sep = "")
  invisible(x)
}

.r4vn_view_escape <- function(x) {
  x <- as.character(x)
  x[is.na(x)] <- ""
  x <- gsub("&", "&amp;", x, fixed = TRUE)
  x <- gsub("<", "&lt;", x, fixed = TRUE)
  x <- gsub(">", "&gt;", x, fixed = TRUE)
  x <- gsub('"', "&quot;", x, fixed = TRUE)
  x
}

.r4vn_view_p_is_sig <- function(x) {
  z <- trimws(as.character(x))
  out <- rep(FALSE, length(z))
  less <- grepl("^<", z)
  if (any(less)) {
    v <- suppressWarnings(as.numeric(sub("^<\\s*", "", z[less])))
    out[less] <- is.finite(v) & v <= 0.05
  }
  if (any(!less)) {
    v <- suppressWarnings(as.numeric(z[!less]))
    out[!less] <- is.finite(v) & v < 0.05
  }
  out
}

.r4vn_view_table <- function(x, row_names = NULL, compact = FALSE) {
  if (is.null(x)) return("")
  if (is.table(x)) x <- as.matrix(x)
  if (is.matrix(x)) {
    rn <- rownames(x)
    x <- as.data.frame(x, stringsAsFactors = FALSE, check.names = FALSE)
    if (is.null(row_names)) row_names <- !is.null(rn)
    if (isTRUE(row_names) && !is.null(rn)) x <- data.frame(` ` = rn, x, stringsAsFactors = FALSE, check.names = FALSE)
  } else if (!is.data.frame(x)) {
    x <- as.data.frame(x, stringsAsFactors = FALSE, check.names = FALSE)
  }
  if (!nrow(x) && !ncol(x)) return('<div class="r4vn-empty">No results.</div>')

  cn <- names(x)

  # Keep table headers aligned with the cells underneath:
  # the first (stub/label) column is left-aligned and all result
  # columns are right-aligned.
  head <- vapply(
    seq_along(cn),
    function(j) {
      cls <- if (j == 1L) "stub" else "value"
      paste0(
        '<th class="',
        cls,
        '">',
        .r4vn_view_escape(cn[j]),
        "</th>"
      )
    },
    character(1)
  )
  head <- paste0(head, collapse = "")

  pcols <- grepl("(^p$|(^|[._ -])p($|[._ -])|p[._ -]?value|^pr\\(|prob)", tolower(cn))
  rows <- character(nrow(x))
  for (i in seq_len(nrow(x))) {
    cells <- character(ncol(x))
    for (j in seq_len(ncol(x))) {
      val0 <- x[[j]][i]
      val <- if (is.na(val0)) "" else as.character(val0)
      esc <- .r4vn_view_escape(val)
      if (pcols[j] && .r4vn_view_p_is_sig(val)) esc <- paste0("<strong>", esc, "</strong>")
      cls <- if (j == 1L) "stub" else "value"
      cells[j] <- paste0('<td class="', cls, '">', esc, "</td>")
    }
    rows[i] <- paste0("<tr>", paste0(cells, collapse = ""), "</tr>")
  }
  paste0('<div class="r4vn-table-scroll"><table class="r4vn-result-table', if (compact) ' compact' else '', '"><thead><tr>', head,
         '</tr></thead><tbody>', paste0(rows, collapse = ""), '</tbody></table></div>')
}

.r4vn_view_section <- function(title, object, class = NULL, description = NULL) {
  ttl <- if (is.null(title) || !nzchar(title)) "" else paste0('<h2>', .r4vn_view_escape(title), '</h2>')
  desc <- if (is.null(description) || !nzchar(description)) "" else paste0('<div class="section-note">', .r4vn_view_escape(description), '</div>')
  paste0('<section class="r4vn-section ', .r4vn_view_escape(class %||% ""), '">', ttl, desc, .r4vn_view_table(object), '</section>')
}

.r4vn_view_key_values <- function(x) {
  if (is.null(x) || !is.data.frame(x) || ncol(x) < 2L) return(.r4vn_view_table(x))
  a <- as.character(x[[1L]]); b <- as.character(x[[2L]])
  items <- paste0('<div class="kv-item"><div class="kv-key">', .r4vn_view_escape(a), '</div><div class="kv-value">', .r4vn_view_escape(b), '</div></div>')
  paste0('<div class="kv-grid">', paste0(items, collapse = ""), '</div>')
}

.r4vn_view_notes <- function(notes) {
  if (is.null(notes) || !length(notes)) return("")
  notes <- notes[nzchar(as.character(notes))]
  if (!length(notes)) return("")
  paste0('<section class="r4vn-notes"><h2>Notes</h2><ul>',
         paste0('<li>', .r4vn_view_escape(notes), '</li>', collapse = ""), '</ul></section>')
}

.r4vn_view_css <- function() {
  paste0(
    'html,body{margin:0;padding:0;background:#f6f8fb;color:#1f2937;font-family:-apple-system,BlinkMacSystemFont,"Segoe UI",Roboto,Arial,sans-serif;}',
    '.r4vn-page{max-width:1180px;margin:0 auto;padding:28px 30px 46px;}',
    '.r4vn-brand{font-size:12px;font-weight:700;letter-spacing:.12em;text-transform:uppercase;color:#64748b;margin-bottom:8px;}',
    'h1{font-size:26px;line-height:1.2;margin:0 0 7px;color:#0f172a;font-weight:700;}',
    '.r4vn-subtitle{font-size:13px;color:#64748b;margin-bottom:22px;}',
    '.r4vn-section{background:#fff;border:1px solid #e2e8f0;border-radius:10px;padding:18px 20px;margin:0 0 16px;box-shadow:0 1px 2px rgba(15,23,42,.03);}',
    '.r4vn-section h2,.r4vn-notes h2{font-size:16px;margin:0 0 12px;color:#0f172a;font-weight:650;}',
    '.section-note{font-size:12px;color:#64748b;margin:-5px 0 12px;}',
    '.r4vn-table-scroll{overflow-x:auto;}',
    '.r4vn-result-table{width:100%;border-collapse:collapse;font-size:13px;line-height:1.42;}',
    '.r4vn-result-table th{background:#f8fafc;color:#334155;padding:8px 10px;border-top:1px solid #cbd5e1;border-bottom:1px solid #cbd5e1;font-weight:650;white-space:nowrap;}',
    '.r4vn-result-table th.stub{text-align:left;}',
    '.r4vn-result-table th.value{text-align:right;font-variant-numeric:tabular-nums;}',
    '.r4vn-result-table td{padding:7px 10px;border-bottom:1px solid #eef2f7;vertical-align:top;}',
    '.r4vn-result-table td.value{text-align:right;font-variant-numeric:tabular-nums;}',
    '.r4vn-result-table td.stub{text-align:left;}',
    '.r4vn-result-table.compact td,.r4vn-result-table.compact th{padding:6px 8px;}',
    '.r4vn-result-table tbody tr:last-child td{border-bottom:0;}',
    '.r4vn-result-table strong{font-weight:700;color:#0f172a;}',
    '.kv-grid{display:grid;grid-template-columns:repeat(auto-fit,minmax(180px,1fr));gap:1px;background:#e2e8f0;border:1px solid #e2e8f0;border-radius:8px;overflow:hidden;}',
    '.kv-item{background:#fff;padding:9px 11px;min-width:0;}',
    '.kv-key{font-size:11px;color:#64748b;margin-bottom:3px;}',
    '.kv-value{font-size:13px;color:#0f172a;font-weight:600;overflow-wrap:anywhere;}',
    '.r4vn-notes{font-size:12px;color:#475569;background:#fff;border:1px solid #e2e8f0;border-radius:10px;padding:15px 20px;}',
    '.r4vn-notes ul{margin:0;padding-left:20px;}.r4vn-notes li{margin:4px 0;}',
    '.r4vn-empty{font-size:13px;color:#64748b;font-style:italic;}',
    '.r4vn-meta{display:flex;gap:12px;flex-wrap:wrap;margin:0 0 18px;}.r4vn-chip{font-size:12px;background:#eef2f7;border-radius:999px;padding:5px 9px;color:#475569;}',
    '@media(max-width:700px){.r4vn-page{padding:18px 12px 32px}.r4vn-section{padding:14px 12px}h1{font-size:22px}}'
  )
}

.r4vn_view_document <- function(title, body, notes = NULL, subtitle = NULL, prefix = "r4vn-result-") {
  html <- paste0('<!DOCTYPE html><html><head><meta charset="UTF-8"><meta name="viewport" content="width=device-width,initial-scale=1">',
                 '<style>', .r4vn_view_css(), '</style></head><body><main class="r4vn-page">',
                 '<div class="r4vn-brand">R4VN</div><h1>', .r4vn_view_escape(title), '</h1>',
                 if (is.null(subtitle) || !nzchar(subtitle)) '' else paste0('<div class="r4vn-subtitle">', .r4vn_view_escape(subtitle), '</div>'),
                 body, .r4vn_view_notes(notes), '</main></body></html>')
  file <- tempfile(pattern = prefix, fileext = ".html")
  writeLines(html, file, useBytes = TRUE)
  list(html = html, file = file)
}

.r4vn_view_open <- function(file) {
  # Never launch an external Viewer/browser during R CMD check, examples,
  # knitr rendering, or other non-interactive package work.
  if (!interactive()) return(invisible(file))
  normalized <- normalizePath(file, winslash = "/", mustWork = FALSE)
  viewer <- getOption("viewer")
  if (is.function(viewer)) viewer(normalized) else utils::browseURL(normalized)
  invisible(normalized)
}

.r4vn_stat_default_viewer <- function(x) {
  blocks <- paste0(vapply(names(x$sections), function(nm) {
    .r4vn_view_section(nm, x$sections[[nm]])
  }, character(1)), collapse = "")
  .r4vn_view_document(x$title, blocks, notes = x$notes, subtitle = "R4VN statistical result")
}

.r4vn_show <- function(x, show = TRUE,
                       console = get0("console", envir = parent.frame(), inherits = FALSE, ifnotfound = FALSE),
                       renderer = NULL) {
  # Any R4VN result that carries an underlying fitted model becomes the
  # active model for postestimation. This is deliberately independent of
  # Viewer/Console display so `show = FALSE` still updates model state.
  active_fit <- NULL
  if (is.list(x) && is.list(x$raw) && !is.null(x$raw$model)) active_fit <- x$raw$model
  # tabsurv()/cox() use a richer survival object rather than r4vn_stat. Prefer
  # the final multivariable Cox fit when present so margins/predict/lincom can
  # follow a direct cox() call just like they follow logistic()/poisson().
  if (is.null(active_fit) && inherits(x, "r4vn_surv") && is.list(x$cox) && !is.null(x$cox$multi_fit)) {
    active_fit <- x$cox$multi_fit
  }
  if (!is.null(active_fit) && exists(".r4vn_set_active_model", mode = "function", inherits = TRUE)) {
    try(.r4vn_set_active_model(active_fit, result = x), silent = TRUE)
  }
  if (isTRUE(show)) {
    if (!is.function(renderer)) {
      fn <- NULL
      if (!is.null(x$call) && length(x$call)) {
        head <- x$call[[1L]]
        fn <- tryCatch({
          if (is.symbol(head)) as.character(head)
          else if (is.call(head) && length(head) >= 3L && as.character(head[[1L]]) %in% c("::", ":::")) as.character(head[[3L]])
          else { z <- as.character(head); if (length(z)) tail(z, 1L) else NULL }
        }, error = function(e) NULL)
      }
      if (!is.null(fn) && length(fn) == 1L) {
        candidate <- paste0(".r4vn_viewer_", fn)
        if (exists(candidate, mode = "function", inherits = TRUE)) renderer <- get(candidate, mode = "function", inherits = TRUE)
      }
      if (!is.function(renderer) && inherits(x, "r4vn_stat")) renderer <- .r4vn_stat_default_viewer
    }
    if (is.function(renderer)) {
      rendered <- renderer(x)
      if (is.list(rendered)) {
        if (!is.null(rendered$html)) x$html <- rendered$html
        if (!is.null(rendered$file)) x$file <- rendered$file
        if (!is.null(rendered$table_html)) x$table_html <- rendered$table_html
        if (!is.null(rendered$file)) .r4vn_view_open(rendered$file)
      }
    }
  }
  if (isTRUE(console)) print(x)
  invisible(x)
}

.r4vn_check_counts <- function(x, arg = "counts") {
  if (!is.numeric(x) || any(!is.finite(x)) || any(x < 0) || any(abs(x - round(x)) > sqrt(.Machine$double.eps))) {
    stop(sprintf("`%s` must contain non-negative integer counts.", arg), call. = FALSE)
  }
  invisible(TRUE)
}

# ============================================================================
# Consolidated from: 010-tab-core.R
# ============================================================================
# Internal R4VN helpers. Not exported.

.r4vn_matrix_input <- function(..., row.names = NULL, col.names = NULL) {
  z <- list(...)
  if (length(z) == 1L && (is.matrix(z[[1L]]) || is.table(z[[1L]]))) {
    m <- as.matrix(z[[1L]])
  } else {
    if (!length(z) || !all(vapply(z, is.numeric, logical(1)))) stop("Supply numeric row vectors or one matrix/table.", call. = FALSE)
    lens <- vapply(z, length, integer(1))
    if (length(unique(lens)) != 1L) stop("All row vectors must have the same length.", call. = FALSE)
    m <- do.call(rbind, z)
  }
  .r4vn_check_counts(m)
  storage.mode(m) <- "double"
  if (is.null(rownames(m))) rownames(m) <- row.names %||% paste0("Row ", seq_len(nrow(m)))
  if (is.null(colnames(m))) colnames(m) <- col.names %||% paste0("Col ", seq_len(ncol(m)))
  m
}

.r4vn_percent_option <- function(percent, row = FALSE, col = FALSE, cell = FALSE, total = FALSE) {
  flags <- c(row = isTRUE(row), col = isTRUE(col), total = isTRUE(cell) || isTRUE(total))
  if (sum(flags) > 1L) stop("Only one of `row`, `col`, `cell`/`total` may be TRUE.", call. = FALSE)
  if (sum(flags) == 1L) return(names(flags)[which(flags)])
  match.arg(percent, c("none", "row", "col", "total"))
}

.r4vn_tab_calc <- function(m, percent = c("none", "row", "col", "total"), exp = FALSE,
                           chi = TRUE, fisher = FALSE, lr = FALSE, residual = FALSE,
                           adjresidual = FALSE, correct = FALSE, digits = 1L,
                           p_digits = 3L, workspace = 2e5) {
  percent <- match.arg(percent)
  chi_obj <- if (nrow(m) > 1L && ncol(m) > 1L) suppressWarnings(stats::chisq.test(m, correct = correct)) else NULL
  pct <- switch(percent,
    none = matrix(NA_real_, nrow(m), ncol(m)),
    row = 100 * prop.table(m, 1L),
    col = 100 * prop.table(m, 2L),
    total = 100 * prop.table(m)
  )
  body <- matrix("", nrow(m), ncol(m), dimnames = dimnames(m))
  for (i in seq_len(nrow(m))) for (j in seq_len(ncol(m))) {
    body[i, j] <- if (percent == "none") format(m[i, j], trim = TRUE, scientific = FALSE) else paste0(format(m[i, j], trim = TRUE, scientific = FALSE), " (", .r4vn_num(pct[i, j], digits), "%)")
  }
  body <- cbind(body, Total = format(rowSums(m), trim = TRUE, scientific = FALSE))
  body <- rbind(body, Total = c(format(colSums(m), trim = TRUE, scientific = FALSE), format(sum(m), trim = TRUE, scientific = FALSE)))
  tests <- list()
  if (isTRUE(chi) && !is.null(chi_obj)) tests[[length(tests) + 1L]] <- data.frame(Test = if (correct && all(dim(m) == 2L)) "Pearson chi-square with Yates correction" else "Pearson chi-square", Statistic = .r4vn_num(unname(chi_obj$statistic), 3), df = .r4vn_num(unname(chi_obj$parameter), 0), p = .r4vn_p(chi_obj$p.value, p_digits), stringsAsFactors = FALSE)
  if (isTRUE(lr) && !is.null(chi_obj)) {
    e <- chi_obj$expected
    pos <- m > 0 & e > 0
    g2 <- 2 * sum(m[pos] * log(m[pos] / e[pos]))
    df <- (nrow(m) - 1L) * (ncol(m) - 1L)
    tests[[length(tests) + 1L]] <- data.frame(Test = "Likelihood-ratio chi-square", Statistic = .r4vn_num(g2, 3), df = df, p = .r4vn_p(stats::pchisq(g2, df, lower.tail = FALSE), p_digits), stringsAsFactors = FALSE)
  }
  if (isTRUE(fisher) && nrow(m) > 1L && ncol(m) > 1L) {
    ft <- tryCatch(stats::fisher.test(m, workspace = workspace), error = function(e) NULL)
    tests[[length(tests) + 1L]] <- data.frame(Test = "Fisher exact", Statistic = "", df = "", p = if (is.null(ft)) "not computed" else .r4vn_p(ft$p.value, p_digits), stringsAsFactors = FALSE)
  }
  if (!is.null(chi_obj)) {
    n <- sum(m); k <- min(nrow(m) - 1L, ncol(m) - 1L)
    v <- if (n > 0 && k > 0) sqrt(unname(chi_obj$statistic) / (n * k)) else NA_real_
    phi <- if (all(dim(m) == 2L) && n > 0) sqrt(unname(chi_obj$statistic) / n) else NA_real_
    tests[[length(tests) + 1L]] <- data.frame(Test = if (all(dim(m) == 2L)) "Phi / Cramer's V" else "Cramer's V", Statistic = .r4vn_num(if (all(dim(m) == 2L)) phi else v, 3), df = "", p = "", stringsAsFactors = FALSE)
  }
  sections <- list("Observed" = body)
  if (isTRUE(exp) && !is.null(chi_obj)) sections[["Expected counts"]] <- apply(chi_obj$expected, 2L, .r4vn_num, digits = digits)
  if (isTRUE(residual) && !is.null(chi_obj)) sections[["Pearson residuals"]] <- apply(chi_obj$residuals, 2L, .r4vn_num, digits = 2L)
  if (isTRUE(adjresidual) && !is.null(chi_obj)) {
    rs <- rowSums(m); cs <- colSums(m); n <- sum(m)
    den <- sqrt(chi_obj$expected * outer(1 - rs / n, 1 - cs / n))
    ar <- (m - chi_obj$expected) / den
    sections[["Adjusted residuals"]] <- apply(ar, 2L, .r4vn_num, digits = 2L)
  }
  if (length(tests)) sections[["Tests and association"]] <- do.call(rbind, tests)
  list(sections = sections, raw = list(counts = m, percentages = pct, expected = if (is.null(chi_obj)) NULL else chi_obj$expected, chi_square = chi_obj))
}

.r4vn_tab_data <- function(row_expr, col_expr = NULL, data = NULL, env = parent.frame(),
                           percent = c("none", "row", "col", "total"),
                           row = FALSE, col = FALSE, cell = FALSE, total = FALSE, exp = FALSE,
                           chi = TRUE, fisher = FALSE, lr = FALSE, residual = FALSE,
                           adjresidual = FALSE, correct = FALSE, missing = c("no", "ifany", "always"),
                           digits = 1, p_digits = 3, workspace = 2e5, show = TRUE) {
  data <- .r4vn_stat_data(data); missing <- match.arg(missing)
  percent <- .r4vn_percent_option(percent, row, col, cell, total)
  x <- .r4vn_eval_var(row_expr, data, env, "row variable")
  if (missing == "no") x <- ifelse(is.na(x), NA, as.character(x))
  else if (missing == "always" || (missing == "ifany" && anyNA(x))) x <- base::addNA(factor(x), ifany = missing != "always")
  if (is.null(col_expr) || identical(col_expr, quote(NULL))) {
    m <- table(x, useNA = if (missing == "no") "no" else "ifany")
    pct <- 100 * m / sum(m)
    out <- data.frame(Level = names(m), n = as.vector(m), Percent = .r4vn_num(as.vector(pct), digits), stringsAsFactors = FALSE)
    return(.r4vn_show(.r4vn_result("One-way frequency table", list("Observed" = out), raw = list(counts = m), call = match.call()), show))
  }
  y <- .r4vn_eval_var(col_expr, data, env, "column variable")
  keep <- if (missing == "no") !is.na(x) & !is.na(y) else rep(TRUE, length(x))
  if (missing != "no") {
    if (missing == "always" || anyNA(x)) x <- base::addNA(factor(x), ifany = missing != "always")
    if (missing == "always" || anyNA(y)) y <- base::addNA(factor(y), ifany = missing != "always")
  }
  m <- table(x[keep], y[keep], useNA = if (missing == "no") "no" else "ifany")
  z <- .r4vn_tab_calc(m, percent, exp, chi, fisher, lr, residual, adjresidual, correct, digits, p_digits, workspace)
  .r4vn_show(.r4vn_result("Two-way frequency table", z$sections, raw = z$raw, call = match.call()), show)
}

# ============================================================================
# Consolidated from: 020-ttest-core.R
# ============================================================================
# Internal R4VN helpers. Not exported.

.r4vn_t_group <- function(n, mean, sd, level = 0.95) {
  se <- sd / sqrt(n); crit <- stats::qt(1 - (1 - level) / 2, n - 1)
  c(n = n, mean = mean, se = se, sd = sd, lower = mean - crit * se, upper = mean + crit * se)
}

# ============================================================================
# Consolidated from: 030-anova-core.R
# ============================================================================
# Internal R4VN helpers. Not exported.

.r4vn_summary_group <- function(g) {
  if (is.numeric(g) && length(g) == 3L) {
    if (!is.null(names(g)) && all(c("n", "mean", "sd") %in% names(g))) return(unname(g[c("n", "mean", "sd")]))
    return(unname(g))
  }
  if (is.list(g) && all(c("n", "mean", "sd") %in% names(g))) return(unlist(g[c("n", "mean", "sd")], use.names = FALSE))
  stop("Each group must be `c(n, mean, sd)` or a list with n, mean, and sd.", call. = FALSE)
}

.r4vn_anova_data <- function(x_expr, by_expr, data = NULL, env = parent.frame(), bartlett = TRUE,
                             level = 0.95, digits = 3, p_digits = 3, show = TRUE) {
  data <- .r4vn_stat_data(data); x <- .r4vn_eval_var(x_expr, data, env, "outcome"); g <- .r4vn_eval_var(by_expr, data, env, "by")
  if (!is.numeric(x)) stop("The ANOVA outcome must be numeric.", call. = FALSE)
  ok <- stats::complete.cases(x, g); x <- x[ok]; g <- droplevels(factor(g[ok])); if (nlevels(g) < 2L) stop("`by` must have at least two groups.", call. = FALSE)
  s <- split(x, g); args <- lapply(s, function(z) c(n = length(z), mean = mean(z), sd = stats::sd(z)))
  do.call(anovai, c(args, list(group.names = names(s), bartlett = bartlett, level = level, digits = digits, p_digits = p_digits, show = show)))
}

# ============================================================================
# Consolidated from: 040-proportion-core.R
# ============================================================================
# Internal R4VN helpers. Not exported.

.r4vn_prop_ci <- function(events, total, level = 0.95, method = c("exact", "wilson", "wald")) {
  method <- match.arg(method); p <- events / total; alpha <- 1 - level; z <- stats::qnorm(1 - alpha / 2)
  if (method == "exact") return(unname(stats::binom.test(events, total, conf.level = level)$conf.int))
  if (method == "wald") return(pmax(0, pmin(1, p + c(-1, 1) * z * sqrt(p * (1 - p) / total))))
  den <- 1 + z^2 / total; center <- (p + z^2 / (2 * total)) / den; half <- z * sqrt(p * (1 - p) / total + z^2 / (4 * total^2)) / den
  c(center - half, center + half)
}

# ============================================================================
# Consolidated from: 050-epi-core.R
# ============================================================================
# Internal R4VN helpers. Not exported.

.r4vn_as_2x2 <- function(x, name = NULL) {
  supplied_dimnames <- NULL

  if (is.matrix(x) || is.table(x)) {
    m <- as.matrix(x)
    supplied_dimnames <- dimnames(m)
  } else if (is.numeric(x) && length(x) == 4L) {
    m <- matrix(x, 2L, 2L, byrow = TRUE)
  } else {
    stop(
      "Each epidemiological table must be a 2 x 2 matrix or four counts.",
      call. = FALSE
    )
  }

  if (!all(dim(m) == 2L)) {
    stop("A 2 x 2 table is required.", call. = FALSE)
  }

  .r4vn_check_counts(m)

  # Preserve meaningful labels supplied by epi() or by a user-provided
  # matrix/table. Generic epidemiological labels are used only when no
  # usable row/column labels are available.
  row_labels <- if (
    !is.null(supplied_dimnames) &&
      length(supplied_dimnames) >= 1L &&
      !is.null(supplied_dimnames[[1L]]) &&
      length(supplied_dimnames[[1L]]) == 2L
  ) {
    as.character(supplied_dimnames[[1L]])
  } else {
    c("Exposed", "Unexposed")
  }

  col_labels <- if (
    !is.null(supplied_dimnames) &&
      length(supplied_dimnames) >= 2L &&
      !is.null(supplied_dimnames[[2L]]) &&
      length(supplied_dimnames[[2L]]) == 2L
  ) {
    as.character(supplied_dimnames[[2L]])
  } else {
    c("Case", "Noncase")
  }

  dn_names <- if (!is.null(supplied_dimnames)) names(supplied_dimnames) else NULL

  row_dimension <- if (
    !is.null(dn_names) &&
      length(dn_names) >= 1L &&
      !is.na(dn_names[1L]) &&
      nzchar(dn_names[1L])
  ) {
    dn_names[1L]
  } else {
    "Exposure"
  }

  col_dimension <- if (
    !is.null(dn_names) &&
      length(dn_names) >= 2L &&
      !is.na(dn_names[2L]) &&
      nzchar(dn_names[2L])
  ) {
    dn_names[2L]
  } else {
    "Outcome"
  }

  dimnames(m) <- stats::setNames(
    list(row_labels, col_labels),
    c(row_dimension, col_dimension)
  )

  storage.mode(m) <- "double"
  m
}

.r4vn_epi_single <- function(m, level = 0.95, correction = 0.5) {
  a <- m[1,1]; b <- m[1,2]; c <- m[2,1]; d <- m[2,2]; n1 <- a + b; n0 <- c + d; n <- sum(m); z <- stats::qnorm(1 - (1 - level) / 2)
  if (n <= 0 || n1 <= 0 || n0 <= 0 || any(colSums(m) <= 0)) stop("The 2 x 2 table must have positive row and column totals.", call. = FALSE)
  p1 <- a / n1; p0 <- c / n0; pt <- (a + c) / n; or <- (a * d) / (b * c); rr <- p1 / p0; rd <- p1 - p0
  mc <- if (any(m == 0)) m + correction else m; aa <- mc[1,1]; bb <- mc[1,2]; cc <- mc[2,1]; dd <- mc[2,2]
  se_or <- sqrt(1/aa + 1/bb + 1/cc + 1/dd); ci_or <- exp(log((aa*dd)/(bb*cc)) + c(-1,1) * z * se_or)
  se_rr <- sqrt(1/aa - 1/(aa+bb) + 1/cc - 1/(cc+dd)); ci_rr <- exp(log((aa/(aa+bb))/(cc/(cc+dd))) + c(-1,1) * z * se_rr)
  se_rd <- sqrt(p1*(1-p1)/n1 + p0*(1-p0)/n0); ci_rd <- rd + c(-1,1)*z*se_rd
  afe <- if (is.finite(rr) && rr >= 1) (rr - 1) / rr else if (is.finite(rr)) 1 - rr else NA_real_; par <- pt - p0; paf <- if (pt > 0) par / pt else NA_real_; nnt <- if (rd == 0) Inf else 1 / abs(rd)
  chi <- suppressWarnings(stats::chisq.test(m, correct = FALSE)); fisher <- stats::fisher.test(m)
  list(a=a,b=b,c=c,d=d,risk.exposed=p1,risk.unexposed=p0,risk.total=pt,OR=or,OR.ci=ci_or,RR=rr,RR.ci=ci_rr,PR=rr,PR.ci=ci_rr,RD=rd,RD.ci=ci_rd,AF.exposed=afe,PAR=par,PAF=paf,NNT=nnt,chi.square=unname(chi$statistic),chi.p=chi$p.value,fisher.p=fisher$p.value,zero.correction=any(m==0))
}

.r4vn_strata <- function(by) {
  if (is.array(by) && length(dim(by)) == 3L && all(dim(by)[1:2] == 2L)) {
    out <- lapply(seq_len(dim(by)[3]), function(i) .r4vn_as_2x2(by[,,i])); names(out) <- dimnames(by)[[3]] %||% paste0("Stratum ", seq_along(out)); return(out)
  }
  if (!is.list(by)) stop("`by` must be a list of 2 x 2 tables or a 2 x 2 x K array.", call. = FALSE)
  out <- lapply(by, .r4vn_as_2x2); if (is.null(names(out))) names(out) <- paste0("Stratum ", seq_along(out)); out
}

# ============================================================================
# Consolidated from: 060-model-core.R
# ============================================================================
# Internal helpers for correlation and regression models.

.r4vn_deparse1 <- function(x) paste(deparse(x, width.cutoff = 500L), collapse = "")

.r4vn_is_formula_expr <- function(x) {
  inherits(x, "formula") || (is.call(x) && identical(x[[1L]], as.name("~")))
}

.r4vn_build_formula <- function(lhs_expr, rhs_exprs = list(), env = parent.frame(), noconstant = FALSE) {
  lhs_value <- tryCatch(eval(lhs_expr, envir = env), error = function(e) NULL)
  if (.r4vn_is_formula_expr(lhs_expr) || inherits(lhs_value, "formula")) {
    if (length(rhs_exprs)) stop("Do not supply additional unnamed predictors when the first argument is a formula.", call. = FALSE)
    f <- if (inherits(lhs_value, "formula")) lhs_value else if (inherits(lhs_expr, "formula")) lhs_expr else stats::as.formula(lhs_expr, env = env)
    if (noconstant) f <- stats::update.formula(f, . ~ . - 1)
    return(f)
  }
  lhs <- .r4vn_deparse1(lhs_expr)
  rhs <- if (length(rhs_exprs)) vapply(rhs_exprs, .r4vn_deparse1, character(1)) else "1"
  rhs <- paste(rhs, collapse = " + ")
  if (noconstant) rhs <- paste0(rhs, " - 1")
  stats::as.formula(paste(lhs, "~", rhs), env = env)
}

.r4vn_apply_ref <- function(data, ref = NULL) {
  if (is.null(ref)) return(data)
  if (!is.list(ref) || is.null(names(ref)) || any(!nzchar(names(ref)))) stop("`ref` must be a named list, for example `list(sex = \"Female\")`.", call. = FALSE)
  for (nm in names(ref)) {
    if (!nm %in% names(data)) stop(sprintf("Reference variable `%s` was not found in data.", nm), call. = FALSE)
    data[[nm]] <- stats::relevel(factor(data[[nm]]), ref = as.character(ref[[nm]])[1L])
  }
  data
}

.r4vn_eval_subset <- function(expr, data, env) {
  if (is.null(expr) || identical(expr, quote(NULL))) return(rep(TRUE, nrow(data)))
  z <- eval(expr, envir = data, enclos = env)
  if (!is.logical(z) || length(z) != nrow(data)) stop("`subset` must evaluate to one logical value per row.", call. = FALSE)
  z[is.na(z)] <- FALSE
  z
}

.r4vn_eval_optional <- function(expr, data, env, arg) {
  if (is.null(expr) || identical(expr, quote(NULL))) return(NULL)
  .r4vn_eval_var(expr, data, env, arg)
}

.r4vn_prepare_model_data <- function(data = NULL, env = parent.frame(), subset_expr = quote(NULL),
                                     weights_expr = quote(NULL), cluster_expr = quote(NULL),
                                     exposure_expr = quote(NULL), offset_expr = quote(NULL), ref = NULL) {
  # Evaluate all optional vectors against the original data first, then apply
  # the same subset to every vector. This keeps data, weights, clusters,
  # exposure, and offset perfectly aligned, including when the user supplies
  # an external vector rather than a column name.
  d0 <- .r4vn_apply_ref(.r4vn_stat_data(data), ref)
  keep <- .r4vn_eval_subset(subset_expr, d0, env)
  w0 <- .r4vn_eval_optional(weights_expr, d0, env, "weights")
  cl0 <- .r4vn_eval_optional(cluster_expr, d0, env, "cluster")
  ex0 <- .r4vn_eval_optional(exposure_expr, d0, env, "exposure")
  off0 <- .r4vn_eval_optional(offset_expr, d0, env, "offset")

  d <- d0[keep, , drop = FALSE]
  rownames(d) <- seq_len(nrow(d))
  w <- if (is.null(w0)) NULL else w0[keep]
  cl <- if (is.null(cl0)) NULL else cl0[keep]
  ex <- if (is.null(ex0)) NULL else ex0[keep]
  off <- if (is.null(off0)) NULL else off0[keep]

  if (!is.null(w) && (!is.numeric(w) || any(!is.finite(w)) || any(w < 0))) {
    stop("`weights` must be finite and non-negative.", call. = FALSE)
  }
  if (!is.null(ex) && (!is.numeric(ex) || any(!is.finite(ex)) || any(ex <= 0))) {
    stop("`exposure` must be finite and strictly positive.", call. = FALSE)
  }
  if (!is.null(off) && (!is.numeric(off) || any(!is.finite(off)))) {
    stop("`offset` must be finite numeric values.", call. = FALSE)
  }
  list(data = d, weights = w, cluster = cl, exposure = ex, offset = off)
}

.r4vn_model_formula <- function(formula, n, data_names = character(), weights = NULL, offset = NULL) {
  # model.frame() evaluates weights and offsets in the data/formula context.
  # Give the formula a small child environment containing concrete vectors,
  # instead of passing local expressions such as `prep$weights`. This avoids
  # "object 'prep' not found" and does not add temporary columns to `data`, so
  # formulas containing `.` remain correct.
  weight_name <- ".r4vn_internal_weights_7e4f9c"
  offset_name <- ".r4vn_internal_offset_7e4f9c"
  hit <- intersect(c(weight_name, offset_name), data_names)
  if (length(hit)) stop(sprintf("Please rename the reserved internal variable `%s`.", hit[1L]), call. = FALSE)
  old_env <- environment(formula)
  if (is.null(old_env)) old_env <- parent.frame()
  fit_env <- new.env(parent = old_env)
  fit_env[[weight_name]] <- if (is.null(weights)) rep(1, n) else weights
  fit_env[[offset_name]] <- if (is.null(offset)) rep(0, n) else offset
  environment(formula) <- fit_env
  formula
}

.r4vn_used_rows <- function(fit, n) {
  keep <- rep(TRUE, n)
  if (!is.null(fit$na.action)) keep[as.integer(fit$na.action)] <- FALSE
  keep
}

.r4vn_safe_inverse <- function(x) {
  out <- tryCatch(solve(x), error = function(e) NULL)
  if (is.null(out)) out <- tryCatch(qr.solve(x, diag(nrow(x)), tol = 1e-10), error = function(e) NULL)
  if (is.null(out)) stop("The covariance matrix could not be inverted; the model may be singular.", call. = FALSE)
  out
}

.r4vn_model_vcov <- function(fit, vce = c("model", "robust", "cluster"), cluster = NULL) {
  vce <- match.arg(vce)
  b <- stats::coef(fit)
  good <- !is.na(b)
  nm <- names(b)
  full <- matrix(NA_real_, length(b), length(b), dimnames = list(nm, nm))
  if (!any(good)) return(full)
  if (vce == "model") {
    v <- tryCatch(stats::vcov(fit, complete = TRUE), error = function(e) stats::vcov(fit))
    full[rownames(v), colnames(v)] <- v
    return(full)
  }
  X <- stats::model.matrix(fit)[, good, drop = FALSE]
  n <- nrow(X); k <- ncol(X)
  if (n <= k) stop("Robust covariance requires more observations than estimable coefficients.", call. = FALSE)
  if (inherits(fit, "glm")) {
    mu <- fit$fitted.values
    eta <- fit$linear.predictors
    y <- fit$y
    prior <- fit$prior.weights
    mu_eta <- fit$family$mu.eta(eta)
    var_mu <- pmax(fit$family$variance(mu), .Machine$double.eps)
    wi <- prior * (mu_eta^2 / var_mu)
    score_scalar <- prior * (y - mu) * mu_eta / var_mu
    bread <- .r4vn_safe_inverse(crossprod(X, wi * X))
    score <- X * as.numeric(score_scalar)
  } else {
    w <- stats::weights(fit)
    if (is.null(w)) w <- rep(1, n)
    u <- stats::residuals(fit)
    bread <- .r4vn_safe_inverse(crossprod(X, w * X))
    score <- X * as.numeric(w * u)
  }
  if (vce == "robust") {
    meat <- crossprod(score) * n / (n - k)
  } else {
    if (is.null(cluster)) stop("Supply `cluster` when `vce = \"cluster\"`.", call. = FALSE)
    if (anyNA(cluster)) stop("The cluster variable may not contain missing values among fitted observations.", call. = FALSE)
    cluster <- droplevels(factor(cluster))
    if (length(cluster) != n) stop("The cluster variable is not aligned with the fitted observations.", call. = FALSE)
    G <- nlevels(cluster)
    if (G < 2L) stop("At least two clusters are required.", call. = FALSE)
    summed <- rowsum(score, cluster, reorder = FALSE)
    meat <- crossprod(summed) * (G / (G - 1)) * ((n - 1) / (n - k))
  }
  v <- bread %*% meat %*% bread
  full[good, good] <- v
  full
}

.r4vn_coef_raw <- function(fit, vcov, level = 0.95, distribution = c("z", "t")) {
  distribution <- match.arg(distribution)
  b <- stats::coef(fit)
  se <- sqrt(diag(vcov))
  stat <- b / se
  if (distribution == "t") {
    df <- stats::df.residual(fit)
    p <- 2 * stats::pt(abs(stat), df = df, lower.tail = FALSE)
    crit <- stats::qt(1 - (1 - level) / 2, df = df)
  } else {
    p <- 2 * stats::pnorm(abs(stat), lower.tail = FALSE)
    crit <- stats::qnorm(1 - (1 - level) / 2)
  }
  list(term = names(b), estimate = unname(b), se = unname(se), statistic = unname(stat),
       p.value = unname(p), lower = unname(b - crit * se), upper = unname(b + crit * se))
}

.r4vn_coef_table <- function(raw, digits = 3, p_digits = 3, statistic = c("z", "t"),
                             transform = FALSE, label = "Estimate") {
  statistic <- match.arg(statistic)
  est <- raw$estimate; se <- raw$se; lo <- raw$lower; hi <- raw$upper
  if (transform) {
    est0 <- est
    est <- base::exp(est0); se <- base::exp(est0) * se
    lo <- base::exp(lo); hi <- base::exp(hi)
  }
  term <- raw$term
  term[term == "(Intercept)"] <- "_cons"
  out <- data.frame(Term = term, stringsAsFactors = FALSE, check.names = FALSE)
  out[[label]] <- .r4vn_num(est, digits)
  out[["Std. err."]] <- .r4vn_num(se, digits)
  out[[statistic]] <- .r4vn_num(raw$statistic, 2)
  out[[if (statistic == "t") "P>|t|" else "P>|z|"]] <- .r4vn_p(raw$p.value, p_digits)
  out[["Lower"]] <- .r4vn_num(lo, digits)
  out[["Upper"]] <- .r4vn_num(hi, digits)
  out
}

.r4vn_wald_overall <- function(fit, vcov, linear = FALSE) {
  b <- stats::coef(fit)
  use <- !is.na(b) & names(b) != "(Intercept)"
  q <- sum(use)
  if (!q) return(list(statistic = NA_real_, df1 = 0L, df2 = stats::df.residual(fit), p.value = NA_real_))
  V <- vcov[use, use, drop = FALSE]
  W <- as.numeric(t(b[use]) %*% .r4vn_safe_inverse(V) %*% b[use])
  if (linear) {
    F <- W / q
    df2 <- stats::df.residual(fit)
    list(statistic = F, df1 = q, df2 = df2, p.value = stats::pf(F, q, df2, lower.tail = FALSE))
  } else list(statistic = W, df1 = q, df2 = NA_integer_, p.value = stats::pchisq(W, q, lower.tail = FALSE))
}

.r4vn_vif <- function(fit, digits = 3) {
  X <- stats::model.matrix(fit)
  if ("(Intercept)" %in% colnames(X)) X <- X[, colnames(X) != "(Intercept)", drop = FALSE]
  if (ncol(X) < 2L) return(NULL)
  out <- vapply(seq_len(ncol(X)), function(j) {
    y <- X[, j]; z <- X[, -j, drop = FALSE]
    if (stats::var(y) == 0) return(NA_real_)
    fj <- stats::lm.fit(cbind(1, z), y)
    tss <- sum((y - mean(y))^2)
    r2 <- if (tss > 0) 1 - sum(fj$residuals^2) / tss else NA_real_
    1 / (1 - r2)
  }, numeric(1))
  data.frame(Term = colnames(X), VIF = .r4vn_num(out, digits), stringsAsFactors = FALSE)
}

.r4vn_glm_null <- function(fit) {
  y <- fit$y
  w <- fit$prior.weights
  off <- fit$offset
  tryCatch(stats::glm(y ~ 1, family = fit$family, weights = w, offset = off,
                      na.action = stats::na.omit, model = TRUE, x = TRUE, y = TRUE), error = function(e) NULL)
}

.r4vn_logistic_gof <- function(fit, groups = 10L, digits = 3, p_digits = 3) {
  y <- fit$y; p <- fit$fitted.values
  uq <- unique(stats::quantile(p, probs = seq(0, 1, length.out = groups + 1L), na.rm = TRUE, type = 2))
  if (length(uq) < 4L) return(data.frame(Test = "Hosmer-Lemeshow", Chi.square = "", df = "", p = "not computed", stringsAsFactors = FALSE))
  g <- cut(p, breaks = uq, include.lowest = TRUE, labels = FALSE)
  obs1 <- rowsum(y, g, reorder = FALSE); exp1 <- rowsum(p, g, reorder = FALSE); n <- rowsum(rep(1, length(y)), g, reorder = FALSE)
  obs0 <- n - obs1; exp0 <- n - exp1
  stat <- sum((obs1 - exp1)^2 / pmax(exp1, .Machine$double.eps) + (obs0 - exp0)^2 / pmax(exp0, .Machine$double.eps))
  df <- length(obs1) - 2L
  data.frame(Test = "Hosmer-Lemeshow", Chi.square = .r4vn_num(stat, digits), df = df,
             p = .r4vn_p(stats::pchisq(stat, df, lower.tail = FALSE), p_digits), stringsAsFactors = FALSE)
}

.r4vn_classification <- function(fit, cutoff = 0.5, digits = 1) {
  y <- as.integer(fit$y > 0.5); pred <- as.integer(fit$fitted.values >= cutoff)
  tp <- sum(pred == 1 & y == 1); fp <- sum(pred == 1 & y == 0); tn <- sum(pred == 0 & y == 0); fn <- sum(pred == 0 & y == 1)
  tab <- matrix(c(tp, fp, fn, tn), nrow = 2, byrow = TRUE,
                dimnames = list("Predicted" = c("Event", "Nonevent"), "Observed" = c("Event", "Nonevent")))
  metric <- data.frame(Measure = c("Sensitivity", "Specificity", "Positive predictive value", "Negative predictive value", "Accuracy"),
                       Percent = .r4vn_num(100 * c(tp / max(1, tp + fn), tn / max(1, tn + fp), tp / max(1, tp + fp), tn / max(1, tn + fn), (tp + tn) / max(1, tp + fp + tn + fn)), digits),
                       stringsAsFactors = FALSE)
  list(table = tab, metrics = metric)
}