R/zzz-r4vn-inference.R

In R4VN: Health Data Analysis and Publication-Ready Reporting

# ==========================================================================
# R4VN 1.4+ inference enhancements
# ==========================================================================
# This file is deliberately loaded late. It wraps the established R4VN 1.4
# implementations so existing calls retain their behaviour while newer
# post-hoc, effect-size, and hierarchical-by features are added consistently.

.r4vn_ttest_legacy <- ttest
.r4vn_anova_legacy <- anova
.r4vn_anovai_legacy <- anovai
.r4vn_kwallis_legacy <- kwallis
.r4vn_prtest_legacy <- prtest
.r4vn_sdtest_legacy <- sdtest
.r4vn_ztest_legacy <- ztest
.r4vn_esize_legacy <- esize
.r4vn_ranksum_legacy <- ranksum

.r4vn_call_preserve <- function(fun, call, env, drop = character()) {
  cl <- call
  cl[[1L]] <- quote(.r4vn_wrapped_function)
  aa <- as.list(cl)
  if (length(drop)) aa[intersect(names(aa), drop)] <- NULL
  cl <- as.call(aa)
  ee <- new.env(parent = env)
  ee$.r4vn_wrapped_function <- fun
  eval(cl, envir = ee)
}

.r4vn_t_effect <- function(x, y = NULL, group = NULL, mu = 0, paired = FALSE,
                           level = 0.95, digits = 3) {
  alpha <- 1 - level
  zcrit <- stats::qnorm(1 - alpha / 2)
  rows <- list()
  add <- function(measure, estimate, lower = NA_real_, upper = NA_real_) {
    rows[[length(rows) + 1L]] <<- data.frame(
      Measure = measure,
      Estimate = .r4vn_num(estimate, digits),
      Lower = .r4vn_num(lower, digits),
      Upper = .r4vn_num(upper, digits),
      stringsAsFactors = FALSE
    )
  }

  if (!is.null(group)) {
    ok <- is.finite(x) & !is.na(group)
    s <- split(x[ok], droplevels(factor(group[ok])), drop = TRUE)
    if (length(s) != 2L || any(lengths(s) < 2L)) return(NULL)
    n1 <- length(s[[1L]]); n2 <- length(s[[2L]])
    m1 <- mean(s[[1L]]); m2 <- mean(s[[2L]])
    sd1 <- stats::sd(s[[1L]]); sd2 <- stats::sd(s[[2L]])
    df <- n1 + n2 - 2
    sp <- sqrt(((n1 - 1) * sd1^2 + (n2 - 1) * sd2^2) / df)
    if (!is.finite(sp) || sp <= 0) return(NULL)
    d <- (m1 - m2) / sp
    se_d <- sqrt((n1 + n2) / (n1 * n2) + d^2 / (2 * df))
    ci_d <- d + c(-1, 1) * zcrit * se_d
    J <- 1 - 3 / (4 * df - 1)
    add("Cohen d", d, ci_d[1L], ci_d[2L])
    add("Hedges g", J * d, J * ci_d[1L], J * ci_d[2L])
    if (is.finite(sd2) && sd2 > 0) add(paste0("Glass delta (SD: ", names(s)[2L], ")"), (m1 - m2) / sd2)
  } else if (!is.null(y)) {
    if (isTRUE(paired)) {
      ok <- is.finite(x) & is.finite(y)
      dxy <- x[ok] - y[ok]
      n <- length(dxy)
      if (n < 2L || !is.finite(stats::sd(dxy)) || stats::sd(dxy) <= 0) return(NULL)
      d <- mean(dxy) / stats::sd(dxy)
      df <- n - 1
      se_d <- sqrt(1 / n + d^2 / (2 * df))
      ci_d <- d + c(-1, 1) * zcrit * se_d
      J <- 1 - 3 / (4 * df - 1)
      add("Cohen dz (paired)", d, ci_d[1L], ci_d[2L])
      add("Hedges gz (paired)", J * d, J * ci_d[1L], J * ci_d[2L])
    } else {
      x <- x[is.finite(x)]; y <- y[is.finite(y)]
      if (length(x) < 2L || length(y) < 2L) return(NULL)
      return(.r4vn_t_effect(c(x, y), group = factor(c(rep("x", length(x)), rep("y", length(y)))), level = level, digits = digits))
    }
  } else {
    x <- x[is.finite(x)]
    n <- length(x)
    sx <- stats::sd(x)
    if (n < 2L || !is.finite(sx) || sx <= 0) return(NULL)
    d <- (mean(x) - mu) / sx
    df <- n - 1
    se_d <- sqrt(1 / n + d^2 / (2 * df))
    ci_d <- d + c(-1, 1) * zcrit * se_d
    J <- 1 - 3 / (4 * df - 1)
    add("Cohen d", d, ci_d[1L], ci_d[2L])
    add("Hedges g", J * d, J * ci_d[1L], J * ci_d[2L])
  }
  if (!length(rows)) NULL else do.call(rbind, rows)
}

