R/tab.R

In R4VN: Health Data Analysis and Publication-Ready Reporting

#============================================================
# Internal helpers for superby tables
#============================================================
.r4vn_superby_escape <- function(x) {
  x <- as.character(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)
  gsub("'", "&#39;", x, fixed = TRUE)
}

.r4vn_superby_levels <- function(x) {
  observed <- x[!is.na(x)]
  values <- if (is.factor(x)) levels(x) else if (is.logical(x)) c(FALSE, TRUE) else {
    z <- unique(observed)
    if (is.numeric(z)) sort(z) else z
  }
  values[as.character(values) %in% as.character(observed)]
}

.r4vn_superby_row_keys <- function(rows) {
  vapply(rows, function(z) paste(z$variable, z$type, z$item, z$kind, sep = "\034"), character(1))
}

.r4vn_superby_robust_vcov <- function(fit) {
  X <- stats::model.matrix(fit)
  mu <- stats::fitted(fit)
  score_residual <- fit$y - mu
  bread <- tryCatch(solve(crossprod(X, X * as.vector(mu))), error = function(e) NULL)
  if (is.null(bread)) return(NULL)
  meat <- crossprod(X, X * as.vector(score_residual^2))
  output <- bread %*% meat %*% bread
  dimnames(output) <- list(colnames(X), colnames(X))
  output
}

.r4vn_superby_prepare_value <- function(x, type) {
  if (identical(type, "categorical")) factor(x, levels = .r4vn_superby_levels(x)) else suppressWarnings(as.numeric(x))
}

.r4vn_superby_interaction_one <- function(data, output, superby_name, variable) {
  meta <- output$metadata
  meta_index <- match(variable, meta$variable)
  if (is.na(meta_index)) return(NA_real_)
  variable_type <- meta$type[meta_index]
  outcome_name <- output$by
  continuous_outcome <- identical(output$outcome_type, "continuous")
  if (is.null(outcome_name) || !outcome_name %in% names(data)) return(NA_real_)

  multi_meta <- output$multi
  adjusted_meta <- output$adjusted
  if (!is.null(multi_meta) && nrow(multi_meta) && variable %in% multi_meta$variable) {
    covariate_meta <- multi_meta[multi_meta$variable != variable, , drop = FALSE]
  } else if (!is.null(adjusted_meta) && nrow(adjusted_meta)) {
    covariate_meta <- adjusted_meta[adjusted_meta$variable != variable, , drop = FALSE]
  } else {
    covariate_meta <- meta[0, , drop = FALSE]
  }
  covariate_meta <- covariate_meta[!covariate_meta$variable %in% c(outcome_name, superby_name, variable), , drop = FALSE]
  if (nrow(covariate_meta)) covariate_meta <- covariate_meta[!duplicated(covariate_meta$variable), , drop = FALSE]

  if (continuous_outcome) {
    model_data <- data.frame(.outcome = suppressWarnings(as.numeric(data[[outcome_name]])), stringsAsFactors = FALSE)
  } else {
    outcome_levels <- output$by_levels
    if (length(outcome_levels) != 2L) return(NA_real_)
    event_level <- output$event
    model_data <- data.frame(.outcome = as.integer(as.character(data[[outcome_name]]) == as.character(event_level)), stringsAsFactors = FALSE)
  }
  super_levels <- .r4vn_superby_levels(data[[superby_name]])
  if (length(super_levels) < 2L) return(NA_real_)
  model_data$.super <- factor(data[[superby_name]], levels = super_levels)
  model_data$.x <- .r4vn_superby_prepare_value(data[[variable]], variable_type)

  z_names <- character()
  if (nrow(covariate_meta)) {
    for (j in seq_len(nrow(covariate_meta))) {
      z_name <- paste0(".z", j)
      z_names <- c(z_names, z_name)
      model_data[[z_name]] <- .r4vn_superby_prepare_value(data[[covariate_meta$variable[j]]], covariate_meta$type[j])
    }
  }
  keep <- stats::complete.cases(model_data)
  if (continuous_outcome) keep <- keep & is.finite(model_data$.outcome)
  if (!identical(variable_type, "categorical")) keep <- keep & is.finite(model_data$.x)
  if (length(z_names)) {
    for (z in z_names) if (is.numeric(model_data[[z]])) keep <- keep & is.finite(model_data[[z]])
  }
  model_data <- model_data[keep, , drop = FALSE]
  if (nrow(model_data) < 10L || nlevels(droplevels(model_data$.super)) < 2L) return(NA_real_)
  if (identical(variable_type, "categorical")) {
    model_data$.x <- droplevels(model_data$.x)
    if (nlevels(model_data$.x) < 2L) return(NA_real_)
  } else if (!is.finite(stats::sd(model_data$.x)) || stats::sd(model_data$.x) == 0) return(NA_real_)
  if (!continuous_outcome && length(unique(model_data$.outcome)) < 2L) return(NA_real_)

  reduced_terms <- c(".x", ".super", z_names)
  full_terms <- c(".x * .super", z_names)
  reduced_formula <- stats::as.formula(paste(".outcome ~", paste(reduced_terms, collapse = " + ")))
  full_formula <- stats::as.formula(paste(".outcome ~", paste(full_terms, collapse = " + ")))

  if (continuous_outcome) {
    reduced_fit <- tryCatch(stats::lm(reduced_formula, data = model_data), error = function(e) NULL)
    full_fit <- tryCatch(stats::lm(full_formula, data = model_data), error = function(e) NULL)
    if (is.null(reduced_fit) || is.null(full_fit)) return(NA_real_)
    comparison <- tryCatch(stats::anova(reduced_fit, full_fit), error = function(e) NULL)
    if (is.null(comparison) || nrow(comparison) < 2L || !"Pr(>F)" %in% names(comparison)) return(NA_real_)
    return(unname(comparison[["Pr(>F)"]][2L]))
  }

  family <- if (isTRUE(output$effect_type %in% c("RR", "PR"))) stats::poisson("log") else stats::binomial("logit")
  full_fit <- tryCatch(stats::glm(full_formula, family = family, data = model_data, y = TRUE), error = function(e) NULL)
  if (is.null(full_fit)) return(NA_real_)
  covariance <- if (isTRUE(output$effect_type %in% c("RR", "PR"))) .r4vn_superby_robust_vcov(full_fit) else tryCatch(stats::vcov(full_fit), error = function(e) NULL)
  if (is.null(covariance)) return(NA_real_)
  X <- stats::model.matrix(full_fit)
  assignment <- attr(X, "assign")
  term_labels <- attr(stats::terms(full_fit), "term.labels")
  interaction_position <- which(term_labels %in% c(".x:.super", ".super:.x"))
  if (!length(interaction_position)) return(NA_real_)
  interaction_terms <- colnames(X)[assignment %in% interaction_position]
  beta <- stats::coef(full_fit)
  interaction_terms <- intersect(interaction_terms, names(beta))
  interaction_terms <- interaction_terms[is.finite(beta[interaction_terms])]
  if (!length(interaction_terms)) return(NA_real_)
  V <- covariance[interaction_terms, interaction_terms, drop = FALSE]
  b <- beta[interaction_terms]
  valid <- is.finite(b) & is.finite(diag(V)) & diag(V) > 0
  b <- b[valid]
  V <- V[valid, valid, drop = FALSE]
  if (!length(b)) return(NA_real_)
  inverse <- tryCatch(solve(V), error = function(e) tryCatch(qr.solve(V), error = function(e2) NULL))
  if (is.null(inverse)) return(NA_real_)
  df <- qr(V)$rank
  if (df < 1L) return(NA_real_)
  statistic <- as.numeric(t(b) %*% inverse %*% b)
  if (!is.finite(statistic)) return(NA_real_)
  stats::pchisq(statistic, df = df, lower.tail = FALSE)
}

.r4vn_superby_interactions <- function(data, output, superby_name, enabled = TRUE) {
  result <- stats::setNames(rep(NA_real_, nrow(output$metadata)), output$metadata$variable)
  if (!isTRUE(enabled) || is.null(output$by)) return(list(show = FALSE, p = result, note = ""))
  continuous_outcome <- identical(output$outcome_type, "continuous")
  if (!continuous_outcome && length(output$by_levels) != 2L) {
    return(list(show = FALSE, p = result, note = "Interaction p-values are available only for binary or continuous outcomes."))
  }
  for (variable in names(result)) result[variable] <- .r4vn_superby_interaction_one(data, output, superby_name, variable)
  method <- if (continuous_outcome) "nested linear-model F tests" else if (isTRUE(output$effect_type %in% c("RR", "PR"))) "joint robust Wald tests" else "joint Wald tests"
  note <- paste0("Interaction p-values test each predictor by ", superby_name, " interaction using ", method, ". Covariates follow the multivariable model when the predictor is included in `multi`; otherwise they follow `adjusted`, and otherwise the interaction is unadjusted.")
  list(show = TRUE, p = result, note = note)
}