.r4vn_adjustment_label <- function(method) {
  key <- as.character(method)[1L]
  switch(
    tolower(key),
    bonferroni = "Bonferroni",
    holm = "Holm",
    hochberg = "Hochberg",
    hommel = "Hommel",
    bh = "Benjamini-Hochberg",
    by = "Benjamini-Yekutieli",
    fdr = "false-discovery-rate",
    none = "no multiplicity",
    key
  )
}

.r4vn_posthoc_section_title <- function(method, adjust = "holm") {
  adjustment <- .r4vn_adjustment_label(adjust)
  switch(
    method,
    tukey = "Tukey honestly significant difference test",
    `games-howell` = "Games-Howell multiple comparisons",
    scheffe = "Scheffe multiple comparisons",
    bonferroni = "Bonferroni-adjusted pairwise t-tests",
    pairwise = paste0("Pairwise t-tests (", adjustment, "-adjusted p-values)"),
    dunn = paste0("Dunn's test (", adjustment, "-adjusted p-values)"),
    wilcoxon = paste0("Pairwise Wilcoxon rank-sum tests (", adjustment,
                      "-adjusted p-values)"),
    method
  )
}

.r4vn_posthoc_summary <- function(mat, group.names = NULL,
                                  method = c("tukey", "games-howell", "scheffe",
                                             "bonferroni", "pairwise"),
                                  adjust = "holm", level = 0.95,
                                  digits = 3, p_digits = 3) {
  method <- match.arg(method)
  mat <- as.matrix(mat)
  if (nrow(mat) < 2L) return(NULL)
  colnames(mat) <- c("n", "mean", "sd")
  gn <- group.names %||% rownames(mat) %||% paste0("Group ", seq_len(nrow(mat)))
  k <- nrow(mat); N <- sum(mat[, "n"]); dfw <- N - k
  mse <- sum((mat[, "n"] - 1) * mat[, "sd"]^2) / dfw
  prs <- utils::combn(seq_len(k), 2L, simplify = FALSE)
  alpha <- 1 - level
  raw_p <- numeric(length(prs)); out <- vector("list", length(prs))

  for (ii in seq_along(prs)) {
    i <- prs[[ii]][1L]; j <- prs[[ii]][2L]
    ni <- mat[i, "n"]; nj <- mat[j, "n"]
    mi <- mat[i, "mean"]; mj <- mat[j, "mean"]
    si <- mat[i, "sd"]; sj <- mat[j, "sd"]
    dif <- mi - mj
    if (method == "games-howell") {
      se <- sqrt(si^2 / ni + sj^2 / nj)
      df <- (si^2 / ni + sj^2 / nj)^2 / ((si^2 / ni)^2 / (ni - 1) + (sj^2 / nj)^2 / (nj - 1))
      stat <- dif / se
      qv <- sqrt(2) * abs(stat)
      pp <- stats::ptukey(qv, nmeans = k, df = df, lower.tail = FALSE)
      crit <- stats::qtukey(1 - alpha, nmeans = k, df = df) / sqrt(2)
      ci <- dif + c(-1, 1) * crit * se
      stat_name <- "t (Welch)"
    } else {
      se <- sqrt(mse * (1 / ni + 1 / nj))
      stat <- dif / se
      df <- dfw
      if (method == "tukey") {
        qv <- sqrt(2) * abs(stat)
        pp <- stats::ptukey(qv, nmeans = k, df = dfw, lower.tail = FALSE)
        crit <- stats::qtukey(1 - alpha, nmeans = k, df = dfw) / sqrt(2)
        ci <- dif + c(-1, 1) * crit * se
        stat_name <- "t"
      } else if (method == "scheffe") {
        fv <- stat^2 / (k - 1)
        pp <- stats::pf(fv, k - 1, dfw, lower.tail = FALSE)
        crit <- sqrt((k - 1) * stats::qf(1 - alpha, k - 1, dfw))
        ci <- dif + c(-1, 1) * crit * se
        stat_name <- "t"
      } else {
        pp <- 2 * stats::pt(abs(stat), dfw, lower.tail = FALSE)
        adjust_used <- if (method == "bonferroni") "bonferroni" else adjust
        ci_alpha <- if (tolower(adjust_used) == "bonferroni") {
          alpha / length(prs)
        } else {
          alpha
        }
        crit <- stats::qt(1 - ci_alpha / 2, dfw)
        ci <- dif + c(-1, 1) * crit * se
        stat_name <- "t"
      }
    }
    raw_p[ii] <- pp
    out[[ii]] <- data.frame(
      Comparison = paste0(gn[i], " - ", gn[j]),
      Difference = .r4vn_num(dif, digits),
      SE = .r4vn_num(se, digits),
      Statistic = .r4vn_num(stat, digits),
      df = .r4vn_num(df, 2),
      Lower = .r4vn_num(ci[1L], digits),
      Upper = .r4vn_num(ci[2L], digits),
      p = .r4vn_p(pp, p_digits),
      stringsAsFactors = FALSE,
      check.names = FALSE
    )
  }
  if (method %in% c("bonferroni", "pairwise")) {
    adjust_used <- if (method == "bonferroni") "bonferroni" else adjust
    padj <- stats::p.adjust(raw_p, method = adjust_used)
    for (ii in seq_along(out)) out[[ii]]$p.adjusted <- .r4vn_p(padj[ii], p_digits)
  }
  ans <- do.call(rbind, out); rownames(ans) <- NULL
  attr(ans, "method") <- method
  attr(ans, "adjust") <- if (method %in% c("bonferroni", "pairwise")) adjust_used else NULL
  ans
}

.r4vn_kw_posthoc <- function(x, g, method = c("dunn", "wilcoxon"), adjust = "holm",
                             digits = 3, p_digits = 3) {
  method <- match.arg(method)
  ok <- is.finite(x) & !is.na(g); x <- x[ok]; g <- droplevels(factor(g[ok]))
  prs <- utils::combn(levels(g), 2L, simplify = FALSE)
  if (!length(prs)) return(NULL)
  if (method == "wilcoxon") {
    pp <- vapply(prs, function(z) {
      tryCatch(stats::wilcox.test(x[g == z[1L]], x[g == z[2L]], exact = FALSE)$p.value, error = function(e) NA_real_)
    }, numeric(1))
    pa <- stats::p.adjust(pp, method = adjust)
    return(data.frame(
      Comparison = vapply(prs, paste, collapse = " - ", character(1)),
      p = .r4vn_p(pp, p_digits), p.adjusted = .r4vn_p(pa, p_digits),
      Method = "Pairwise Wilcoxon rank-sum test", stringsAsFactors = FALSE
    ))
  }

  r <- rank(x, ties.method = "average")
  N <- length(r)
  tabtie <- table(x)
  C <- if (N > 1L) 1 - sum(tabtie^3 - tabtie) / (N^3 - N) else 1
  means <- tapply(r, g, mean); ns <- table(g)
  zval <- vapply(prs, function(z) {
    den <- sqrt((N * (N + 1) / 12) * C * (1 / ns[[z[1L]]] + 1 / ns[[z[2L]]]))
    (means[[z[1L]]] - means[[z[2L]]]) / den
  }, numeric(1))
  pp <- 2 * stats::pnorm(abs(zval), lower.tail = FALSE)
  pa <- stats::p.adjust(pp, method = adjust)
  data.frame(
    Comparison = vapply(prs, paste, collapse = " - ", character(1)),
    Mean.rank.difference = .r4vn_num(vapply(prs, function(z) means[[z[1L]]] - means[[z[2L]]], numeric(1)), digits),
    z = .r4vn_num(zval, digits), p = .r4vn_p(pp, p_digits),
    p.adjusted = .r4vn_p(pa, p_digits), Method = "Dunn's test",
    stringsAsFactors = FALSE, check.names = FALSE
  )
}