.r4vn_tab_superby <- function(data, vars, superby_name, user_call, caller_env,
                              interaction = TRUE, p_digit = 3, bold_p = TRUE,
                              p_bold = 0.05, template = "journal", append = NULL,
                              file = NULL, raw = FALSE, name = FALSE,
                              title = NULL, show = TRUE) {
  if (!superby_name %in% names(data)) stop(sprintf("`superby` variable `%s` was not found in `data`.", superby_name), call. = FALSE)
  if (superby_name %in% vars$variable) stop("The `superby` variable cannot also be included in `vars`.", call. = FALSE)
  super_vector <- data[[superby_name]]
  super_levels <- .r4vn_superby_levels(super_vector)
  if (length(super_levels) < 2L) stop("`superby` must contain at least two observed groups.", call. = FALSE)
  if (is.numeric(super_vector) && length(super_levels) > 20L && is.null(attr(super_vector, "labels", exact = TRUE))) {
    stop("`superby` appears continuous. Convert it to a factor or grouped variable before using it.", call. = FALSE)
  }

  run_child <- function(child_data) {
    child_call <- user_call
    child_call[[1L]] <- quote(tab)
    child_call$data <- quote(.r4vn_child_data)
    child_call$vars <- quote(.r4vn_child_vars)
    child_call$superby <- NULL
    child_call$interaction <- NULL
    child_call$append <- NULL
    child_call$file <- quote(.r4vn_child_file)
    child_call$title <- NULL
    child_call$show <- FALSE
    env <- new.env(parent = caller_env)
    env$.r4vn_child_data <- child_data
    env$.r4vn_child_vars <- vars
    env$.r4vn_child_file <- tempfile(pattern = "r4vn-superby-child-", fileext = ".html")
    eval(child_call, envir = env)
  }

  overall_output <- run_child(data)
  if (identical(overall_output$by, superby_name)) stop("`by` and `superby` must be different variables.", call. = FALSE)
  if ((!is.null(overall_output$adjusted) && superby_name %in% overall_output$adjusted$variable) ||
      (!is.null(overall_output$multi) && superby_name %in% overall_output$multi$variable)) {
    stop("Do not include the `superby` variable in `adjusted` or `multi`; it is handled automatically in interaction models.", call. = FALSE)
  }
  master_keys <- .r4vn_superby_row_keys(overall_output$rows)

  prepare_child_data <- function(child_data) {
    if (!identical(overall_output$outcome_type, "continuous") && !is.null(overall_output$by) && length(overall_output$by_levels)) {
      old_label <- attr(child_data[[overall_output$by]], "label", exact = TRUE)
      child_data[[overall_output$by]] <- factor(child_data[[overall_output$by]], levels = overall_output$by_levels)
      if (!is.null(old_label)) attr(child_data[[overall_output$by]], "label") <- old_label
    }
    categorical_meta <- rbind(overall_output$metadata, overall_output$adjusted, overall_output$multi)
    categorical_meta <- categorical_meta[!duplicated(categorical_meta$variable), , drop = FALSE]
    categorical_meta <- categorical_meta[categorical_meta$type == "categorical" & categorical_meta$variable %in% names(child_data), , drop = FALSE]
    if (nrow(categorical_meta)) {
      for (z in categorical_meta$variable) {
        old_label <- attr(child_data[[z]], "label", exact = TRUE)
        child_data[[z]] <- factor(child_data[[z]], levels = .r4vn_superby_levels(data[[z]]))
        if (!is.null(old_label)) attr(child_data[[z]], "label") <- old_label
      }
    }
    child_data
  }

  child_outputs <- vector("list", length(super_levels))
  child_errors <- character(length(super_levels))
  for (j in seq_along(super_levels)) {
    level <- super_levels[j]
    mask <- !is.na(super_vector) & as.character(super_vector) == as.character(level)
    child_data <- prepare_child_data(data[mask, , drop = FALSE])
    child_outputs[[j]] <- tryCatch(run_child(child_data), error = function(e) {
      child_errors[j] <<- conditionMessage(e)
      NULL
    })
  }

  cluster_labels <- c("Overall", as.character(super_levels))
  cluster_counts <- c(nrow(data), vapply(super_levels, function(level) sum(!is.na(super_vector) & as.character(super_vector) == as.character(level)), numeric(1)))
  all_outputs <- c(list(overall_output), child_outputs)
  master_n <- length(master_keys)

  align_block <- function(output) {
    if (is.null(output)) {
      values <- matrix("", nrow = master_n, ncol = ncol(overall_output$data) - 1L)
      headers_text <- names(overall_output$data)[-1L]
      headers_html <- overall_output$column_headers_html
      return(list(values = values, headers_text = headers_text, headers_html = headers_html, note = ""))
    }
    child_keys <- .r4vn_superby_row_keys(output$rows)
    index <- match(master_keys, child_keys)
    values <- matrix("", nrow = master_n, ncol = max(0L, ncol(output$data) - 1L))
    if (ncol(values) && any(!is.na(index))) values[!is.na(index), ] <- as.matrix(output$data[index[!is.na(index)], -1L, drop = FALSE])
    headers_text <- names(output$data)[-1L]
    headers_html <- output$column_headers_html
    if (length(headers_html) != length(headers_text)) headers_html <- .r4vn_superby_escape(headers_text)
    list(values = values, headers_text = headers_text, headers_html = headers_html, note = output$note_html)
  }
  blocks_aligned <- lapply(all_outputs, align_block)
  if (!sum(vapply(blocks_aligned, function(z) ncol(z$values), integer(1)))) stop("No result columns were produced for the superby table.", call. = FALSE)

  interaction_result <- .r4vn_superby_interactions(data, overall_output, superby_name, interaction)
  interaction_by_row <- rep(NA_real_, master_n)
  if (interaction_result$show) {
    for (i in seq_along(overall_output$rows)) {
      current <- overall_output$rows[[i]]
      if (isTRUE(current$variable_start) && current$variable %in% names(interaction_result$p)) interaction_by_row[i] <- interaction_result$p[current$variable]
    }
  }

  format_p_text <- function(p) {
    if (!length(p) || is.na(p) || !is.finite(p)) return("")
    limit <- 10^(-p_digit)
    text <- if (p < limit) paste0("&lt;", formatC(limit, format = "f", digits = p_digit)) else formatC(p, format = "f", digits = p_digit)
    if (isTRUE(bold_p) && p < p_bold) paste0("<strong>", text, "</strong>") else text
  }
  cell_html <- function(text, column_name) {
    text <- as.character(text)
    if (!length(text) || is.na(text) || !nzchar(text)) return("")
    escaped <- .r4vn_superby_escape(text)
    is_p_column <- grepl("(^|[ -])p($|[ -])|Test p", column_name, ignore.case = TRUE)
    if (is_p_column && isTRUE(bold_p)) {
      numeric_text <- sub("^<", "", trimws(text))
      numeric_text <- sub("[^0-9.].*$", "", numeric_text)
      p <- suppressWarnings(as.numeric(numeric_text))
      if (is.finite(p) && p < p_bold) escaped <- paste0("<strong>", escaped, "</strong>")
    }
    escaped
  }
  characteristic_html <- vapply(overall_output$rows, function(current) {
    if (current$type %in% c("categorical_header", "numeric_header", "numeric")) {
      paste0("<span class=\"variable-name\">", current$label, "</span>")
    } else {
      paste0("<span class=\"level-name", if (current$type == "missing") " missing-name" else "", "\">", .r4vn_superby_escape(current$item), "</span>")
    }
  }, character(1))

  top_headers <- character(length(blocks_aligned))
  second_headers <- character(length(blocks_aligned))
  for (j in seq_along(blocks_aligned)) {
    colspan <- ncol(blocks_aligned[[j]]$values)
    top_headers[j] <- paste0("<th class=\"superby-group\" colspan=\"", colspan, "\">", .r4vn_superby_escape(cluster_labels[j]), "<span class=\"header-n\">n = ", format(cluster_counts[j], big.mark = ",", scientific = FALSE, trim = TRUE), "</span></th>")
    second_headers[j] <- paste0("<th class=\"result-head\">", blocks_aligned[[j]]$headers_html, "</th>", collapse = "")
  }
  header_html <- paste0("<thead><tr class=\"superby-title\"><th rowspan=\"2\">Characteristic</th>", paste(top_headers, collapse = ""),
    if (interaction_result$show) "<th rowspan=\"2\">Interaction p</th>" else "", "</tr><tr>", paste(second_headers, collapse = ""), "</tr></thead>")

  body_html <- character(master_n)
  for (i in seq_len(master_n)) {
    cells <- character()
    for (j in seq_along(blocks_aligned)) {
      block <- blocks_aligned[[j]]
      if (ncol(block$values)) {
        for (k in seq_len(ncol(block$values))) cells <- c(cells, paste0("<td class=\"result\">", cell_html(block$values[i, k], block$headers_text[k]), "</td>"))
      }
    }
    interaction_cell <- if (interaction_result$show) paste0("<td class=\"interaction-p\">", format_p_text(interaction_by_row[i]), "</td>") else ""
    current <- overall_output$rows[[i]]
    row_class <- paste0("row-", gsub("_", "-", current$type, fixed = TRUE), if (isTRUE(current$variable_start)) " variable-start" else "")
    body_html[i] <- paste0("<tr class=\"", row_class, "\"><td>", characteristic_html[i], "</td>", paste(cells, collapse = ""), interaction_cell, "</tr>")
  }

  note_blocks <- character(length(blocks_aligned))
  for (j in seq_along(blocks_aligned)) {
    if (nzchar(blocks_aligned[[j]]$note)) note_blocks[j] <- paste0("<div class=\"superby-note\"><strong>", .r4vn_superby_escape(cluster_labels[j]), ":</strong>", blocks_aligned[[j]]$note, "</div>")
  }
  error_note <- ""
  if (any(nzchar(child_errors))) {
    messages <- paste0(.r4vn_superby_escape(as.character(super_levels[nzchar(child_errors)])), ": ", .r4vn_superby_escape(child_errors[nzchar(child_errors)]))
    error_note <- paste0("<div class=\"model-note\">Some subgroup blocks could not be estimated and were left blank: ", paste(messages, collapse = "; "), ".</div>")
  }
  missing_super <- sum(is.na(super_vector))
  super_note <- paste0("<div class=\"table-note\">Columns are shown for the full dataset and separately for each observed level of ", .r4vn_superby_escape(superby_name), ".", if (missing_super > 0L) paste0(" ", missing_super, " observation(s) with missing `superby` are included in Overall but excluded from subgroup blocks and interaction models.") else "", "</div>")
  interaction_note <- if (interaction_result$show) paste0("<div class=\"effect-note\">", .r4vn_superby_escape(interaction_result$note), "</div>") else if (nzchar(interaction_result$note)) paste0("<div class=\"effect-note\">", .r4vn_superby_escape(interaction_result$note), "</div>") else ""
  note_html <- paste0(super_note, paste(note_blocks, collapse = ""), interaction_note, error_note)
  title_html <- if (is.null(title) || !nzchar(as.character(title)[1L])) "" else paste0("<div class=\"table-title\">", .r4vn_superby_escape(as.character(title)[1L]), "</div>")

  css <- overall_output$css
  common_css <- paste0(overall_output$common_css, ".superby-title th,.superby-group{text-align:center}.interaction-p{text-align:right;white-space:nowrap}.superby-note{margin-top:8px;padding-top:4px;border-top:1px solid #ddd}.superby-note>.table-note,.superby-note>.test-note,.superby-note>.effect-note,.superby-note>.model-note{margin-left:12px}")
  table_block <- paste0("<section class=\"r4vn-table r4vn-superby-table\">", title_html, "<table>", header_html, "<tbody>", paste(body_html, collapse = ""), "</tbody></table>", note_html, "</section>")
  blocks <- table_block
  if (inherits(append, "r4vn_tab")) blocks <- c(append$blocks, table_block)
  document <- paste0("<!DOCTYPE html><html><head><meta charset=\"UTF-8\"><meta name=\"viewport\" content=\"width=device-width,initial-scale=1\"><style>", css, common_css, "</style></head><body><div class=\"table-wrapper\">", paste(blocks, collapse = "<div class=\"table-separator\"></div>"), "</div></body></html>")
  if (is.null(file)) file <- tempfile(pattern = "r4vn-superby-", fileext = ".html")
  if (!is.character(file) || length(file) != 1L || !nzchar(file)) stop("`file` must be a single valid file path.", call. = FALSE)
  if (is.character(append) && length(append) == 1L && file.exists(append)) {
    old <- paste(readLines(append, warn = FALSE, encoding = "UTF-8"), collapse = "\n")
    if (grepl("</body>", old, fixed = TRUE)) {
      document <- sub("</body>", paste0("<div class=\"table-separator\"></div>", table_block, "</body>"), old, fixed = TRUE)
      file <- append
    }
  }
  writeLines(enc2utf8(document), file, useBytes = TRUE)

  export_blocks <- vector("list", length(blocks_aligned))
  for (j in seq_along(blocks_aligned)) {
    block <- as.data.frame(blocks_aligned[[j]]$values, stringsAsFactors = FALSE, check.names = FALSE)
    names(block) <- paste0(cluster_labels[j], " | ", blocks_aligned[[j]]$headers_text)
    export_blocks[[j]] <- block
  }
  characteristic <- overall_output$data[[1L]]
  table_df <- data.frame(Characteristic = characteristic, stringsAsFactors = FALSE, check.names = FALSE)
  for (block in export_blocks) table_df <- cbind(table_df, block)
  if (interaction_result$show) table_df[["Interaction p"]] <- vapply(interaction_by_row, function(p) {
    if (!length(p) || is.na(p) || !is.finite(p)) return("")
    limit <- 10^(-p_digit)
    if (p < limit) paste0("<", formatC(limit, format = "f", digits = p_digit)) else formatC(p, format = "f", digits = p_digit)
  }, character(1))
  names(table_df) <- make.unique(names(table_df), sep = "_")

  output <- list(data = table_df, rows = overall_output$rows,
    raw = if (isTRUE(raw)) lapply(all_outputs, function(z) if (is.null(z)) NULL else z$raw) else NULL,
    metadata = vars, by = overall_output$by, by_levels = overall_output$by_levels,
    by_specification = overall_output$by_specification, outcome_type = overall_output$outcome_type,
    outcome_summary = overall_output$outcome_summary, effect_type = overall_output$effect_type,
    event = overall_output$event, adjusted = overall_output$adjusted,
    adjusted_all = overall_output$adjusted_all, multi = overall_output$multi,
    multi_model = overall_output$multi_model, multi_diagnostics = overall_output$multi_diagnostics,
    superby = superby_name, superby_levels = as.character(super_levels),
    subgroup_tables = stats::setNames(child_outputs, as.character(super_levels)),
    subgroup_models = stats::setNames(lapply(child_outputs, function(z) if (is.null(z)) NULL else z$multi_model), as.character(super_levels)),
    interaction_p = interaction_result$p, descriptive = overall_output$descriptive,
    html = document, table_html = table_block, blocks = blocks,
    file = normalizePath(file, winslash = "/", mustWork = TRUE), call = user_call)
  class(output) <- "r4vn_tab"
  if (show) {
    viewer <- getOption("viewer")
    if (is.function(viewer)) viewer(output$file) else utils::browseURL(output$file)
  }
  invisible(output)
}


#' Create Descriptive, Comparative, and Regression Tables
#'
#' Creates publication-style tables for descriptive analysis, group comparisons,
#' binary-outcome regression, and continuous-outcome linear regression.
#'
#' @usage
#' tab(
#'   ...,
#'   data = NULL,
#'   vars = NULL,
#'   by = NULL,
#'   superby = NULL,
#'   digit = 1,
#'   p_digit = 3,
#'   effect_digit = 2,
#'   missing = "ifany",
#'   row = FALSE,
#'   col = TRUE,
#'   cell = FALSE,
#'   overall = "first",
#'   descriptive = TRUE,
#'   rvrow = NULL,
#'   rvcol = FALSE,
#'   test = TRUE,
#'   pvalue = TRUE,
#'   bold_p = TRUE,
#'   p_bold = 0.05,
#'   test_note = TRUE,
#'   interaction = TRUE,
#'   or = FALSE,
#'   rr = FALSE,
#'   pr = FALSE,
#'   event = NULL,
#'   adjusted = NULL,
#'   multi = NULL,
#'   effect_ref = NULL,
#'   template = c("journal", "clean", "minimal"),
#'   append = NULL,
#'   file = NULL,
#'   raw = FALSE,
#'   name = FALSE,
#'   title = NULL,
#'   show = TRUE,
#'   mode = c("auto", "console", "table")
#' )
#'
#' @param ... In console mode, one row variable and optionally one column
#'   variable, followed by console options such as `exp`, `chi`, and `fisher`.
#'   In publication mode, legacy positional `data`, `vars`, and `by` arguments
#'   are also accepted.
#' @param data Optional data frame. When omitted or \code{NULL}, the active
#'   data frame set by \code{usedf()} or \code{opendata(..., active = TRUE)}
#'   is used.
#' @param vars A variable specification created by \code{vars()}.
#' @param by Optional grouping or outcome variable supplied without quotation
#'   marks. Leave it empty for an overall descriptive table. Use a regular
#'   variable name for a categorical grouping/outcome variable, \code{c.outcome}
#'   for a continuous outcome summarized by mean (SD), or \code{q.outcome} for
#'   a continuous outcome summarized by median (IQR). R4VN also accepts the
#'   unified hierarchical form \code{by = vars(province, sex, outcome)}:
#'   \code{province} and \code{sex} are nested superby strata, in that order,
#'   and \code{outcome} is the innermost grouping/outcome variable.
#' @param superby Backward-compatible single stratification variable. It may be
#'   combined with hierarchical \code{by = vars(...)} and then becomes the
#'   outermost stratum. New code should normally prefer the unified \code{by}
#'   convention.
#' @param digit Number of decimal places for descriptive statistics.
#' @param p_digit Number of decimal places for p-values.
#' @param effect_digit Number of decimal places for OR, RR, PR, or linear
#'   regression coefficients.
#' @param missing Missing-value display for categorical variables:
#'   \code{"no"}, \code{"ifany"}, or \code{"always"}.
#' @param row Logical. Calculate row percentages when \code{by} is categorical.
#'   When \code{row = TRUE}, \code{col} and \code{cell} are automatically
#'   set to \code{FALSE}.
#' @param col Logical. Calculate column percentages when \code{by} is categorical.
#'   This is the default percentage mode.
#' @param cell Logical. Calculate percentages using the complete table total.
#'   When \code{cell = TRUE} and \code{row = FALSE}, \code{row} and
#'   \code{col} are automatically set to \code{FALSE}. Thus users normally
#'   need to specify only \code{row = TRUE}, \code{cell = TRUE}, or neither
#'   for the default column percentages.
#' @param overall Position of the overall column: \code{"none"},
#'   \code{"first"}, or \code{"last"}. Logical values are accepted for
#'   backward compatibility.
#' @param descriptive Logical. Display descriptive-statistics columns.
#' @param rvrow Categorical variables whose displayed level order should be
#'   reversed. Accepts \code{TRUE}, \code{vars(...)}, \code{c(...)}, a
#'   single variable name, or a character vector. This does not change model
#'   reference categories.
#' @param rvcol Logical. Reverse displayed levels of a categorical \code{by}
#'   variable.
#' @param test Logical. Display traditional omnibus-test p-values.
#' @param pvalue Logical. Display separate p-value columns for model coefficients.
#' @param bold_p Logical. Bold p-values smaller than \code{p_bold}.
#' @param p_bold Significance threshold used when \code{bold_p = TRUE}.
#' @param test_note Logical. Add superscript letters and footnotes identifying
#'   omnibus tests.
#' @param interaction Logical. When \code{superby} is supplied, add one final
#'   interaction p-value column. Interaction tests use predictor-by-superby
#'   terms and follow \code{multi}, then \code{adjusted}, then crude models.
#' @param or Logical. Calculate odds ratios using logistic regression.
#' @param rr Logical. Calculate risk ratios using modified Poisson regression
#'   with robust variance.
#' @param pr Logical. Calculate prevalence ratios using modified Poisson
#'   regression with robust variance.
#' @param event Event level of a binary outcome. The last observed level is used
#'   when omitted.
#' @param adjusted Variables included as adjustment covariates in separate models
#'   for each focal predictor. Prefer \code{vars(c.age, b2.sex, q.bmi)} so
#'   variable types and reference levels remain explicit. Also accepts
#'   \code{c(...)}, a character vector, \code{TRUE}, or \code{"ALL"}.
#' @param multi Variables included together in one final multivariable model.
#'   Prefer \code{vars(c.age, b2.sex, c.bmi)}. \code{TRUE} or
#'   \code{"ALL"} includes every variable listed in \code{vars}.
#' @param effect_ref Optional backward-compatible reference categories. The
#'   \code{b2.}, \code{b3.}, and related prefixes take precedence.
#' @param template HTML style: \code{"journal"}, \code{"clean"}, or
#'   \code{"minimal"}.
#' @param append Optional previous \code{r4vn_tab} object or existing HTML path.
#' @param file Optional output HTML path. A temporary file is created when omitted.
#' @param raw Logical. Retain unformatted results in the returned object.
#' @param name Logical. Display original variable names beside variable labels.
#' @param title Optional table title.
#' @param show Logical. Display the HTML table in the RStudio Viewer or browser.
#' @param mode Dispatch mode. `"auto"` selects publication mode when a
#'   `vars()` specification is supplied and otherwise selects console mode.
#'   Use `"console"` or `"table"` to force a mode.
#'
#' @section Common call patterns:
#' \preformatted{
#' tab(data, vars = vars(...))
#' tab(data, vars = vars(...), by = group)
#' tab(data, vars = vars(...), by = outcome, or = TRUE)
#' tab(data, vars = vars(...), by = c.outcome)
#' tab(data, vars = vars(...), by = q.outcome)
#' tab(data, vars = vars(...), by = outcome, superby = subgroup, or = TRUE)
#' }
#'
#' @details
#' Prefixes used inside \code{vars()} determine descriptive summaries and
#' categorical reference levels:
#' \itemize{
#'   \item no prefix: automatic typing; numeric/integer variables use mean and
#'     standard deviation, while factor/character/logical variables are
#'     categorical with the first observed level as reference;
#'   \item \code{b1.}, \code{b2.}, \code{b3.}, ...: force a categorical variable with the
#'     corresponding observed level as reference;
#'   \item \code{c.}: mean and standard deviation;
#'   \item \code{q.}: median and interquartile range;
#'   \item \code{f.}: mean, median, and range.
#' }
#'
#' For grouped categorical tables, column percentages are the default. Setting
#' \code{row = TRUE} automatically turns \code{col} and \code{cell} off;
#' setting \code{cell = TRUE} automatically turns \code{row} and \code{col}
#' off. Users therefore do not need to manually disable \code{col = TRUE}.
#'
#' With a categorical \code{by} variable, categorical predictors are tested
#' using Pearson's chi-squared test or Fisher's exact test. Variables declared
#' with \code{c.} use a t-test or one-way ANOVA; variables declared with
#' \code{q.} or \code{f.} use the Wilcoxon rank-sum or Kruskal-Wallis test.
#'
#' Binary outcomes can be analyzed with OR, RR, or PR. OR uses logistic
#' regression. RR and PR use modified Poisson regression with robust variance.
#'
#' With \code{by = c.outcome}, the continuous outcome is summarized by mean
#' (SD), categorical predictors use t-tests/ANOVA, and numeric predictors use
#' Pearson correlation tests. With \code{by = q.outcome}, the outcome is
#' summarized by median (IQR), categorical predictors use
#' Wilcoxon/Kruskal-Wallis tests, and numeric predictors use Spearman tests.
#' Both modes report unstandardized beta coefficients from linear regression.
#'
#' \code{adjusted} and \code{multi} have different roles. \code{adjusted}
#' fits a separate adjusted model for each focal predictor. \code{multi} fits
#' one final model containing all specified variables.
#'
#' When \code{superby} is supplied, \code{tab()} first calculates the complete
#' dataset and then repeats the same analysis independently within every level
#' of \code{superby}. The resulting blocks are combined side by side. When
#' possible, one final interaction p-value column tests whether each predictor
#' effect differs across the levels of \code{superby}.
#'
#' @return Invisibly returns an object of class \code{r4vn_tab}. Important
#'   components include \code{data}, \code{file}, \code{html},
#'   \code{table_html}, \code{rows}, \code{multi_model}, and
#'   \code{multi_diagnostics}.
#'
#' @seealso \code{\link{vars}}, \code{\link{tabmulti}}, and
#'   \code{\link{tabexport}}.
#' @aliases tabconti
#' @family R4VN tables
#'
#' @examples
#' set.seed(2026)
#' n <- 180
#' dat <- data.frame(
#'   age = round(rnorm(n, 45, 12)),
#'   sex = factor(sample(c("Female", "Male"), n, TRUE)),
#'   bmi = round(rnorm(n, 23, 3), 1),
#'   smoking = factor(sample(c("No", "Yes"), n, TRUE,
#'                           prob = c(0.70, 0.30))),
#'   education = factor(sample(c("Primary", "Secondary", "College"),
#'                             n, TRUE))
#' )
#' dat$sbp <- round(80 + 0.75 * dat$age + 1.1 * dat$bmi +
#'                  5 * (dat$sex == "Male") +
#'                  4 * (dat$smoking == "Yes") + rnorm(n, 0, 12), 1)
#' lp <- -3.2 + 0.045 * dat$age + 0.10 * (dat$bmi - 23) +
#'       0.45 * (dat$sex == "Male") + 0.65 * (dat$smoking == "Yes")
#' dat$hypertension <- factor(
#'   rbinom(n, 1, plogis(lp)),
#'   levels = c(0, 1), labels = c("No", "Yes")
#' )
#'
#' tb0 <- tab(dat, vars = vars(c.age, b2.sex, q.bmi, b2.smoking, education),
#'            show = FALSE)
#' head(tb0$data)
#'
#' tb1 <- tab(dat, vars = vars(c.age, b2.sex, c.bmi, b2.smoking, b2.education),
#'            by = hypertension, or = TRUE, event = "Yes",
#'            multi = vars(c.age, b2.sex, c.bmi, b2.smoking), show = FALSE)
#'
#' tb2 <- tab(dat, vars = vars(c.age, b2.sex, c.bmi, b2.smoking),
#'            by = c.sbp, multi = vars(c.age, b2.sex, c.bmi, b2.smoking),
#'            show = FALSE)
#'
#' tb3 <- tab(dat, vars = vars(c.age, c.bmi, b2.smoking, b2.education),
#'            by = hypertension, superby = sex, overall = "none",
#'            or = TRUE, event = "Yes",
#'            multi = vars(c.age, c.bmi, b2.smoking), show = FALSE)
#'
#' # Extended usage examples
#' \donttest{
#' set.seed(2026)
#' n <- 300
#' d <- data.frame(
#'   sex = factor(sample(c("Female", "Male"), n, TRUE)),
#'   age = rnorm(n, 45, 12),
#'   bmi = rnorm(n, 23, 3),
#'   smoking = factor(sample(c("No", "Yes"), n, TRUE)),
#'   region = factor(sample(c("Urban", "Rural"), n, TRUE)),
#'   outcome = factor(rbinom(n, 1, .3), levels = 0:1, labels = c("No", "Yes")),
#'   sbp = rnorm(n, 125, 18)
#' )
#'
#' # Overall descriptive table. Numeric variables without a prefix are
#' # automatically summarized with mean (SD); factors remain categorical.
#' t1_auto <- tab(d, vars = vars(age, sex, bmi, smoking), show = FALSE)
#'
#' # Explicit q. remains available when median (IQR) is preferred.
#' t1 <- tab(d, vars = vars(sex, age, q.bmi, smoking), show = FALSE)
#'
#' # Compare groups, show overall first, tests, and missing values when present
#' t2 <- tab(d, vars = vars(sex, c.age, q.bmi, smoking), by = outcome,
#'           overall = "first", test = TRUE, missing = "ifany", show = FALSE)
#'
#' # Row, column, or cell percentages for categorical variables
#' tab(d, vars = vars(sex, smoking), by = outcome, row = TRUE, show = FALSE)
#' tab(d, vars = vars(sex, smoking), by = outcome, show = FALSE)
#' tab(d, vars = vars(sex, smoking), by = outcome, cell = TRUE, show = FALSE)
#'
#' # Reverse selected row levels or the by-variable columns
#' tab(d, vars = vars(sex, smoking), by = outcome,
#'     rvrow = vars(smoking), rvcol = TRUE, show = FALSE)
#'
#' # Crude odds ratios for a binary outcome
#' tab(d, vars = vars(sex, c.age, smoking), by = outcome,
#'     or = TRUE, event = "Yes", show = FALSE)
#'
#' # Risk ratios or prevalence ratios using modified Poisson models
#' tab(d, vars = vars(sex, c.age, smoking), by = outcome,
#'     rr = TRUE, event = "Yes", show = FALSE)
#' tab(d, vars = vars(sex, c.age, smoking), by = outcome,
#'     pr = TRUE, event = "Yes", show = FALSE)
#'
#' # Separate adjusted models for every focal predictor
#' tab(d, vars = vars(sex, c.age, smoking), by = outcome,
#'     or = TRUE, adjusted = vars(age, sex), event = "Yes", show = FALSE)
#'
#' # One final multivariable model; effects are placed beside their variables
#' tab(d, vars = vars(sex, c.age, smoking, q.bmi), by = outcome,
#'     or = TRUE, multi = vars(sex, age, smoking), event = "Yes", show = FALSE)
#'
#' # Hide descriptive columns and show only model results
#' tab(d, vars = vars(sex, c.age, smoking), by = outcome,
#'     descriptive = FALSE, or = TRUE, multi = TRUE,
#'     event = "Yes", show = FALSE)
#'
#' # Continuous outcome: c. gives parametric methods and beta coefficients
#' tab(d, vars = vars(sex, c.age, smoking, q.bmi), by = c.sbp,
#'     adjusted = vars(age, sex), multi = vars(age, sex, bmi), show = FALSE)
#'
#' # Continuous outcome: q. gives rank-based descriptive comparisons
#' tab(d, vars = vars(sex, c.age, smoking, q.bmi), by = q.sbp,
#'     test = TRUE, show = FALSE)
#'
#' # Supergroup columns plus interaction
#' tab(d, vars = vars(sex, c.age, smoking), by = outcome, superby = region,
#'     interaction = TRUE, overall = "first", show = FALSE)
#'
#' # Templates, titles, raw numerical output, and named variables
#' tab(d, vars = vars(sex, c.age, smoking), by = outcome,
#'     template = "minimal", title = "Participant characteristics",
#'     raw = TRUE, name = TRUE, show = FALSE)
#' }
#' @export
tab <- function(data = NULL, vars = NULL, by = NULL, superby = NULL, digit = 1, p_digit = 3, effect_digit = 2,
                missing = "ifany", row = FALSE, col = TRUE, cell = FALSE,
                overall = "first", descriptive = TRUE, rvrow = NULL, rvcol = FALSE, test = TRUE,
                pvalue = TRUE, bold_p = TRUE, p_bold = 0.05, test_note = TRUE, interaction = TRUE,
                or = FALSE, rr = FALSE, pr = FALSE, event = NULL,
                adjusted = NULL, multi = NULL, effect_ref = NULL,
                template = c("journal", "clean", "minimal"), append = NULL,
                file = NULL, raw = FALSE, name = FALSE, title = NULL,
                show = TRUE) {
  # Use active data when `data` is omitted.
  if (is.null(data)) data <- .r4vn_get_active()

  # Validate general arguments.
  if (!is.data.frame(data)) stop("`data` must be a data frame.", call. = FALSE)
  if (!inherits(vars, "r4vn_vars")) stop("`vars` must be created using `vars()`.", call. = FALSE)
  validate_integer <- function(x, arg) {
    if (!is.numeric(x) || length(x) != 1L || is.na(x) || x < 0 || x != floor(x)) {
      stop(sprintf("`%s` must be a single non-negative integer.", arg), call. = FALSE)
    }
  }
  validate_flag <- function(x, arg) {
    if (!is.logical(x) || length(x) != 1L || is.na(x)) stop(sprintf("`%s` must be TRUE or FALSE.", arg), call. = FALSE)
  }
  validate_integer(digit, "digit")
  validate_integer(p_digit, "p_digit")
  validate_integer(effect_digit, "effect_digit")
  missing <- match.arg(missing, c("no", "ifany", "always"))
  template <- match.arg(template)
  for (arg in c("descriptive", "rvcol", "test", "pvalue", "bold_p", "test_note", "interaction", "or", "rr", "pr", "raw", "name", "show")) {
    validate_flag(get(arg), arg)
  }
  if (!is.numeric(p_bold) || length(p_bold) != 1L || is.na(p_bold) || p_bold < 0 || p_bold > 1) {
    stop("`p_bold` must be between 0 and 1.", call. = FALSE)
  }
  if (is.logical(overall) && length(overall) == 1L && !is.na(overall)) overall <- if (overall) "first" else "none"
  overall <- match.arg(overall, c("none", "first", "last"))
  effect_flags <- c(OR = or, RR = rr, PR = pr)
  if (sum(effect_flags) > 1L) stop("Only one of `or`, `rr`, or `pr` may be TRUE.", call. = FALSE)
  effect_type <- if (any(effect_flags)) names(effect_flags)[which(effect_flags)] else NULL

  # superby = NULL preserves the original engine. A non-NULL superby creates
  # one complete-data block plus one independently calculated block per group.
  superby_expression <- substitute(superby)
  if (!identical(superby_expression, quote(NULL))) {
    if (!is.symbol(superby_expression)) stop("`superby` must be a single variable name.", call. = FALSE)

    # Resolve deferred vars() selectors only after the actual data are known.
    # For broad selectors such as vars(.) or wildcard selectors, structural
    # variables used as by/superby are removed automatically.
    selector_mode <- any(
      vars$variable == "." |
        grepl("*", vars$variable, fixed = TRUE)
    )
    vars <- .r4vn_resolve_vars_input(
      vars,
      data = data,
      arg = "vars",
      default_type = "auto",
      strict = TRUE
    )
    if (isTRUE(selector_mode)) {
      structural_variables <- as.character(superby_expression)
      by_for_selector <- substitute(by)
      if (!identical(by_for_selector, quote(NULL)) && is.symbol(by_for_selector)) {
        structural_variables <- c(
          structural_variables,
          sub("^[cq]\\.", "", as.character(by_for_selector))
        )
      }
      vars <- vars[
        !vars$variable %in% unique(structural_variables),
        ,
        drop = FALSE
      ]
      rownames(vars) <- NULL
      class(vars) <- c("r4vn_vars", "data.frame")
      if (!nrow(vars)) {
        stop("No predictor variables remain after resolving `vars()`.", call. = FALSE)
      }
    }

    return(.r4vn_tab_superby(
      data = data, vars = vars, superby_name = as.character(superby_expression),
      user_call = match.call(), caller_env = parent.frame(), interaction = interaction,
      p_digit = p_digit, bold_p = bold_p, p_bold = p_bold, template = template,
      append = append, file = file, raw = raw, name = name, title = title, show = show
    ))
  }

  # Route continuous outcomes to the internal linear-regression engine.
  # Users still call only tab(), for example by = c.sbp or by = q.sbp.
  by_expression <- substitute(by)
  if (!identical(by_expression, quote(NULL)) && is.symbol(by_expression)) {
    by_specification <- as.character(by_expression)
    if (grepl("^[cq]\\.", by_specification)) {
      if (!exists(".tab_continuous", mode = "function", inherits = TRUE)) {
        stop("The continuous-outcome engine was not found. Add `tabconti.R` to the package R/ folder or source it before using `by = c.outcome`/`q.outcome`.", call. = FALSE)
      }

      selector_mode <- any(
        vars$variable == "." |
          grepl("*", vars$variable, fixed = TRUE)
      )
      vars <- .r4vn_resolve_vars_input(
        vars,
        data = data,
        arg = "vars",
        default_type = "auto",
        strict = TRUE
      )
      if (isTRUE(selector_mode)) {
        outcome_name_for_selector <- sub("^[cq]\\.", "", by_specification)
        vars <- vars[
          vars$variable != outcome_name_for_selector,
          ,
          drop = FALSE
        ]
        rownames(vars) <- NULL
        class(vars) <- c("r4vn_vars", "data.frame")
        if (!nrow(vars)) {
          stop("No predictor variables remain after resolving `vars()`.", call. = FALSE)
        }
      }

      return(.tab_continuous(
        data = data, vars = vars, outcome_spec = by_specification,
        digit = digit, p_digit = p_digit, effect_digit = effect_digit,
        missing = missing, overall = overall, descriptive = descriptive,
        rvrow_expr = substitute(rvrow), test = test, pvalue = pvalue,
        bold_p = bold_p, p_bold = p_bold, test_note = test_note,
        or = or, rr = rr, pr = pr, event = event,
        adjusted_expr = substitute(adjusted), multi_expr = substitute(multi),
        effect_ref = effect_ref, template = template, append = append,
        file = file, raw = raw, name = name, title = title, show = show,
        caller_env = parent.frame(), user_call = match.call()
      ))
    }
  }

  # Capture the grouping/outcome variable.
  by_expression <- substitute(by)
  has_by <- !identical(by_expression, quote(NULL))
  by_name <- if (has_by) {
    if (!is.symbol(by_expression)) stop("`by` must be a single variable name.", call. = FALSE)
    as.character(by_expression)
  } else NULL
  if (!has_by && !is.null(effect_type)) stop("OR, RR, and PR require a binary `by` outcome.", call. = FALSE)

  # Resolve vars(.), wildcard selectors, and exclusions against the actual data.
  # When a broad selector is used with `by`, the outcome/grouping variable is
  # removed automatically so it is not analysed as its own predictor.
  selector_mode <- any(
    vars$variable == "." |
      grepl("*", vars$variable, fixed = TRUE)
  )
  vars <- .r4vn_resolve_vars_input(
    vars,
    data = data,
    arg = "vars",
    default_type = "auto",
    strict = TRUE
  )
  if (isTRUE(selector_mode) && has_by && by_name %in% vars$variable) {
    vars <- vars[
      vars$variable != by_name,
      ,
      drop = FALSE
    ]
    rownames(vars) <- NULL
    class(vars) <- c("r4vn_vars", "data.frame")
    if (!nrow(vars)) {
      stop("No predictor variables remain after resolving `vars()`.", call. = FALSE)
    }
  }

  # Normalize grouped percentage options.
  #
  # `col = TRUE` is the publication-table default. Therefore a user should be
  # able to request row percentages simply with `row = TRUE`, without also
  # having to write `col = FALSE`. Likewise, `cell = TRUE` automatically
  # switches off the default column percentages.
  #
  # Precedence is intentionally simple:
  #   row = TRUE  -> row percentages;  col = FALSE; cell = FALSE
  #   cell = TRUE -> cell percentages; row = FALSE; col = FALSE
  #   otherwise   -> column percentages when col = TRUE
  #
  # If all three are FALSE for a grouped table, keep the explicit error because
  # no percentage denominator has been requested.
  if (has_by) {
    if (isTRUE(row)) {
      col <- FALSE
      cell <- FALSE
    } else if (isTRUE(cell)) {
      row <- FALSE
      col <- FALSE
    } else if (isTRUE(col)) {
      row <- FALSE
      cell <- FALSE
    }

    if (sum(c(row = isTRUE(row), col = isTRUE(col), cell = isTRUE(cell))) != 1L) {
      stop(
        "When `by` is supplied, use one percentage mode: `row = TRUE`, `cell = TRUE`, or the default `col = TRUE`.",
        call. = FALSE
      )
    }
  }

  # Parse rvrow while preserving unevaluated variable names.
  parse_simple_variable_list <- function(expression, available = NULL) {
    if (identical(expression, quote(NULL))) return(character())
    if (is.logical(expression) && length(expression) == 1L) {
      return(if (isTRUE(expression) && !is.null(available)) available else character())
    }
    if (is.character(expression)) return(expression)
    if (is.symbol(expression)) {
      object_name <- as.character(expression)
      value <- tryCatch(get(object_name, envir = parent.frame(2L)), error = function(e) NULL)
      if (is.logical(value) && length(value) == 1L) return(if (isTRUE(value) && !is.null(available)) available else character())
      if (inherits(value, "r4vn_vars")) return(value$variable)
      if (is.character(value)) return(value)
      return(object_name)
    }
    if (is.call(expression) && as.character(expression[[1L]]) %in% c("c", "vars")) {
      items <- as.list(expression)[-1L]
      if (all(vapply(items, is.symbol, logical(1)))) return(vapply(items, as.character, character(1)))
      value <- tryCatch(eval(expression, envir = parent.frame(2L)), error = function(e) NULL)
      if (inherits(value, "r4vn_vars")) return(value$variable)
      if (is.character(value)) return(value)
    }
    value <- tryCatch(eval(expression, envir = parent.frame(2L)), error = function(e) NULL)
    if (inherits(value, "r4vn_vars")) return(value$variable)
    if (is.character(value)) return(value)
    stop("The variable list could not be parsed.", call. = FALSE)
  }

  # Infer model metadata when adjusted variables are supplied without vars().
  infer_adjustment_metadata <- function(variable_names) {
    variable_names <- unique(as.character(variable_names))
    if (!length(variable_names)) {
      output <- data.frame(variable = character(), type = character(), specification = character(), reference_index = integer(), stringsAsFactors = FALSE)
      class(output) <- c("r4vn_vars", "data.frame")
      return(output)
    }
    output <- vector("list", length(variable_names))
    for (i in seq_along(variable_names)) {
      variable <- variable_names[i]
      if (!variable %in% names(data)) stop(sprintf("Adjustment variable `%s` was not found in `data`.", variable), call. = FALSE)
      table_index <- match(variable, vars$variable)
      if (!is.na(table_index)) {
        output[[i]] <- vars[table_index, , drop = FALSE]
      } else {
        x <- data[[variable]]
        categorical <- is.factor(x) || is.character(x) || is.logical(x)
        output[[i]] <- data.frame(
          variable = variable,
          type = if (categorical) "categorical" else "mean",
          specification = variable,
          reference_index = if (categorical) 1L else NA_integer_,
          stringsAsFactors = FALSE
        )
      }
    }
    output <- do.call(rbind, output)
    rownames(output) <- NULL
    class(output) <- c("r4vn_vars", "data.frame")
    output
  }

  # Parse adjusted. vars() preserves c./q./f. types and bN. references.
  parse_adjusted_metadata <- function(expression) {
    empty <- infer_adjustment_metadata(character())
    if (identical(expression, quote(NULL))) return(list(all = FALSE, meta = empty))
    if (is.logical(expression) && length(expression) == 1L) {
      return(list(all = isTRUE(expression), meta = if (isTRUE(expression)) vars else empty))
    }
    if (is.character(expression)) {
      if (length(expression) == 1L && toupper(expression) == "ALL") return(list(all = TRUE, meta = vars))
      return(list(all = FALSE, meta = infer_adjustment_metadata(expression)))
    }
    if (is.symbol(expression)) {
      object_name <- as.character(expression)
      if (toupper(object_name) == "ALL") return(list(all = TRUE, meta = vars))
      value <- tryCatch(get(object_name, envir = parent.frame(2L)), error = function(e) NULL)
      if (inherits(value, "r4vn_vars")) return(list(all = FALSE, meta = value))
      if (is.logical(value) && length(value) == 1L && isTRUE(value)) return(list(all = TRUE, meta = vars))
      if (is.character(value)) {
        if (length(value) == 1L && toupper(value) == "ALL") return(list(all = TRUE, meta = vars))
        return(list(all = FALSE, meta = infer_adjustment_metadata(value)))
      }
      return(list(all = FALSE, meta = infer_adjustment_metadata(object_name)))
    }
    if (is.call(expression) && identical(as.character(expression[[1L]]), "vars")) {
      value <- eval(expression, envir = parent.frame(2L))
      if (!inherits(value, "r4vn_vars")) stop("`adjusted = vars(...)` did not create a valid variable specification.", call. = FALSE)
      return(list(all = FALSE, meta = value))
    }
    if (is.call(expression) && identical(as.character(expression[[1L]]), "c")) {
      items <- as.list(expression)[-1L]
      if (all(vapply(items, is.symbol, logical(1)))) {
        return(list(all = FALSE, meta = infer_adjustment_metadata(vapply(items, as.character, character(1)))))
      }
      value <- eval(expression, envir = parent.frame(2L))
      if (is.character(value)) return(list(all = FALSE, meta = infer_adjustment_metadata(value)))
    }
    value <- tryCatch(eval(expression, envir = parent.frame(2L)), error = function(e) NULL)
    if (inherits(value, "r4vn_vars")) return(list(all = FALSE, meta = value))
    if (is.character(value)) return(list(all = FALSE, meta = infer_adjustment_metadata(value)))
    stop("`adjusted` must be NULL, TRUE, 'ALL', vars(...), c(...), or a character vector.", call. = FALSE)
  }

  categorical_variables <- vars$variable[vars$type == "categorical"]
  reverse_rows <- parse_simple_variable_list(substitute(rvrow), categorical_variables)
  invalid_reverse_rows <- setdiff(reverse_rows, categorical_variables)
  if (length(invalid_reverse_rows)) stop(sprintf("Invalid `rvrow` variables: %s.", paste(invalid_reverse_rows, collapse = ", ")), call. = FALSE)

  resolve_model_meta <- function(meta, arg) {
    if (is.null(meta) || !nrow(meta) || !inherits(meta, "r4vn_vars")) return(meta)
    needs_resolution <- any(
      meta$type == "default" | meta$variable == "." | grepl("*", meta$variable, fixed = TRUE)
    )
    if (!isTRUE(needs_resolution)) return(meta)
    .r4vn_resolve_vars_input(meta, data = data, arg = arg, default_type = "auto", strict = TRUE)
  }

  adjustment <- parse_adjusted_metadata(substitute(adjusted))
  adjust_all <- isTRUE(adjustment$all)
  adjusted_meta <- resolve_model_meta(adjustment$meta, "adjusted")
  if (anyDuplicated(adjusted_meta$variable)) adjusted_meta <- adjusted_meta[!duplicated(adjusted_meta$variable), , drop = FALSE]
  absent_adjustment <- setdiff(adjusted_meta$variable, names(data))
  if (length(absent_adjustment)) stop(sprintf("Adjustment variables not found in `data`: %s.", paste(absent_adjustment, collapse = ", ")), call. = FALSE)
  if (!is.null(by_name) && length(by_name) == 1L && by_name %in% adjusted_meta$variable) {
    stop("The outcome variable cannot be included in `adjusted`.", call. = FALSE)
  }
  has_adjusted <- !is.null(effect_type) && nrow(adjusted_meta) > 0L

  # Parse the variables used in one final multivariable model.
  multi_specification <- parse_adjusted_metadata(substitute(multi))
  multi_all <- isTRUE(multi_specification$all)
  multi_meta <- resolve_model_meta(multi_specification$meta, "multi")
  if (multi_all) multi_meta <- vars
  if (anyDuplicated(multi_meta$variable)) multi_meta <- multi_meta[!duplicated(multi_meta$variable), , drop = FALSE]
  if (nrow(multi_meta)) {
    if (is.null(effect_type)) stop("`multi` requires one of `or`, `rr`, or `pr`.", call. = FALSE)
    if (!has_by) stop("`multi` requires a binary `by` outcome.", call. = FALSE)
    absent_multi <- setdiff(multi_meta$variable, names(data))
    if (length(absent_multi)) stop(sprintf("Variables in `multi` were not found in `data`: %s.", paste(absent_multi, collapse = ", ")), call. = FALSE)
    outside_table <- setdiff(multi_meta$variable, vars$variable)
    if (length(outside_table)) stop(sprintf("Variables in `multi` must also appear in `vars`: %s.", paste(outside_table, collapse = ", ")), call. = FALSE)
    if (!is.null(by_name) && by_name %in% multi_meta$variable) stop("The outcome variable cannot be included in `multi`.", call. = FALSE)
  }
  has_multi <- !is.null(effect_type) && nrow(multi_meta) > 0L
  # Check requested variables.
  duplicated_variables <- unique(vars$variable[duplicated(vars$variable)])
  if (length(duplicated_variables)) stop(sprintf("Duplicated variable specification: %s.", paste(duplicated_variables, collapse = ", ")), call. = FALSE)
  absent_variables <- setdiff(unique(c(vars$variable, by_name)), names(data))
  if (length(absent_variables)) stop(sprintf("Variables not found in `data`: %s.", paste(absent_variables, collapse = ", ")), call. = FALSE)

  # Formatting helpers.
  escape_html <- function(x) {
    x <- as.character(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)
    gsub("'", "&#39;", x, fixed = TRUE)
  }
  format_number <- function(x, digits = digit) {
    if (!length(x) || is.na(x) || !is.finite(x)) return("")
    formatC(x, format = "f", digits = digits, big.mark = ",")
  }
  format_count <- function(x) {
    if (!length(x) || is.na(x)) return("")
    format(x, big.mark = ",", scientific = FALSE, trim = TRUE)
  }
  format_p <- function(x) {
    if (!length(x) || is.na(x) || !is.finite(x)) return("")
    limit <- 10^(-p_digit)
    text <- if (x < limit) paste0("&lt;", formatC(limit, format = "f", digits = p_digit)) else formatC(x, format = "f", digits = p_digit)
    if (isTRUE(bold_p) && x < p_bold) paste0("<strong>", text, "</strong>") else text
  }
  get_label <- function(x, variable) {
    label <- attr(x, "label", exact = TRUE)
    if (is.null(label) || !length(label) || is.na(label[1L]) || !nzchar(as.character(label[1L]))) variable else as.character(label[1L])
  }
  display_label <- function(label, variable) {
    if (!isTRUE(name)) return(escape_html(label))
    paste0(escape_html(label), " <span class=\"variable-code\">[", escape_html(variable), "]</span>")
  }
  get_levels <- function(x, reverse = FALSE) {
    observed <- x[!is.na(x)]
    values <- if (is.factor(x)) levels(x) else if (is.logical(x)) c(FALSE, TRUE) else {
      z <- unique(observed)
      if (is.numeric(z)) sort(z) else z
    }
    values <- values[values %in% observed]
    if (reverse) rev(values) else values
  }
  compact_value <- function(z, kind) {
    if (is.null(z) || !length(z)) return("")
    if (kind == "categorical") return(if (nzchar(z["first"])) paste0(z["first"], " (", z["second"], ")") else "")
    if (kind %in% c("mean", "median")) return(if (nzchar(z["first"])) paste0(z["first"], " (", z["second"], ")") else "")
    if (kind == "range") return(if (nzchar(z["first"])) paste0(z["first"], " - ", z["second"]) else "")
    ""
  }

  # Prepare the grouping/outcome variable.
  if (has_by) {
    by_vector <- data[[by_name]]
    by_label <- get_label(by_vector, by_name)
    original_by_levels <- get_levels(by_vector, FALSE)
    by_levels <- if (isTRUE(rvcol)) rev(original_by_levels) else original_by_levels
    if (length(by_levels) < 2L) stop("The grouping variable must contain at least two observed levels.", call. = FALSE)
    by_factor <- factor(by_vector, levels = original_by_levels)
    valid_by <- !is.na(by_vector)
    group_total <- sum(valid_by)
    group_counts <- vapply(by_levels, function(level) sum(by_vector == level, na.rm = TRUE), numeric(1))
    group_percent <- if (group_total > 0L) 100 * group_counts / group_total else rep(NA_real_, length(group_counts))
    if (!is.null(effect_type) && length(original_by_levels) != 2L) stop("OR, RR, and PR require a binary outcome.", call. = FALSE)
    event_level <- if (is.null(event)) utils::tail(original_by_levels, 1L) else as.character(event)[1L]
    if (!event_level %in% original_by_levels) stop("`event` is not an outcome level.", call. = FALSE)
  } else {
    by_label <- NULL
    original_by_levels <- "Overall"
    by_levels <- "Overall"
    by_factor <- factor(rep("Overall", nrow(data)), levels = "Overall")
    valid_by <- rep(TRUE, nrow(data))
    group_total <- nrow(data)
    group_counts <- nrow(data)
    group_percent <- 100
    event_level <- NULL
  }

  # Descriptive summaries.
  summarize_mean <- function(x) {
    x <- x[!is.na(x)]
    if (!length(x)) return(c(first = "", second = ""))
    c(first = format_number(mean(x)), second = format_number(if (length(x) > 1L) stats::sd(x) else NA_real_))
  }
  summarize_median <- function(x) {
    x <- x[!is.na(x)]
    if (!length(x)) return(c(first = "", second = ""))
    q <- stats::quantile(x, c(.25, .5, .75), names = FALSE, type = 7)
    c(first = format_number(q[2L]), second = paste0(format_number(q[1L]), " - ", format_number(q[3L])))
  }
  summarize_range <- function(x) {
    x <- x[!is.na(x)]
    if (!length(x)) return(c(first = "", second = ""))
    c(first = format_number(min(x)), second = format_number(max(x)))
  }

  # Traditional omnibus tests.
  categorical_test <- function(x, g) {
    complete <- !is.na(x) & !is.na(g)
    tabulation <- table(x[complete], g[complete])
    if (nrow(tabulation) < 2L || ncol(tabulation) < 2L) return(list(p = NA_real_, method = NULL))
    chi <- suppressWarnings(stats::chisq.test(tabulation, correct = FALSE))
    if (any(chi$expected < 5)) {
      # Exact Fisher computation for large R x C tables can take minutes or
      # effectively never finish (a common accidental case is a continuous
      # numeric variable supplied without the c./q./f. prefix). Keep the exact
      # test for the canonical 2 x 2 table. For larger sparse tables use the
      # Fisher-Freeman-Halton Monte Carlo form, which preserves the conditional
      # margins but has bounded computation time. The RNG state is restored so
      # tab() never changes the user's random-number stream.
      if (identical(dim(tabulation), c(2L, 2L))) {
        fisher <- tryCatch(stats::fisher.test(tabulation), error = function(e) NULL)
        if (!is.null(fisher)) return(list(p = unname(fisher$p.value), method = "Fisher's exact test"))
      } else {
        had_seed <- exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE)
        if (had_seed) old_seed <- get(".Random.seed", envir = .GlobalEnv, inherits = FALSE)
        on.exit({
          if (had_seed) {
            assign(".Random.seed", old_seed, envir = .GlobalEnv)
          } else if (exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE)) {
            rm(".Random.seed", envir = .GlobalEnv)
          }
        }, add = TRUE)
        fisher <- tryCatch(
          stats::fisher.test(tabulation, simulate.p.value = TRUE, B = 5000L),
          error = function(e) NULL
        )
        if (!is.null(fisher)) {
          return(list(
            p = unname(fisher$p.value),
            method = "Fisher-Freeman-Halton test (Monte Carlo, B = 5000)"
          ))
        }
      }
    }
    list(p = unname(chi$p.value), method = "Pearson's chi-squared test")
  }
  mean_test <- function(x, g) {
    complete <- !is.na(x) & !is.na(g)
    x <- x[complete]
    g <- droplevels(factor(g[complete]))
    if (nlevels(g) < 2L) return(list(p = NA_real_, method = NULL))
    if (nlevels(g) == 2L) {
      split_x <- split(x, g)
      equal_variance <- FALSE
      if (all(vapply(split_x, length, integer(1)) >= 2L)) {
        variance_result <- tryCatch(stats::var.test(split_x[[1L]], split_x[[2L]]), error = function(e) NULL)
        equal_variance <- !is.null(variance_result) && is.finite(variance_result$p.value) && variance_result$p.value >= .05
      }
      result <- tryCatch(stats::t.test(x ~ g, var.equal = equal_variance), error = function(e) NULL)
      return(list(p = if (is.null(result)) NA_real_ else unname(result$p.value), method = if (equal_variance) "Student's t-test (equal variances)" else "Welch's t-test"))
    }
    p <- tryCatch(summary(stats::aov(x ~ g))[[1L]][["Pr(>F)"]][1L], error = function(e) NA_real_)
    list(p = unname(p), method = "One-way analysis of variance")
  }
  median_test <- function(x, g) {
    complete <- !is.na(x) & !is.na(g)
    x <- x[complete]
    g <- droplevels(factor(g[complete]))
    if (nlevels(g) < 2L) return(list(p = NA_real_, method = NULL))
    if (nlevels(g) == 2L) {
      result <- tryCatch(stats::wilcox.test(x ~ g, exact = FALSE), error = function(e) NULL)
      return(list(p = if (is.null(result)) NA_real_ else unname(result$p.value), method = "Wilcoxon rank-sum test"))
    }
    result <- tryCatch(stats::kruskal.test(x ~ g), error = function(e) NULL)
    list(p = if (is.null(result)) NA_real_ else unname(result$p.value), method = "Kruskal-Wallis test")
  }

  # Resolve categorical references from b-prefixes and backward-compatible effect_ref.
  reference_for <- function(variable, levels_original, index_from_prefix) {
    if (!is.null(effect_ref)) {
      if (is.list(effect_ref) && !is.null(names(effect_ref)) && variable %in% names(effect_ref)) {
        candidate <- as.character(effect_ref[[variable]])[1L]
        if (candidate %in% levels_original) return(candidate)
      }
      if (is.atomic(effect_ref) && !is.null(names(effect_ref)) && variable %in% names(effect_ref)) {
        candidate <- as.character(effect_ref[[variable]])[1L]
        if (candidate %in% levels_original) return(candidate)
      }
      if (length(effect_ref) == 1L && as.character(effect_ref)[1L] %in% levels_original) return(as.character(effect_ref)[1L])
    }
    if (is.na(index_from_prefix) || index_from_prefix < 1L || index_from_prefix > length(levels_original)) {
      stop(sprintf("The reference index is invalid for variable `%s`, which has %s observed levels.", variable, length(levels_original)), call. = FALSE)
    }
    as.character(levels_original[index_from_prefix])
  }

  # Robust covariance for modified Poisson regression.
  robust_vcov_poisson <- function(fit) {
    X <- stats::model.matrix(fit)
    mu <- stats::fitted(fit)
    score_residual <- fit$y - mu
    bread <- tryCatch(solve(crossprod(X, X * as.vector(mu))), error = function(e) NULL)
    if (is.null(bread)) return(NULL)
    meat <- crossprod(X, X * as.vector(score_residual^2))
    output <- bread %*% meat %*% bread
    dimnames(output) <- list(colnames(X), colnames(X))
    output
  }

  # Combine table and adjustment metadata, preserving explicit adjusted specs.
  combined_model_metadata <- function(focal_variable = NULL) {
    meta <- rbind(vars, adjusted_meta)
    meta <- meta[!duplicated(meta$variable, fromLast = TRUE), , drop = FALSE]
    if (!is.null(focal_variable)) {
      focal_index <- match(focal_variable, vars$variable)
      if (!is.na(focal_index)) {
        meta <- meta[meta$variable != focal_variable, , drop = FALSE]
        meta <- rbind(vars[focal_index, , drop = FALSE], meta)
      }
    }
    rownames(meta) <- NULL
    meta
  }

  # Prepare model variables and preserve categorical references from vars().
  prepare_model_data <- function(model_variables, model_meta) {
    model_data <- data[c(by_name, model_variables)]
    model_data <- model_data[stats::complete.cases(model_data), , drop = FALSE]
    if (!nrow(model_data)) return(NULL)
    model_data$.outcome <- as.integer(as.character(model_data[[by_name]]) == event_level)
    if (length(unique(model_data$.outcome)) < 2L) return(NULL)
    for (z in model_variables) {
      meta_index <- match(z, model_meta$variable)
      declared_type <- if (is.na(meta_index)) NULL else model_meta$type[meta_index]
      if (!is.null(declared_type) && declared_type == "categorical") {
        original <- get_levels(data[[z]], FALSE)
        ref <- reference_for(z, original, model_meta$reference_index[meta_index])
        model_data[[z]] <- factor(model_data[[z]], levels = original)
        model_data[[z]] <- stats::relevel(model_data[[z]], ref = ref)
      } else if (!is.null(declared_type) && declared_type %in% c("mean", "median", "full")) {
        model_data[[z]] <- as.numeric(model_data[[z]])
      } else if (is.factor(model_data[[z]]) || is.character(model_data[[z]]) || is.logical(model_data[[z]])) {
        model_data[[z]] <- factor(model_data[[z]])
      } else {
        model_data[[z]] <- as.numeric(model_data[[z]])
      }
    }
    model_data
  }

  # Track logistic-regression estimation problems that should be reported
  # publication-style as NA rather than as unstable Wald estimates.
  logistic_separation_detected <- FALSE

  mark_logistic_separation <- function() {
    logistic_separation_detected <<- TRUE
    invisible(NULL)
  }

  separation_effects <- function() {
    structure(list(), r4vn_logistic_separation = TRUE)
  }

  has_separation_effects <- function(x) {
    isTRUE(attr(x, "r4vn_logistic_separation", exact = TRUE))
  }

  # Fit a regression model while capturing the two standard glm() warnings
  # produced by complete or quasi-complete separation. These warnings are
  # intentionally muffled because tab() converts them into an informative
  # publication note and displays NA for the unreliable effect estimate.
  fit_effect_model <- function(formula, model_data) {
    if (identical(effect_type, "OR")) {
      fit_warnings <- character()
      fit <- withCallingHandlers(
        tryCatch(
          stats::glm(
            formula,
            family = stats::binomial("logit"),
            data = model_data,
            y = TRUE
          ),
          error = function(e) NULL
        ),
        warning = function(w) {
          fit_warnings <<- c(fit_warnings, conditionMessage(w))
          invokeRestart("muffleWarning")
        }
      )

      separation_warning <- any(grepl(
        "algorithm did not converge|fitted probabilities numerically 0 or 1 occurred",
        fit_warnings,
        ignore.case = TRUE
      ))

      separation <- !is.null(fit) && (
        !isTRUE(fit$converged) || separation_warning
      )

      return(list(
        fit = fit,
        separation = separation,
        warnings = fit_warnings
      ))
    }

    fit <- tryCatch(
      stats::glm(
        formula,
        family = stats::poisson("log"),
        data = model_data,
        y = TRUE
      ),
      error = function(e) NULL
    )

    list(
      fit = fit,
      separation = FALSE,
      warnings = character()
    )
  }

  # Fit one crude or adjusted model and return terms for the focal predictor.
  model_effects <- function(predictor_name, categorical, reference = NULL, adjustment_metadata = NULL) {
    if (is.null(adjustment_metadata)) adjustment_metadata <- adjusted_meta[0, , drop = FALSE]
    adjustment_variables <- adjustment_metadata$variable
    model_variables <- unique(c(predictor_name, adjustment_variables))
    model_variables <- setdiff(model_variables, by_name)
    model_meta <- combined_model_metadata(predictor_name)
    model_data <- prepare_model_data(model_variables, model_meta)
    if (is.null(model_data)) return(list())
    if (!categorical) {
      model_data[[predictor_name]] <- as.numeric(model_data[[predictor_name]])
      if (!is.finite(stats::sd(model_data[[predictor_name]])) || stats::sd(model_data[[predictor_name]]) == 0) return(list())
    } else {
      original <- get_levels(data[[predictor_name]], FALSE)
      model_data[[predictor_name]] <- factor(model_data[[predictor_name]], levels = original)
      model_data[[predictor_name]] <- stats::relevel(model_data[[predictor_name]], ref = reference)
    }
    formula <- stats::reformulate(model_variables, response = ".outcome")
    fit_result <- fit_effect_model(formula, model_data)
    fit <- fit_result$fit
    if (is.null(fit)) return(list())

    if (identical(effect_type, "OR") && isTRUE(fit_result$separation)) {
      mark_logistic_separation()
      return(separation_effects())
    }

    beta <- stats::coef(fit)
    covariance <- if (effect_type == "OR") tryCatch(stats::vcov(fit), error = function(e) NULL) else robust_vcov_poisson(fit)
    if (is.null(covariance)) return(list())
    term_names <- names(beta)
    selected <- if (categorical) startsWith(term_names, predictor_name) & term_names != predictor_name else term_names == predictor_name
    selected <- selected & term_names != "(Intercept)" & is.finite(beta)
    if (!any(selected)) return(list())
    beta <- beta[selected]
    term_names <- term_names[selected]
    covariance <- covariance[term_names, term_names, drop = FALSE]
    se <- sqrt(diag(covariance))
    valid <- is.finite(beta) & is.finite(se) & se > 0
    beta <- beta[valid]
    se <- se[valid]
    term_names <- term_names[valid]
    if (!length(beta)) return(list())
    estimate <- exp(beta)
    lower <- exp(beta - 1.96 * se)
    upper <- exp(beta + 1.96 * se)
    p <- 2 * stats::pnorm(abs(beta / se), lower.tail = FALSE)
    output <- vector("list", length(beta))
    names(output) <- term_names
    for (j in seq_along(beta)) {
      output[[j]] <- list(
        estimate = unname(estimate[j]),
        lower = unname(lower[j]),
        upper = unname(upper[j]),
        p = unname(p[j]),
        text = paste0(format_number(estimate[j], effect_digit), " (", format_number(lower[j], effect_digit), " - ", format_number(upper[j], effect_digit), ")")
      )
    }
    output
  }

  # Calculate coefficient-level variance inflation factors from the model matrix.
  model_vif_range <- function(fit) {
    X <- stats::model.matrix(fit)
    if (ncol(X) <= 2L) return(c(min = 1, max = 1))
    X <- X[, colnames(X) != "(Intercept)", drop = FALSE]
    variable_ok <- apply(X, 2L, function(z) is.finite(stats::sd(z)) && stats::sd(z) > 0)
    X <- X[, variable_ok, drop = FALSE]
    if (ncol(X) < 2L) return(c(min = 1, max = 1))
    correlation <- suppressWarnings(stats::cor(X))
    inverse <- tryCatch(solve(correlation), error = function(e) tryCatch(qr.solve(correlation), error = function(e2) NULL))
    if (is.null(inverse)) return(c(min = NA_real_, max = NA_real_))
    values <- diag(inverse)
    values <- values[is.finite(values) & values >= 1]
    if (!length(values)) return(c(min = NA_real_, max = NA_real_))
    c(min = min(values), max = max(values))
  }

  # Hosmer-Lemeshow goodness-of-fit test implemented without external packages.
  hosmer_lemeshow_test <- function(y, fitted, groups = 10L) {
    complete <- is.finite(y) & is.finite(fitted)
    y <- y[complete]
    fitted <- fitted[complete]
    if (length(y) < 20L || length(unique(fitted)) < 3L) return(list(statistic = NA_real_, df = NA_integer_, p = NA_real_, groups = NA_integer_))
    groups <- min(as.integer(groups), max(2L, floor(length(y) / 5L)))
    breaks <- unique(stats::quantile(fitted, probs = seq(0, 1, length.out = groups + 1L), na.rm = TRUE, names = FALSE))
    if (length(breaks) < 3L) return(list(statistic = NA_real_, df = NA_integer_, p = NA_real_, groups = NA_integer_))
    group <- cut(fitted, breaks = breaks, include.lowest = TRUE, labels = FALSE)
    observed <- rowsum(y, group, reorder = FALSE)
    expected <- rowsum(fitted, group, reorder = FALSE)
    number <- as.numeric(table(group))
    observed <- as.numeric(observed)
    expected <- as.numeric(expected)
    denominator_event <- pmax(expected, .Machine$double.eps)
    denominator_nonevent <- pmax(number - expected, .Machine$double.eps)
    statistic <- sum((observed - expected)^2 / denominator_event + ((number - observed) - (number - expected))^2 / denominator_nonevent)
    df <- max(1L, length(number) - 2L)
    list(statistic = statistic, df = df, p = stats::pchisq(statistic, df = df, lower.tail = FALSE), groups = length(number))
  }

  # Fit the single final multivariable model requested by multi().
  fit_multi_model <- function() {
    if (!has_multi) return(NULL)
    model_variables <- multi_meta$variable
    model_meta <- rbind(vars, multi_meta)
    model_meta <- model_meta[!duplicated(model_meta$variable, fromLast = TRUE), , drop = FALSE]
    model_data <- prepare_model_data(model_variables, model_meta)
    if (is.null(model_data)) return(NULL)
    formula <- stats::reformulate(model_variables, response = ".outcome")
    fit_result <- fit_effect_model(formula, model_data)
    fit <- fit_result$fit
    if (is.null(fit)) return(NULL)

    if (identical(effect_type, "OR") && isTRUE(fit_result$separation)) {
      mark_logistic_separation()
      return(list(
        fit = fit,
        covariance = NULL,
        effects = separation_effects(),
        diagnostics = list(
          r2_name = "Nagelkerke R2",
          r2 = NA_real_,
          gof_name = "Hosmer-Lemeshow",
          gof_statistic = NA_real_,
          gof_df = NA_integer_,
          gof_p = NA_real_,
          vif = c(min = NA_real_, max = NA_real_),
          n = stats::nobs(fit)
        ),
        metadata = multi_meta,
        separation = TRUE
      ))
    }

    covariance <- if (effect_type == "OR") tryCatch(stats::vcov(fit), error = function(e) NULL) else robust_vcov_poisson(fit)
    if (is.null(covariance)) return(NULL)
    beta <- stats::coef(fit)
    term_names <- names(beta)
    selected <- term_names != "(Intercept)" & is.finite(beta)
    beta_selected <- beta[selected]
    terms_selected <- term_names[selected]
    covariance_selected <- covariance[terms_selected, terms_selected, drop = FALSE]
    se <- sqrt(diag(covariance_selected))
    valid <- is.finite(beta_selected) & is.finite(se) & se > 0
    beta_selected <- beta_selected[valid]
    se <- se[valid]
    terms_selected <- terms_selected[valid]
    effects <- vector("list", length(beta_selected))
    names(effects) <- terms_selected
    if (length(beta_selected)) {
      estimate <- exp(beta_selected)
      lower <- exp(beta_selected - 1.96 * se)
      upper <- exp(beta_selected + 1.96 * se)
      p <- 2 * stats::pnorm(abs(beta_selected / se), lower.tail = FALSE)
      for (j in seq_along(beta_selected)) {
        effects[[j]] <- list(
          estimate = unname(estimate[j]), lower = unname(lower[j]), upper = unname(upper[j]), p = unname(p[j]),
          text = paste0(format_number(estimate[j], effect_digit), " (", format_number(lower[j], effect_digit), " - ", format_number(upper[j], effect_digit), ")")
        )
      }
    }
    vif <- model_vif_range(fit)
    if (effect_type == "OR") {
      null_fit <- tryCatch(stats::glm(.outcome ~ 1, family = stats::binomial("logit"), data = model_data), error = function(e) NULL)
      r2 <- NA_real_
      if (!is.null(null_fit)) {
        ll_model <- as.numeric(stats::logLik(fit))
        ll_null <- as.numeric(stats::logLik(null_fit))
        n <- stats::nobs(fit)
        denominator <- 1 - exp(2 * ll_null / n)
        if (is.finite(denominator) && denominator != 0) r2 <- (1 - exp(2 * (ll_null - ll_model) / n)) / denominator
      }
      gof <- hosmer_lemeshow_test(model_data$.outcome, stats::fitted(fit), groups = 10L)
      diagnostics <- list(
        r2_name = "Nagelkerke R2", r2 = r2,
        gof_name = "Hosmer-Lemeshow", gof_statistic = gof$statistic,
        gof_df = gof$df, gof_p = gof$p, vif = vif, n = stats::nobs(fit)
      )
    } else {
      r2 <- if (is.finite(fit$null.deviance) && fit$null.deviance > 0) 1 - fit$deviance / fit$null.deviance else NA_real_
      pearson <- sum(stats::residuals(fit, type = "pearson")^2, na.rm = TRUE)
      pearson_df <- stats::df.residual(fit)
      diagnostics <- list(
        r2_name = "Deviance R2", r2 = r2,
        gof_name = "Pearson goodness-of-fit", gof_statistic = pearson,
        gof_df = pearson_df, gof_p = if (pearson_df > 0) stats::pchisq(pearson, df = pearson_df, lower.tail = FALSE) else NA_real_,
        vif = vif, n = stats::nobs(fit)
      )
    }
    list(
      fit = fit,
      covariance = covariance,
      effects = effects,
      diagnostics = diagnostics,
      metadata = multi_meta,
      separation = FALSE
    )
  }

  multi_model <- fit_multi_model()

  # Match a displayed categorical level to its regression term.
  find_category_effect <- function(effect_list, variable, level) {
    if (!length(effect_list)) return(NULL)
    candidate <- paste0(variable, level)
    if (candidate %in% names(effect_list)) return(effect_list[[candidate]])
    matches <- names(effect_list)[endsWith(names(effect_list), as.character(level))]
    if (length(matches) == 1L) effect_list[[matches]] else NULL
  }

  # Standard internal row.
  make_row <- function(variable, label, item, type, kind, group_values,
                       overall_value = NULL, test_p = NA_real_, test_method = NULL,
                       crude = "", crude_p = NA_real_, adjusted = "",
                       adjusted_p = NA_real_, multi = "", multi_p = NA_real_,
                       variable_start = FALSE, raw_values = NULL) {
    list(
      variable = variable, label = label, item = item, type = type, kind = kind,
      group_values = group_values, overall_value = overall_value,
      test_p = test_p, test_method = test_method, crude = crude,
      crude_p = crude_p, adjusted = adjusted, adjusted_p = adjusted_p,
      multi = multi, multi_p = multi_p,
      variable_start = variable_start, raw = raw_values
    )
  }

  calculate_percent <- function(level_mask, group_mask, variable_observed) {
    numerator <- sum(level_mask & group_mask & valid_by)
    denominator <- if (!has_by) {
      sum(variable_observed & valid_by)
    } else if (isTRUE(row)) {
      sum(level_mask & valid_by)
    } else if (isTRUE(col)) {
      sum(group_mask & variable_observed & valid_by)
    } else {
      sum(variable_observed & valid_by)
    }
    if (denominator > 0L) 100 * numerator / denominator else NA_real_
  }

  # Calculate rows, tests, crude models, and adjusted models.
  table_rows <- list()
  for (i in seq_len(nrow(vars))) {
    variable <- vars$variable[i]
    summary_type <- vars$type[i]
    x <- data[[variable]]
    base_label <- get_label(x, variable)
    label <- display_label(base_label, variable)
    omnibus <- list(p = NA_real_, method = NULL)
    if (has_by && isTRUE(test)) {
      omnibus <- if (summary_type == "categorical") categorical_test(x, by_factor) else if (summary_type == "mean") mean_test(x, by_factor) else median_test(x, by_factor)
    }

    # Select adjustment variables separately for each focal predictor.
    focal_adjusted_meta <- adjusted_meta[adjusted_meta$variable != variable, , drop = FALSE]
    if (!is.null(by_name)) focal_adjusted_meta <- focal_adjusted_meta[focal_adjusted_meta$variable != by_name, , drop = FALSE]
    variable_in_multi <- has_multi && variable %in% multi_meta$variable && !is.null(multi_model)

    # Numeric variables.
    if (summary_type %in% c("mean", "median", "full")) {
      if (!is.numeric(x)) stop(sprintf("Variable `%s` must be numeric.", variable), call. = FALSE)
      crude_effects <- if (!is.null(effect_type)) model_effects(variable, FALSE, adjustment_metadata = adjusted_meta[0, , drop = FALSE]) else list()
      adjusted_effects <- if (has_adjusted) model_effects(variable, FALSE, adjustment_metadata = focal_adjusted_meta) else list()

      crude_entry <- if (has_separation_effects(crude_effects)) {
        list(text = "NA", p = NA_real_)
      } else if (length(crude_effects)) {
        crude_effects[[1L]]
      } else NULL

      adjusted_entry <- if (has_separation_effects(adjusted_effects)) {
        list(text = "NA", p = NA_real_)
      } else if (length(adjusted_effects)) {
        adjusted_effects[[1L]]
      } else NULL

      multi_entry <- if (variable_in_multi && isTRUE(multi_model$separation)) {
        list(text = "NA", p = NA_real_)
      } else if (variable_in_multi && variable %in% names(multi_model$effects)) {
        multi_model$effects[[variable]]
      } else NULL
      if (summary_type == "full") {
        table_rows[[length(table_rows) + 1L]] <- make_row(
          variable, label, "", "numeric_header", "", vector("list", length(by_levels)),
          test_p = omnibus$p, test_method = omnibus$method,
          crude = if (is.null(crude_entry)) "" else crude_entry$text,
          crude_p = if (is.null(crude_entry)) NA_real_ else crude_entry$p,
          adjusted = if (is.null(adjusted_entry)) "" else adjusted_entry$text,
          adjusted_p = if (is.null(adjusted_entry)) NA_real_ else adjusted_entry$p,
          multi = if (is.null(multi_entry)) "" else multi_entry$text,
          multi_p = if (is.null(multi_entry)) NA_real_ else multi_entry$p,
          variable_start = TRUE
        )
        methods <- c("mean", "median", "range")
      } else {
        methods <- summary_type
      }
      for (method in methods) {
        item <- if (summary_type == "full") switch(method, mean = "Mean (SD)", median = "Median (IQR)", range = "Range") else ""
        shown_label <- if (summary_type == "mean") {
          display_label(paste0(base_label, ", M (SD)"), variable)
        } else if (summary_type == "median") {
          display_label(paste0(base_label, ", Median (IQR)"), variable)
        } else {
          label
        }
        group_values <- vector("list", length(by_levels))
        raw_values <- vector("list", length(by_levels))
        for (g_index in seq_along(by_levels)) {
          subset_x <- x[by_factor == by_levels[g_index] & valid_by]
          group_values[[g_index]] <- switch(method, mean = summarize_mean(subset_x), median = summarize_median(subset_x), range = summarize_range(subset_x))
          raw_values[[g_index]] <- list(values = subset_x)
        }
        overall_value <- switch(method, mean = summarize_mean(x[valid_by]), median = summarize_median(x[valid_by]), range = summarize_range(x[valid_by]))
        attach_model <- summary_type != "full"
        table_rows[[length(table_rows) + 1L]] <- make_row(
          variable, shown_label, item,
          if (summary_type == "full") "numeric_detail" else "numeric",
          method, group_values, overall_value,
          test_p = if (attach_model) omnibus$p else NA_real_,
          test_method = if (attach_model) omnibus$method else NULL,
          crude = if (attach_model && !is.null(crude_entry)) crude_entry$text else "",
          crude_p = if (attach_model && !is.null(crude_entry)) crude_entry$p else NA_real_,
          adjusted = if (attach_model && !is.null(adjusted_entry)) adjusted_entry$text else "",
          adjusted_p = if (attach_model && !is.null(adjusted_entry)) adjusted_entry$p else NA_real_,
          multi = if (attach_model && !is.null(multi_entry)) multi_entry$text else "",
          multi_p = if (attach_model && !is.null(multi_entry)) multi_entry$p else NA_real_,
          variable_start = attach_model, raw_values = raw_values
        )
      }
      next
    }

    # Categorical variables.
    levels_original <- get_levels(x, FALSE)
    reference <- reference_for(variable, levels_original, vars$reference_index[i])
    multi_reference <- NULL
    if (variable_in_multi) {
      multi_index <- match(variable, multi_meta$variable)
      multi_reference <- reference_for(variable, levels_original, multi_meta$reference_index[multi_index])
    }
    levels_to_show <- get_levels(x, variable %in% reverse_rows)
    crude_effects <- if (!is.null(effect_type)) model_effects(variable, TRUE, reference, adjusted_meta[0, , drop = FALSE]) else list()
    adjusted_effects <- if (has_adjusted) model_effects(variable, TRUE, reference, focal_adjusted_meta) else list()
    multi_effects <- if (variable_in_multi) multi_model$effects else list()

    table_rows[[length(table_rows) + 1L]] <- make_row(
      variable, label, "", "categorical_header", "", vector("list", length(by_levels)),
      test_p = omnibus$p, test_method = omnibus$method, variable_start = TRUE
    )
    variable_observed <- !is.na(x)
    for (level_value in levels_to_show) {
      level_mask <- !is.na(x) & x == level_value
      group_values <- vector("list", length(by_levels))
      raw_values <- vector("list", length(by_levels))
      for (g_index in seq_along(by_levels)) {
        group_mask <- by_factor == by_levels[g_index] & valid_by
        count <- sum(level_mask & group_mask & valid_by)
        percent <- calculate_percent(level_mask, group_mask, variable_observed)
        group_values[[g_index]] <- c(first = format_count(count), second = format_number(percent))
        raw_values[[g_index]] <- list(n = count, percent = percent)
      }
      overall_count <- sum(level_mask & valid_by)
      overall_denominator <- sum(variable_observed & valid_by)
      overall_value <- c(first = format_count(overall_count), second = format_number(if (overall_denominator > 0L) 100 * overall_count / overall_denominator else NA_real_))
      crude_text <- adjusted_text <- multi_text <- ""
      crude_p <- adjusted_p <- multi_p <- NA_real_
      if (!is.null(effect_type)) {
        if (as.character(level_value) == reference) {
          crude_text <- "Ref"
          if (has_adjusted) adjusted_text <- "Ref"
        } else {
          if (has_separation_effects(crude_effects)) {
            crude_text <- "NA"
            crude_p <- NA_real_
          } else {
            crude_entry <- find_category_effect(crude_effects, variable, level_value)
            if (!is.null(crude_entry)) {
              crude_text <- crude_entry$text
              crude_p <- crude_entry$p
            }
          }

          if (has_separation_effects(adjusted_effects)) {
            adjusted_text <- "NA"
            adjusted_p <- NA_real_
          } else {
            adjusted_entry <- find_category_effect(adjusted_effects, variable, level_value)
            if (!is.null(adjusted_entry)) {
              adjusted_text <- adjusted_entry$text
              adjusted_p <- adjusted_entry$p
            }
          }
        }
        if (variable_in_multi) {
          if (as.character(level_value) == multi_reference) {
            multi_text <- "Ref"
          } else if (isTRUE(multi_model$separation)) {
            multi_text <- "NA"
            multi_p <- NA_real_
          } else {
            multi_category_entry <- find_category_effect(multi_effects, variable, level_value)
            if (!is.null(multi_category_entry)) {
              multi_text <- multi_category_entry$text
              multi_p <- multi_category_entry$p
            }
          }
        }
      }
      table_rows[[length(table_rows) + 1L]] <- make_row(
        variable, label, as.character(level_value), "categorical_level", "categorical",
        group_values, overall_value, crude = crude_text, crude_p = crude_p,
        adjusted = adjusted_text, adjusted_p = adjusted_p,
        multi = multi_text, multi_p = multi_p, raw_values = raw_values
      )
    }

    # Missing category.
    missing_count <- sum(is.na(x) & valid_by)
    include_missing <- identical(missing, "always") || (identical(missing, "ifany") && missing_count > 0L)
    if (include_missing) {
      group_values <- vector("list", length(by_levels))
      raw_values <- vector("list", length(by_levels))
      for (g_index in seq_along(by_levels)) {
        group_mask <- by_factor == by_levels[g_index] & valid_by
        count <- sum(is.na(x) & group_mask)
        denominator <- if (!has_by) sum(valid_by) else if (isTRUE(row)) missing_count else if (isTRUE(col)) sum(group_mask & valid_by) else sum(valid_by)
        percent <- if (denominator > 0L) 100 * count / denominator else NA_real_
        group_values[[g_index]] <- c(first = format_count(count), second = format_number(percent))
        raw_values[[g_index]] <- list(n = count, percent = percent)
      }
      overall_value <- c(first = format_count(missing_count), second = format_number(if (sum(valid_by) > 0L) 100 * missing_count / sum(valid_by) else NA_real_))
      table_rows[[length(table_rows) + 1L]] <- make_row(
        variable, label, "Missing", "missing", "categorical",
        group_values, overall_value, raw_values = raw_values
      )
    }
  }

  # Superscript letters for traditional tests.
  tests_used <- unique(vapply(table_rows, function(z) if (is.null(z$test_method)) "" else z$test_method, character(1)))
  tests_used <- tests_used[nzchar(tests_used)]
  test_letters <- stats::setNames(c(letters, paste0("a", letters))[seq_along(tests_used)], tests_used)

  # Descriptive result columns and headers.
  result_columns <- list()
  if (isTRUE(descriptive)) {
    overall_header <- paste0("Overall<span class=\"header-n\">n = ", format_count(sum(valid_by)), "; ", format_number(100), "%</span>")
    if (has_by) {
      if (overall == "first") result_columns <- c(result_columns, list(list(type = "overall", label = "Overall", html = overall_header)))
      for (j in seq_along(by_levels)) {
        header <- paste0(escape_html(by_levels[j]), "<span class=\"header-n\">n = ", format_count(group_counts[j]), "; ", format_number(group_percent[j]), "%</span>")
        result_columns <- c(result_columns, list(list(type = "group", label = as.character(by_levels[j]), html = header)))
      }
      if (overall == "last") result_columns <- c(result_columns, list(list(type = "overall", label = "Overall", html = overall_header)))
    } else {
      result_columns <- list(list(type = "overall", label = "Overall", html = overall_header))
    }
  }

  result_header_cells <- paste0("<th class=\"result-head\">", vapply(result_columns, `[[`, character(1), "html"), "</th>", collapse = "")
  crude_heading <- if (!is.null(effect_type)) paste0("Crude ", effect_type, " (95% CI)") else NULL
  adjusted_heading <- if (has_adjusted) paste0("Adjusted ", effect_type, " (95% CI)") else NULL
  multi_heading <- if (has_multi) paste0("Multivariable ", effect_type, " (95% CI)") else NULL
  if (!has_by) {
    header_html <- paste0("<thead><tr><th>Characteristic</th>", result_header_cells, "</tr></thead>")
  } else {
    extra <- as.integer(test) + as.integer(!is.null(effect_type)) + as.integer(!is.null(effect_type) && pvalue) + as.integer(has_adjusted) + as.integer(has_adjusted && pvalue) + as.integer(has_multi) + as.integer(has_multi && pvalue)
    header_html <- paste0(
      "<thead><tr class=\"by-title\"><th></th><th colspan=\"", length(result_columns) + extra, "\">",
      escape_html(by_label), if (name) paste0(" <span class=\"variable-code\">[", escape_html(by_name), "]</span>") else "",
      "</th></tr><tr><th>Characteristic</th>", result_header_cells,
      if (test) "<th>Test p</th>" else "",
      if (!is.null(effect_type)) paste0("<th>", crude_heading, "</th>") else "",
      if (!is.null(effect_type) && pvalue) "<th>p-value</th>" else "",
      if (has_adjusted) paste0("<th>", adjusted_heading, "</th>") else "",
      if (has_adjusted && pvalue) "<th>p-value</th>" else "",
      if (has_multi) paste0("<th>", multi_heading, "</th>") else "",
      if (has_multi && pvalue) "<th>p-value</th>" else "",
      "</tr></thead>"
    )
  }

  # Render rows.
  body_html <- character(length(table_rows))
  for (i in seq_along(table_rows)) {
    current <- table_rows[[i]]
    characteristic <- if (current$type %in% c("categorical_header", "numeric_header", "numeric")) {
      paste0("<span class=\"variable-name\">", current$label, "</span>")
    } else {
      paste0("<span class=\"level-name", if (current$type == "missing") " missing-name" else "", "\">", escape_html(current$item), "</span>")
    }
    value_cells <- vapply(result_columns, function(column) {
      value <- if (column$type == "overall") current$overall_value else current$group_values[[match(column$label, by_levels)]]
      text <- if (is.null(value) || current$type %in% c("categorical_header", "numeric_header")) "" else compact_value(value, current$kind)
      paste0("<td class=\"result\">", escape_html(text), "</td>")
    }, character(1))
    test_cell <- if (has_by && test) {
      p <- format_p(current$test_p)
      if (nzchar(p) && test_note && !is.null(current$test_method)) p <- paste0(p, "<sup>", test_letters[[current$test_method]], "</sup>")
      paste0("<td class=\"test-p\">", p, "</td>")
    } else ""
    crude_cell <- if (!is.null(effect_type)) paste0("<td class=\"effect\">", escape_html(current$crude), "</td>") else ""
    crude_p_cell <- if (!is.null(effect_type) && pvalue) paste0("<td class=\"effect-p\">", format_p(current$crude_p), "</td>") else ""
    adjusted_cell <- if (has_adjusted) paste0("<td class=\"effect\">", escape_html(current$adjusted), "</td>") else ""
    adjusted_p_cell <- if (has_adjusted && pvalue) paste0("<td class=\"effect-p\">", format_p(current$adjusted_p), "</td>") else ""
    multi_cell <- if (has_multi) paste0("<td class=\"effect\">", escape_html(current$multi), "</td>") else ""
    multi_p_cell <- if (has_multi && pvalue) paste0("<td class=\"effect-p\">", format_p(current$multi_p), "</td>") else ""
    row_class <- paste0("row-", gsub("_", "-", current$type, fixed = TRUE), if (current$variable_start) " variable-start" else "")
    body_html[i] <- paste0("<tr class=\"", row_class, "\"><td>", characteristic, "</td>", paste(value_cells, collapse = ""), test_cell, crude_cell, crude_p_cell, adjusted_cell, adjusted_p_cell, multi_cell, multi_p_cell, "</tr>")
  }

  # Notes and footnotes.
  percentage_note <- if (!has_by) {
    "Percentages use non-missing observations as the denominator."
  } else if (row) {
    "Categorical percentages are calculated by row."
  } else if (col) {
    "Categorical percentages are calculated by column."
  } else {
    "Categorical percentages are calculated using the complete table total."
  }
  test_footnotes <- if (test && test_note && length(tests_used)) {
    entries <- paste0("<sup>", unname(test_letters[tests_used]), "</sup> ", escape_html(tests_used))
    paste0("<div class=\"test-note\">", paste(entries, collapse = "; "), ".</div>")
  } else ""
  adjustment_note <- if (has_adjusted) {
    if (adjust_all) "Adjusted estimates include all other variables listed in the table." else paste0("Adjusted estimates include: ", paste(escape_html(adjusted_meta$specification), collapse = ", "), ".")
  } else ""
  effect_note <- if (!is.null(effect_type)) {
    paste0(
      "<div class=\"effect-note\">Event = ", escape_html(event_level),
      ". Categorical reference categories are marked Ref. Continuous estimates are reported per one-unit increase. ",
      if (effect_type == "OR") "Logistic regression was used." else "Modified Poisson regression with robust variance was used.",
      if (nzchar(adjustment_note)) paste0(" ", adjustment_note) else "", "</div>"
    )
  } else ""

  estimation_note <- if (identical(effect_type, "OR") && isTRUE(logistic_separation_detected)) {
    paste0(
      "<div class=\"effect-note\">",
      "NA: estimate not available. Logistic regression estimates are reported as NA when complete or quasi-complete separation prevents reliable estimation.",
      "</div>"
    )
  } else ""

  descriptive_note <- if (isTRUE(descriptive)) paste0("<div class=\"table-note\">M: mean; SD: standard deviation; IQR: interquartile range. ", percentage_note, "</div>") else ""
  multi_note <- ""
  if (has_multi && !is.null(multi_model) && !isTRUE(multi_model$separation)) {
    diagnostics <- multi_model$diagnostics
    vif_text <- if (all(is.finite(diagnostics$vif))) paste0(format_number(diagnostics$vif[1L], 2L), " - ", format_number(diagnostics$vif[2L], 2L)) else "not estimable"
    r2_text <- if (is.finite(diagnostics$r2)) format_number(diagnostics$r2, 3L) else "not estimable"
    gof_text <- if (is.finite(diagnostics$gof_p)) paste0("p = ", if (diagnostics$gof_p < 10^(-p_digit)) paste0("&lt;", formatC(10^(-p_digit), format = "f", digits = p_digit)) else formatC(diagnostics$gof_p, format = "f", digits = p_digit)) else "p not estimable"
    multi_note <- paste0(
      "<div class=\"model-note\">Final multivariable model: ", paste(escape_html(multi_meta$specification), collapse = ", "),
      ". n = ", format_count(diagnostics$n), "; ", diagnostics$r2_name, " = ", r2_text,
      "; ", diagnostics$gof_name, " (", if (is.finite(diagnostics$gof_df)) paste0("df = ", diagnostics$gof_df, ", ") else "", gof_text,
      "); coefficient-level VIF range = ", vif_text, ".</div>"
    )
  }
  note_html <- paste0(descriptive_note, test_footnotes, effect_note, estimation_note, multi_note)
  title_html <- if (is.null(title) || !nzchar(as.character(title)[1L])) "" else paste0("<div class=\"table-title\">", escape_html(as.character(title)[1L]), "</div>")

  # CSS templates.
  css <- switch(
    template,
    journal = "body{font-family:'Times New Roman',Times,serif;background:#fff;color:#111;margin:18px}.table-title{font-size:18px;font-weight:700;margin:0 0 8px}table{border-collapse:collapse;width:auto;min-width:820px;border-top:2px solid #111;border-bottom:2px solid #111}th{padding:4px 9px;text-align:right;border-bottom:1.5px solid #111;font-weight:700;white-space:nowrap;background:#fff}th:first-child{text-align:left;min-width:260px}td{padding:3px 9px;vertical-align:top;border:0}tr.variable-start td{border-top:1px solid #aaa}td.result,td.test-p,td.effect,td.effect-p{text-align:right;white-space:nowrap}.by-title th{text-align:center;border-bottom:1px solid #777}",
    clean = "body{font-family:Arial,Helvetica,sans-serif;background:#fff;color:#111;margin:18px}.table-title{font-size:18px;font-weight:700;margin:0 0 8px}table{border-collapse:collapse;width:auto;min-width:820px;border-top:2px solid #222;border-bottom:2px solid #222}th{padding:6px 10px;text-align:right;border-bottom:1.5px solid #222;font-weight:700;white-space:nowrap;background:#f3f3f3}th:first-child{text-align:left;min-width:260px}td{padding:4px 10px;vertical-align:top;border-bottom:1px solid #ddd}tr.variable-start td{border-top:1px solid #999}td.result,td.test-p,td.effect,td.effect-p{text-align:right;white-space:nowrap}.by-title th{text-align:center;background:#fff;border-bottom:1px solid #999}",
    minimal = "body{font-family:Arial,Helvetica,sans-serif;background:#fff;color:#111;margin:18px}.table-title{font-size:17px;font-weight:700;margin:0 0 7px}table{border-collapse:collapse;width:auto;min-width:780px;border-top:1.5px solid #222;border-bottom:1.5px solid #222}th{padding:4px 8px;text-align:right;border-bottom:1px solid #555;font-weight:700;white-space:nowrap;background:#fff}th:first-child{text-align:left;min-width:240px}td{padding:3px 8px;vertical-align:top;border:0}tr.variable-start td{border-top:1px solid #ddd}td.result,td.test-p,td.effect,td.effect-p{text-align:right;white-space:nowrap}.by-title th{text-align:center;border-bottom:1px solid #aaa}"
  )
  common_css <- ".table-wrapper{display:inline-block;max-width:100%;overflow-x:auto}.variable-name{font-weight:700}.variable-code{font-family:Consolas,monospace;font-size:.78em;color:#666;font-weight:400}.level-name{display:inline-block;padding-left:22px;white-space:nowrap}.missing-name{font-style:italic}.header-n{display:block;font-size:.78em;font-weight:400;text-align:right;margin-top:1px}.result-head{text-align:right}.table-note,.test-note,.effect-note,.model-note{font-size:12px;color:#333;margin-top:6px;line-height:1.35}.table-separator{height:24px}sup{font-size:.72em;vertical-align:super;margin-left:1px}"

  # Build and optionally append the HTML document.
  table_block <- paste0("<section class=\"r4vn-table\">", title_html, "<table>", header_html, "<tbody>", paste(body_html, collapse = ""), "</tbody></table>", note_html, "</section>")
  blocks <- table_block
  if (inherits(append, "r4vn_tab")) blocks <- c(append$blocks, table_block)
  document <- paste0("<!DOCTYPE html><html><head><meta charset=\"UTF-8\"><meta name=\"viewport\" content=\"width=device-width,initial-scale=1\"><style>", css, common_css, "</style></head><body><div class=\"table-wrapper\">", paste(blocks, collapse = "<div class=\"table-separator\"></div>"), "</div></body></html>")
  if (is.null(file)) file <- tempfile(pattern = "r4vn-tab-", fileext = ".html")
  if (!is.character(file) || length(file) != 1L || !nzchar(file)) stop("`file` must be a single valid file path.", call. = FALSE)
  if (is.character(append) && length(append) == 1L && file.exists(append)) {
    old <- paste(readLines(append, warn = FALSE, encoding = "UTF-8"), collapse = "\n")
    if (grepl("</body>", old, fixed = TRUE)) {
      document <- sub("</body>", paste0("<div class=\"table-separator\"></div>", table_block, "</body>"), old, fixed = TRUE)
      file <- append
    }
  }
  writeLines(enc2utf8(document), file, useBytes = TRUE)

  raw_output <- if (isTRUE(raw)) lapply(table_rows, function(z) {
    list(
      variable = z$variable, label = z$label, item = z$item, type = z$type,
      group_values = z$raw, overall = z$overall_value, test_p = z$test_p,
      test_method = z$test_method, crude = z$crude, crude_p = z$crude_p,
      adjusted = z$adjusted, adjusted_p = z$adjusted_p,
      multi = z$multi, multi_p = z$multi_p
    )
  }) else NULL
  # Build a flat data frame for Word/Excel directly from the internal rows.
  # This block is intentionally local so the package does not need a separate
  # tabledata()/build_table_dataframe() helper.
  clean_export_text <- function(x) {
    x <- as.character(x)
    x <- gsub("<br\\s*/?>", " ", x, ignore.case = TRUE)
    x <- gsub("<[^>]+>", "", x)
    x <- gsub("&lt;", "<", x, fixed = TRUE)
    x <- gsub("&gt;", ">", x, fixed = TRUE)
    x <- gsub("&quot;", "\"", x, fixed = TRUE)
    x <- gsub("&#39;", "'", x, fixed = TRUE)
    x <- gsub("&nbsp;", " ", x, fixed = TRUE)
    x <- gsub("&amp;", "&", x, fixed = TRUE)
    x <- gsub("[\r\n\t]+", " ", x)
    x <- gsub("\\s+", " ", x)
    trimws(x)
  }
  export_result_names <- if (length(result_columns)) {
    vapply(result_columns, function(z) as.character(z$label), character(1))
  } else character()
  export_column_names <- c("Characteristic", export_result_names)
  if (has_by && test) export_column_names <- c(export_column_names, "Test p")
  if (!is.null(effect_type)) export_column_names <- c(export_column_names, paste0("Crude ", effect_type, " (95% CI)"))
  if (!is.null(effect_type) && pvalue) export_column_names <- c(export_column_names, "Crude p-value")
  if (has_adjusted) export_column_names <- c(export_column_names, paste0("Adjusted ", effect_type, " (95% CI)"))
  if (has_adjusted && pvalue) export_column_names <- c(export_column_names, "Adjusted p-value")
  if (has_multi) export_column_names <- c(export_column_names, paste0("Multivariable ", effect_type, " (95% CI)"))
  if (has_multi && pvalue) export_column_names <- c(export_column_names, "Multivariable p-value")

  export_rows <- lapply(table_rows, function(current) {
    characteristic <- if (current$type %in% c("categorical_header", "numeric_header", "numeric")) {
      clean_export_text(current$label)
    } else {
      paste0("  ", clean_export_text(current$item))
    }
    values <- if (length(result_columns)) {
      vapply(result_columns, function(column) {
        value <- if (identical(column$type, "overall")) {
          current$overall_value
        } else {
          index <- match(as.character(column$label), as.character(by_levels))
          if (is.na(index) || index > length(current$group_values)) NULL else current$group_values[[index]]
        }
        if (is.null(value) || current$type %in% c("categorical_header", "numeric_header")) "" else compact_value(value, current$kind)
      }, character(1))
    } else character()

    output_row <- c(characteristic, values)
    if (has_by && test) {
      test_text <- clean_export_text(format_p(current$test_p))
      if (nzchar(test_text) && test_note && !is.null(current$test_method) && length(test_letters)) {
        letter <- unname(test_letters[current$test_method])
        if (length(letter) && !is.na(letter) && nzchar(letter)) test_text <- paste0(test_text, " (", letter, ")")
      }
      output_row <- c(output_row, test_text)
    }
    if (!is.null(effect_type)) output_row <- c(output_row, clean_export_text(current$crude))
    if (!is.null(effect_type) && pvalue) output_row <- c(output_row, clean_export_text(format_p(current$crude_p)))
    if (has_adjusted) output_row <- c(output_row, clean_export_text(current$adjusted))
    if (has_adjusted && pvalue) output_row <- c(output_row, clean_export_text(format_p(current$adjusted_p)))
    if (has_multi) output_row <- c(output_row, clean_export_text(current$multi))
    if (has_multi && pvalue) output_row <- c(output_row, clean_export_text(format_p(current$multi_p)))
    output_row
  })

  if (length(export_rows)) {
    table_df <- as.data.frame(do.call(rbind, export_rows), stringsAsFactors = FALSE, check.names = FALSE)
  } else {
    table_df <- as.data.frame(matrix(character(), nrow = 0L, ncol = length(export_column_names)),
                              stringsAsFactors = FALSE, check.names = FALSE)
  }
  names(table_df) <- make.unique(export_column_names, sep = "_")
  rownames(table_df) <- NULL
  column_headers_html <- c(
    if (length(result_columns)) vapply(result_columns, function(z) z$html, character(1)) else character(),
    if (has_by && test) "Test p" else character(),
    if (!is.null(effect_type)) crude_heading else character(),
    if (!is.null(effect_type) && pvalue) "p-value" else character(),
    if (has_adjusted) adjusted_heading else character(),
    if (has_adjusted && pvalue) "p-value" else character(),
    if (has_multi) multi_heading else character(),
    if (has_multi && pvalue) "p-value" else character()
  )
  row_keys <- .r4vn_superby_row_keys(table_rows)
  output <- list(
    data = table_df, rows = table_rows, row_keys = row_keys, raw = raw_output, metadata = vars, by = by_name,
    by_levels = if (has_by) as.character(original_by_levels) else character(), effect_type = effect_type, event = event_level,
    adjusted = adjusted_meta, adjusted_all = adjust_all,
    multi = multi_meta, multi_model = if (is.null(multi_model)) NULL else multi_model$fit,
    multi_diagnostics = if (is.null(multi_model)) NULL else multi_model$diagnostics,
    descriptive = descriptive, html = document, table_html = table_block, note_html = note_html,
    column_headers_html = column_headers_html, css = css, common_css = common_css,
    blocks = blocks, file = normalizePath(file, winslash = "/", mustWork = TRUE),
    call = match.call()
  )
  class(output) <- "r4vn_tab"
  if (show) {
    viewer <- getOption("viewer")
    if (is.function(viewer)) viewer(output$file) else utils::browseURL(output$file)
  }
  invisible(output)
}