.r4vn_kw_effect <- function(fit, n, k, digits = 3) {
  H <- unname(fit$statistic)
  eps <- max(0, (H - k + 1) / (n - k))
  eta <- max(0, (H - k + 1) / (n - 1))
  data.frame(Measure = c("Epsilon-squared", "Eta-squared (H)"),
             Estimate = .r4vn_num(c(eps, eta), digits), stringsAsFactors = FALSE)
}

# Student t tests with optional effect sizes and hierarchical grouping
# @inheritParams ttest
# @param effect Add standardized effect sizes (Cohen d and Hedges g). The
#   default is `FALSE` for backward compatibility.
# @export
ttest <- function(x, y = NULL, by = NULL, data = NULL, mu = 0, equal = FALSE, paired = FALSE,
                  alternative = c("two.sided", "less", "greater"), level = 0.95,
                  effect = FALSE, digits = 3, p_digits = 3, show = TRUE, console = FALSE) {
  call <- match.call(); env <- parent.frame(); d <- .r4vn_stat_data(data)
  by_expr <- substitute(by)
  spec <- .r4vn_by_spec(by_expr, d, env, allow_null = TRUE)
  xnm <- tryCatch(.r4vn_resolve_name_spec(substitute(x), d, env, "x", multiple = FALSE), error = function(e) NULL)

  if (length(spec$all) && .r4vn_is_vars_spec_expr(by_expr, d, env) && !is.null(xnm)) {
    ids <- .r4vn_strata_indices(d, spec$strata); results <- list(); labs <- character()
    for (idx in ids) {
      dd <- d[idx, , drop = FALSE]
      if (!nrow(dd)) next
      z <- do.call(.r4vn_ttest_legacy, list(x = xnm, by = spec$by, data = dd, mu = mu, equal = equal,
                  paired = FALSE, alternative = match.arg(alternative), level = level,
                  digits = digits, p_digits = p_digits, show = FALSE, console = FALSE))
      if (isTRUE(effect)) {
        et <- .r4vn_t_effect(dd[[xnm]], group = dd[[spec$by]], level = level, digits = digits)
        if (!is.null(et)) z$sections[["Effect size"]] <- et
      }
      results[[length(results) + 1L]] <- z
      labs <- c(labs, .r4vn_stratum_label(d, spec$strata, idx))
    }
    if (!length(spec$strata)) {
      out <- results[[1L]]; out$call <- call
    } else {
      out <- .r4vn_stat_collection("t test by hierarchical strata", results, labels = labs,
        notes = .r4vn_hierarchy_note(d, spec), call = call)
    }
    return(.r4vn_show(out, show = show, console = console))
  }

  call_base <- call
  if (isTRUE(effect)) call_base$show <- FALSE
  out <- .r4vn_call_preserve(.r4vn_ttest_legacy, call_base, env, drop = "effect")
  if (isTRUE(effect)) {
    xv <- .r4vn_eval_var(substitute(x), d, env, "x")
    et <- NULL
    if (length(spec$all)) {
      et <- .r4vn_t_effect(xv, group = d[[spec$by]], level = level, digits = digits)
    } else if (!missing(y) && !identical(substitute(y), quote(NULL))) {
      yv <- .r4vn_eval_var(substitute(y), d, env, "y")
      et <- .r4vn_t_effect(xv, y = yv, paired = paired, level = level, digits = digits)
    } else et <- .r4vn_t_effect(xv, mu = mu, level = level, digits = digits)
    if (!is.null(et)) out$sections[["Effect size"]] <- et
    out$call <- call
    return(.r4vn_show(out, show = show, console = console))
  }
  invisible(out)
}