#' Print or Reopen an R4VN Table
#'
#' Opens the HTML file stored in an \code{r4vn_tab} object.
#'
#' @param x An object created by \code{tab()}.
#' @param ... Additional arguments currently ignored.
#' @return The input object, invisibly.
#' @keywords internal
#' @method print r4vn_tab
#' @export
print.r4vn_tab <- function(x, ...) {
  viewer <- getOption("viewer")
  if (is.function(viewer)) viewer(x$file) else utils::browseURL(x$file)
  invisible(x)
}

# ============================================================================
# Console/publication dispatcher (formerly zzz-tab-console.R)
# ============================================================================
# Non-breaking dispatcher for R4VN::tab().
# This block is intentionally placed after the publication-table implementation
# in the same file, so the original engine is captured before tab() is redefined.

if (exists("tab", mode = "function", inherits = FALSE)) {
  .r4vn_tab_publication_engine <- tab

  tab <- function(..., data = NULL, vars = NULL, by = NULL,
                  superby = NULL, digit = 1, p_digit = 3,
                  effect_digit = 2, missing = "ifany", row = FALSE,
                  col = TRUE, cell = FALSE, overall = "first",
                  descriptive = TRUE, rvrow = NULL, rvcol = FALSE,
                  test = TRUE, pvalue = TRUE, bold_p = TRUE,
                  p_bold = 0.05, test_note = TRUE,
                  interaction = TRUE, or = FALSE, rr = FALSE,
                  pr = FALSE, event = NULL, adjusted = NULL,
                  multi = NULL, effect_ref = NULL,
                  template = c("journal", "clean", "minimal"),
                  append = NULL, file = NULL, raw = FALSE,
                  name = FALSE, title = NULL, show = TRUE,
                  mode = c("auto", "console", "table")) {
    mode <- match.arg(mode)
    mc <- match.call(expand.dots = FALSE)
    dots <- as.list(mc$...)
    dot_names <- names(dots)
    if (is.null(dot_names)) dot_names <- rep("", length(dots))
    dot_names[is.na(dot_names)] <- ""
    unnamed <- which(dot_names == "")

    publication <- mode == "table" || !base::missing(vars)
    if (mode == "auto" && !publication && length(unnamed)) {
      candidates <- lapply(unnamed, function(i) {
        tryCatch(eval(dots[[i]], envir = parent.frame()),
                 error = function(e) NULL)
      })
      publication <- any(vapply(candidates, inherits, logical(1),
                                what = "r4vn_vars"))
    }

    if (publication) {
      publication_names <- c(
        "data", "vars", "by", "superby", "digit", "p_digit",
        "effect_digit", "missing", "row", "col", "cell", "overall",
        "descriptive", "rvrow", "rvcol", "test", "pvalue", "bold_p",
        "p_bold", "test_note", "interaction", "or", "rr", "pr",
        "event", "adjusted", "multi", "effect_ref", "template",
        "append", "file", "raw", "name", "title", "show"
      )
      args <- list()
      supplied <- intersect(publication_names, names(mc))
      for (nm in supplied) args[[nm]] <- mc[[nm]]

      consumed <- integer()
      vars_index <- integer()
      if (!"vars" %in% supplied && length(unnamed)) {
        for (i in unnamed) {
          value <- tryCatch(eval(dots[[i]], envir = parent.frame()),
                            error = function(e) NULL)
          if (inherits(value, "r4vn_vars")) {
            vars_index <- i
            args$vars <- dots[[i]]
            consumed <- c(consumed, i)
            break
          }
        }
      }

      if (!"data" %in% supplied) {
        before_vars <- if (length(vars_index)) {
          unnamed[unnamed < vars_index]
        } else unnamed
        before_vars <- setdiff(before_vars, consumed)
        if (length(before_vars)) {
          args$data <- dots[[before_vars[1L]]]
          consumed <- c(consumed, before_vars[1L])
        }
      }

      if (!"vars" %in% supplied && is.null(args$vars)) {
        remaining <- setdiff(unnamed, consumed)
        if (length(remaining)) {
          args$vars <- dots[[remaining[1L]]]
          consumed <- c(consumed, remaining[1L])
        }
      }

      if (!"by" %in% supplied) {
        remaining <- setdiff(unnamed, consumed)
        if (length(remaining)) {
          args$by <- dots[[remaining[1L]]]
          consumed <- c(consumed, remaining[1L])
        }
      }

      remaining_named <- which(nzchar(dot_names) &
                                 !seq_along(dots) %in% consumed)
      for (i in remaining_named) args[[dot_names[i]]] <- dots[[i]]
      remaining_unnamed <- setdiff(unnamed, consumed)
      if (length(remaining_unnamed)) {
        stop("Too many unnamed arguments for publication `tab()`.",
             call. = FALSE)
      }

      # Normalize percentage mode at the dispatcher level as well.
      # This makes the public interface intuitive even though `col = TRUE`
      # is the publication default:
      #   row = TRUE  -> col = FALSE, cell = FALSE
      #   cell = TRUE -> row = FALSE, col = FALSE
      # Otherwise the default remains column percentages.
      row_requested <- "row" %in% names(mc) &&
        isTRUE(eval(mc$row, envir = parent.frame()))
      cell_requested <- "cell" %in% names(mc) &&
        isTRUE(eval(mc$cell, envir = parent.frame()))

      if (row_requested) {
        args$row <- TRUE
        args$col <- FALSE
        args$cell <- FALSE
      } else if (cell_requested) {
        args$row <- FALSE
        args$col <- FALSE
        args$cell <- TRUE
      }

      call <- as.call(c(list(.r4vn_tab_publication_engine), args))
      return(eval(call, envir = parent.frame()))
    }

    if (!length(unnamed)) {
      stop("Console `tab()` requires a row variable.", call. = FALSE)
    }
    if (length(unnamed) > 2L) {
      stop("Console `tab()` accepts at most two unnamed variables.",
           call. = FALSE)
    }

    row_expr <- dots[[unnamed[1L]]]
    col_expr <- if (length(unnamed) >= 2L) {
      dots[[unnamed[2L]]]
    } else if ("by" %in% names(mc)) {
      mc$by
    } else NULL

    opts <- list()
    named_dots <- which(nzchar(dot_names))
    for (i in named_dots) {
      opts[[dot_names[i]]] <- eval(dots[[i]], envir = parent.frame())
    }

    # These options are shared by the publication and console interfaces.
    # Only explicitly supplied values are forwarded, so publication defaults
    # such as col = TRUE do not alter console defaults.
    shared <- c("row", "col", "cell", "missing", "show")
    for (nm in intersect(shared, names(mc))) {
      opts[[nm]] <- eval(mc[[nm]], envir = parent.frame())
    }
    if ("digit" %in% names(mc)) {
      opts$digits <- eval(mc$digit, envir = parent.frame())
    }
    if ("p_digit" %in% names(mc)) {
      opts$p_digits <- eval(mc$p_digit, envir = parent.frame())
    }

    allowed <- c(
      "percent", "row", "col", "cell", "total", "exp", "chi",
      "fisher", "lr", "residual", "adjresidual", "correct",
      "missing", "digits", "p_digits", "workspace", "show"
    )
    bad <- setdiff(names(opts), allowed)
    if (length(bad)) {
      stop(
        sprintf("Unknown console tab option(s): %s.",
                paste(bad, collapse = ", ")),
        call. = FALSE
      )
    }

    caller_env <- parent.frame()

    d <- if ("data" %in% names(mc)) {
      eval(mc$data, envir = caller_env)
    } else NULL

    # Resolve explicit data or the active data frame before do.call(). Bare
    # variable names must not be evaluated in the caller environment first.
    analysis_data <- .r4vn_stat_data(d)

    resolve_console_variable <- function(expr, arg) {
      if (is.null(expr) || identical(expr, quote(NULL))) return(NULL)

      # gioi -> "gioi"; .r4vn_eval_var() then retrieves the column from data.
      if (is.symbol(expr)) return(as.character(expr))

      # Explicit character column name.
      if (is.character(expr) && length(expr) == 1L && expr %in% names(analysis_data)) {
        return(expr)
      }

      # More complex expressions such as ivf$gioi are resolved here.
      .r4vn_eval_var(
        expr,
        data = analysis_data,
        env = caller_env,
        arg = arg
      )
    }

    row_arg <- resolve_console_variable(row_expr, "row variable")
    col_arg <- resolve_console_variable(col_expr, "column variable")

    do.call(
      .r4vn_tab_data,
      c(
        list(
          row_expr = row_arg,
          col_expr = col_arg,
          data = analysis_data,
          env = caller_env
        ),
        opts
      )
    )
  }
}