# One-way ANOVA with post-hoc comparisons
# @inheritParams anova
# @param posthoc Post-hoc method: `"none"`, `"tukey"`, `"games-howell"`,
#   `"scheffe"`, `"bonferroni"`, or `"pairwise"`.
# @param adjust Multiplicity adjustment used by `posthoc = "pairwise"`.
# @export
anova <- function(object, ..., by = NULL, data = NULL, bartlett = TRUE, level = 0.95,
                  posthoc = c("none", "tukey", "games-howell", "scheffe",
                              "bonferroni", "pairwise"),
                  adjust = "holm", digits = 3, p_digits = 3, show = TRUE, console = FALSE) {
  call <- match.call(); env <- parent.frame(); posthoc <- match.arg(posthoc)
  # Preserve stats::anova-style model comparisons exactly.
  d <- tryCatch(.r4vn_stat_data(data), error = function(e) NULL)
  if (is.null(d) || missing(by) || identical(substitute(by), quote(NULL))) {
    return(.r4vn_call_preserve(.r4vn_anova_legacy, call, env, drop = c("posthoc", "adjust")))
  }
  spec <- .r4vn_by_spec(substitute(by), d, env, allow_null = FALSE)
  xnm <- .r4vn_resolve_name_spec(substitute(object), d, env, "object", multiple = FALSE)
  ids <- .r4vn_strata_indices(d, spec$strata)
  results <- list(); labs <- character()
  for (idx in ids) {
    dd <- d[idx, , drop = FALSE]
    z <- do.call(.r4vn_anova_legacy, list(object = xnm, by = spec$by, data = dd, bartlett = bartlett,
                  level = level, digits = digits, p_digits = p_digits, show = FALSE, console = FALSE))
    if (posthoc != "none" && !is.null(z$raw$groups)) {
      section_title <- .r4vn_posthoc_section_title(
        posthoc, if (posthoc == "bonferroni") "bonferroni" else adjust
      )
      z$sections[[section_title]] <- .r4vn_posthoc_summary(z$raw$groups,
        rownames(z$raw$groups) %||% z$sections[["Group statistics"]]$Group,
        method = posthoc, adjust = adjust, level = level, digits = digits, p_digits = p_digits)
    }
    results[[length(results) + 1L]] <- z
    labs <- c(labs, .r4vn_stratum_label(d, spec$strata, idx))
  }
  if (!length(spec$strata)) {
    out <- results[[1L]]; out$call <- call
  } else {
    out <- .r4vn_stat_collection("One-way ANOVA by hierarchical strata", results, labels = labs,
      notes = .r4vn_hierarchy_note(d, spec), call = call)
  }
  .r4vn_show(out, show = show, console = console)
}

# Immediate one-way ANOVA with optional post-hoc comparisons
# @inheritParams anovai
# @param posthoc Post-hoc method: `"none"`, `"tukey"`, `"games-howell"`,
#   `"scheffe"`, `"bonferroni"`, or `"pairwise"`.
# @param adjust Multiplicity adjustment used by `posthoc = "pairwise"`.
# @export
anovai <- function(..., group.names = NULL, bartlett = TRUE, level = 0.95,
                   posthoc = c("none", "tukey", "games-howell", "scheffe",
                               "bonferroni", "pairwise"),
                   adjust = "holm", digits = 3, p_digits = 3, show = TRUE, console = FALSE) {
  posthoc <- match.arg(posthoc); call <- match.call()
  groups <- list(...)
  out <- do.call(.r4vn_anovai_legacy, c(groups, list(group.names = group.names, bartlett = bartlett,
        level = level, digits = digits, p_digits = p_digits, show = FALSE, console = FALSE)))
  if (posthoc != "none") {
    gn <- group.names %||% out$sections[["Group statistics"]]$Group
    section_title <- .r4vn_posthoc_section_title(
      posthoc, if (posthoc == "bonferroni") "bonferroni" else adjust
    )
    out$sections[[section_title]] <- .r4vn_posthoc_summary(out$raw$groups, gn,
      method = posthoc, adjust = adjust, level = level, digits = digits, p_digits = p_digits)
  }
  out$call <- call
  .r4vn_show(out, show = show, console = console)
}

# Kruskal-Wallis test with post-hoc comparisons and effect sizes
# @inheritParams kwallis
# @param posthoc `"none"`, `"dunn"`, or `"wilcoxon"`.
# @param adjust Multiplicity adjustment passed to `p.adjust()`.
# @param effect Add epsilon-squared and eta-squared(H) effect sizes.
# @export
kwallis <- function(x, by, data = NULL, posthoc = c("none", "dunn", "wilcoxon"),
                    adjust = "holm", effect = FALSE, digits = 3, p_digits = 3,
                    show = TRUE, console = FALSE) {
  call <- match.call(); env <- parent.frame(); posthoc <- match.arg(posthoc); d <- .r4vn_stat_data(data)
  spec <- .r4vn_by_spec(substitute(by), d, env, allow_null = FALSE)
  xnm <- .r4vn_resolve_name_spec(substitute(x), d, env, "x", multiple = FALSE)
  ids <- .r4vn_strata_indices(d, spec$strata); results <- list(); labs <- character()
  for (idx in ids) {
    dd <- d[idx, , drop = FALSE]
    z <- do.call(.r4vn_kwallis_legacy, list(x = xnm, by = spec$by, data = dd,
                 digits = digits, p_digits = p_digits, show = FALSE, console = FALSE))
    ok <- is.finite(dd[[xnm]]) & !is.na(dd[[spec$by]])
    xx <- dd[[xnm]][ok]; gg <- droplevels(factor(dd[[spec$by]][ok]))
    if (posthoc != "none") {
      section_title <- .r4vn_posthoc_section_title(posthoc, adjust)
      z$sections[[section_title]] <- .r4vn_kw_posthoc(
        xx, gg, posthoc, adjust, digits, p_digits
      )
    }
    if (isTRUE(effect)) z$sections[["Effect size"]] <- .r4vn_kw_effect(z$raw$test, length(xx), nlevels(gg), digits)
    results[[length(results) + 1L]] <- z; labs <- c(labs, .r4vn_stratum_label(d, spec$strata, idx))
  }
  if (!length(spec$strata)) { out <- results[[1L]]; out$call <- call }
  else out <- .r4vn_stat_collection("Kruskal-Wallis test by hierarchical strata", results, labels = labs,
    notes = .r4vn_hierarchy_note(d, spec), call = call)
  .r4vn_show(out, show = show, console = console)
}

# prtest already existed in R4VN 1.4. This wrapper extends its by convention
# without changing one-sample p0 or ordinary two-group calls.
# @export
prtest <- function(x, by = NULL, data = NULL, event = NULL, p0 = NULL,
                   alternative = c("two.sided", "less", "greater"), level = 0.95,
                   digits = 3, p_digits = 3, show = TRUE, console = FALSE) {
  call <- match.call(); env <- parent.frame(); d <- .r4vn_stat_data(data)
  spec <- .r4vn_by_spec(substitute(by), d, env, allow_null = TRUE)
  if (!length(spec$strata) && !.r4vn_is_vars_spec_expr(substitute(by), d, env)) return(.r4vn_call_preserve(.r4vn_prtest_legacy, call, env))
  if (!length(spec$all)) return(.r4vn_call_preserve(.r4vn_prtest_legacy, call, env))
  xnm <- .r4vn_resolve_name_spec(substitute(x), d, env, "x", multiple = FALSE)
  ids <- .r4vn_strata_indices(d, spec$strata); results <- list(); labs <- character()
  for (idx in ids) {
    dd <- d[idx, , drop = FALSE]
    z <- do.call(.r4vn_prtest_legacy, list(x = xnm, by = spec$by, data = dd, event = event,
              p0 = p0, alternative = match.arg(alternative), level = level, digits = digits,
              p_digits = p_digits, show = FALSE, console = FALSE))
    results[[length(results) + 1L]] <- z; labs <- c(labs, .r4vn_stratum_label(d, spec$strata, idx))
  }
  if (!length(spec$strata)) {
    out <- results[[1L]]; out$call <- call
  } else {
    out <- .r4vn_stat_collection("Proportion z test by hierarchical strata", results, labels = labs,
      notes = .r4vn_hierarchy_note(d, spec), call = call)
  }
  .r4vn_show(out, show = show, console = console)
}


# Additional two-group commands share the same hierarchical by = vars(...)
# convention. Calls with ordinary `by = group` are delegated unchanged to the
# established R4VN 1.4 implementations.
.r4vn_hier_two_group <- function(fun, title, x_expr, by_expr, data, env, args, call,
                                 show = TRUE, console = FALSE) {
  spec <- .r4vn_by_spec(by_expr, data, env, allow_null = FALSE)
  xnm <- .r4vn_resolve_name_spec(x_expr, data, env, "x", multiple = FALSE)
  ids <- .r4vn_strata_indices(data, spec$strata)
  results <- list(); labs <- character()
  for (idx in ids) {
    dd <- data[idx, , drop = FALSE]
    zz <- c(list(x = xnm, by = spec$by, data = dd), args,
            list(show = FALSE, console = FALSE))
    z <- do.call(fun, zz)
    results[[length(results) + 1L]] <- z
    labs <- c(labs, .r4vn_stratum_label(data, spec$strata, idx))
  }
  if (!length(results)) stop("No complete strata are available for analysis.", call. = FALSE)
  if (!length(spec$strata)) {
    out <- results[[1L]]; out$call <- call
  } else {
    out <- .r4vn_stat_collection(title, results, labels = labs,
      notes = .r4vn_hierarchy_note(data, spec), call = call)
  }
  .r4vn_show(out, show = show, console = console)
}

sdtest <- function(x, y = NULL, by = NULL, data = NULL, sd0 = NULL,
                   alternative = c("two.sided", "less", "greater"), level = 0.95,
                   digits = 3, p_digits = 3, show = TRUE, console = FALSE) {
  call <- match.call(); env <- parent.frame(); bx <- substitute(by)
  d <- .r4vn_stat_data(data)
  if (!.r4vn_is_vars_spec_expr(bx, d, env)) return(.r4vn_call_preserve(.r4vn_sdtest_legacy, call, env))
  .r4vn_hier_two_group(.r4vn_sdtest_legacy, "Variance comparison by hierarchical strata",
    substitute(x), bx, d, env,
    list(sd0 = sd0, alternative = match.arg(alternative), level = level,
         digits = digits, p_digits = p_digits), call, show, console)
}

esize <- function(x, by, data = NULL, level = 0.95, digits = 3,
                  show = TRUE, console = FALSE) {
  call <- match.call(); env <- parent.frame(); bx <- substitute(by)
  d <- .r4vn_stat_data(data)
  if (!.r4vn_is_vars_spec_expr(bx, d, env)) return(.r4vn_call_preserve(.r4vn_esize_legacy, call, env))
  .r4vn_hier_two_group(.r4vn_esize_legacy, "Effect size by hierarchical strata",
    substitute(x), bx, d, env,
    list(level = level, digits = digits), call, show, console)
}

ztest <- function(x, sigma1, y = NULL, sigma2 = NULL, by = NULL, data = NULL, mu = 0,
                  alternative = c("two.sided", "less", "greater"), level = 0.95,
                  digits = 3, p_digits = 3, show = TRUE, console = FALSE) {
  call <- match.call(); env <- parent.frame(); bx <- substitute(by)
  d <- .r4vn_stat_data(data)
  if (!.r4vn_is_vars_spec_expr(bx, d, env)) return(.r4vn_call_preserve(.r4vn_ztest_legacy, call, env))
  if (is.null(sigma2)) stop("`sigma2` is required for grouped z tests.", call. = FALSE)
  .r4vn_hier_two_group(.r4vn_ztest_legacy, "z test by hierarchical strata",
    substitute(x), bx, d, env,
    list(sigma1 = sigma1, sigma2 = sigma2, mu = mu,
         alternative = match.arg(alternative), level = level,
         digits = digits, p_digits = p_digits), call, show, console)
}

ranksum <- function(x, y = NULL, by = NULL, data = NULL,
                    alternative = c("two.sided", "less", "greater"),
                    exact = NULL, correct = TRUE, conf.int = TRUE,
                    level = 0.95, digits = 3, p_digits = 3,
                    show = TRUE, console = FALSE) {
  call <- match.call(); env <- parent.frame(); bx <- substitute(by)
  d <- .r4vn_stat_data(data)
  if (!.r4vn_is_vars_spec_expr(bx, d, env)) return(.r4vn_call_preserve(.r4vn_ranksum_legacy, call, env))
  .r4vn_hier_two_group(.r4vn_ranksum_legacy, "Wilcoxon rank-sum test by hierarchical strata",
    substitute(x), bx, d, env,
    list(alternative = match.arg(alternative), exact = exact, correct = correct,
         conf.int = conf.int, level = level, digits = digits, p_digits = p_digits),
    call, show, console)
}