R/nonparametric.R

In R4VN: Health Data Analysis and Publication-Ready Reporting

# R4VN standalone nonparametric statistical commands
# Base R only. These functions use the existing r4vn_stat result and printing
# infrastructure defined in statistics.R/statistics-utils.R.

.r4vn_np_desc <- function(x, group = "All", digits = 3L) {
  x <- x[is.finite(x)]
  if (!length(x)) {
    return(data.frame(Group = group, n = 0L, Median = "", Q1 = "", Q3 = "",
                      IQR = "", Min = "", Max = "", stringsAsFactors = FALSE))
  }
  q <- stats::quantile(x, c(0, .25, .5, .75, 1), names = FALSE, type = 2)
  data.frame(
    Group = group,
    n = length(x),
    Median = .r4vn_num(q[3L], digits),
    Q1 = .r4vn_num(q[2L], digits),
    Q3 = .r4vn_num(q[4L], digits),
    IQR = .r4vn_num(q[4L] - q[2L], digits),
    Min = .r4vn_num(q[1L], digits),
    Max = .r4vn_num(q[5L], digits),
    stringsAsFactors = FALSE,
    check.names = FALSE
  )
}

.r4vn_np_wilcox <- function(args) {
  notes <- character()
  fit <- withCallingHandlers(
    tryCatch(
      do.call(stats::wilcox.test, args),
      error = function(e) {
        if (isTRUE(args$conf.int)) {
          notes <<- c(notes, paste0("Confidence interval was not available: ", conditionMessage(e)))
          args$conf.int <- FALSE
          return(do.call(stats::wilcox.test, args))
        }
        stop(conditionMessage(e), call. = FALSE)
      }
    ),
    warning = function(w) {
      notes <<- c(notes, conditionMessage(w))
      invokeRestart("muffleWarning")
    }
  )
  list(fit = fit, notes = unique(notes))
}

.r4vn_np_wilcox_table <- function(fit, digits = 3L, p_digits = 3L) {
  est <- if (is.null(fit$estimate)) NA_real_ else unname(fit$estimate[1L])
  ci <- if (is.null(fit$conf.int)) c(NA_real_, NA_real_) else unname(fit$conf.int[1:2])
  stat <- if (is.null(fit$statistic)) NA_real_ else unname(fit$statistic[1L])
  stat_name <- names(fit$statistic)
  if (is.null(stat_name) || !length(stat_name) || is.na(stat_name[1L]) || !nzchar(stat_name[1L])) stat_name <- "W"
  data.frame(
    Statistic = stat_name[1L],
    Value = .r4vn_num(stat, digits),
    Estimate = .r4vn_num(est, digits),
    Lower = .r4vn_num(ci[1L], digits),
    Upper = .r4vn_num(ci[2L], digits),
    p = .r4vn_p(fit$p.value, p_digits),
    Alternative = fit$alternative,
    stringsAsFactors = FALSE,
    check.names = FALSE
  )
}

#' Wilcoxon Rank-sum Test
#'
#' Performs the Wilcoxon rank-sum test, also known as the Mann-Whitney U test,
#' for two independent samples. Samples may be supplied as two numeric
#' variables or as one numeric outcome and a two-level grouping variable.
#'
#' @param x Numeric outcome or first numeric sample.
#' @param y Optional second numeric sample.
#' @param by Optional two-level grouping variable. Use either `y` or `by`.
#' @param data Data frame. If `NULL`, the active R4VN data frame is used.
#' @param alternative Alternative hypothesis: `"two.sided"`, `"less"`, or
#'   `"greater"`.
#' @param exact Use an exact p-value when possible. `NULL` lets R decide.
#' @param correct Apply continuity correction for the normal approximation.
#' @param conf.int Report the Hodges-Lehmann location-shift estimate and its
#'   confidence interval when available.
#' @param level Confidence level.
#' @param digits,p_digits Decimal places for estimates and p-values.
#' @param show Logical; open the formatted result in the Viewer. Default `TRUE`.
#' @param console Logical; also print the traditional result in the Console. Default `FALSE`.
#'
#' @details
#' With `by`, the first observed factor level is sample 1 and the second level is
#' sample 2. Use `factor()` or `labvar(..., ref = ...)` to control level order.
#' Missing values are removed independently when `x` and `y` are supplied, and
#' complete cases are used when `by` is supplied.
#'
#' @return Invisibly returns an object of class `r4vn_stat`.
#' @export
#'
#' @examples
#' d <- data.frame(
#'   score = c(12, 15, 11, 19, 18, 21, 14, 17),
#'   group = factor(rep(c("Control", "Intervention"), each = 4)),
#'   score2 = c(10, 13, 12, 14, 19, 20, 18, 22)
#' )
#' ranksum(score, by = group, data = d)
#' ranksum(score, score2, data = d)
#' ranksum(score, by = group, data = d, alternative = "less")
#'
#' # Extended usage examples
#' \donttest{
#' d <- data.frame(
#'   score = c(10, 11, 12, 13, 18, 19, 20, 21),
#'   score2 = c(9, 10, 12, 11, 17, 18, 19, 22),
#'   group = factor(rep(c("Control", "Intervention"), each = 4))
#' )
#'
#' # Two independent variables or one outcome by a two-level group
#' ranksum(score, score2, data = d)
#' ranksum(score, by = group, data = d)
#'
#' # One-sided alternatives, approximation controls, and confidence interval
#' ranksum(score, by = group, data = d, alternative = "less")
#' ranksum(score, by = group, data = d, exact = FALSE, correct = FALSE)
#' ranksum(score, by = group, data = d, conf.int = FALSE)
#'
#' # Active data and hidden console output
#' usedf(d)
#' result <- ranksum(score, by = group, show = FALSE)
#' result$raw$test
#' }
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) {
  env <- parent.frame()
  data <- .r4vn_stat_data(data)
  alternative <- match.arg(alternative)
  x_expr <- substitute(x)
  y_expr <- substitute(y)
  by_expr <- substitute(by)
  has_y <- !missing(y) && !identical(y_expr, quote(NULL))
  has_by <- !missing(by) && !identical(by_expr, quote(NULL))
  if (has_y && has_by) stop("Use either `y` or `by`, not both.", call. = FALSE)
  if (!has_y && !has_by) stop("Supply a second sample in `y` or a two-level `by` variable.", call. = FALSE)

  xv <- .r4vn_eval_var(x_expr, data, env, "x")
  if (!is.numeric(xv)) stop("`x` must be numeric.", call. = FALSE)

  if (has_by) {
    gv <- .r4vn_eval_var(by_expr, data, env, "by")
    ok <- is.finite(xv) & !is.na(gv)
    gf <- droplevels(factor(gv[ok]))
    xv <- xv[ok]
    if (nlevels(gf) != 2L) stop("`by` must have exactly two observed groups.", call. = FALSE)
    samples <- split(xv, gf, drop = TRUE)
    x1 <- samples[[1L]]
    x2 <- samples[[2L]]
    group_names <- names(samples)
  } else {
    yv <- .r4vn_eval_var(y_expr, data, env, "y")
    if (!is.numeric(yv)) stop("`y` must be numeric.", call. = FALSE)
    x1 <- xv[is.finite(xv)]
    x2 <- yv[is.finite(yv)]
    group_names <- c(.r4vn_name(x_expr, "x"), .r4vn_name(y_expr, "y"))
  }
  if (!length(x1) || !length(x2)) stop("Both samples must contain nonmissing observations.", call. = FALSE)

  run <- .r4vn_np_wilcox(list(
    x = x1, y = x2, paired = FALSE, mu = 0,
    alternative = alternative, exact = exact, correct = correct,
    conf.int = conf.int, conf.level = level
  ))
  desc <- rbind(
    .r4vn_np_desc(x1, group_names[1L], digits),
    .r4vn_np_desc(x2, group_names[2L], digits)
  )
  raw <- list(sample1 = x1, sample2 = x2, groups = group_names, test = run$fit)
  result <- .r4vn_result(
    "Wilcoxon rank-sum test",
    list("Descriptive statistics" = desc,
         "Test" = .r4vn_np_wilcox_table(run$fit, digits, p_digits)),
    notes = if (length(run$notes)) run$notes else NULL,
    raw = raw,
    call = match.call()
  )
  .r4vn_show(result, show)
}

#' Wilcoxon Signed-rank Test
#'
#' Performs a one-sample or paired Wilcoxon signed-rank test. With two variables,
#' complete pairs are analyzed and the tested difference is `x - y`.
#'
#' @param x Numeric variable or first paired measurement.
#' @param y Optional second paired measurement.
#' @param data Data frame. If `NULL`, active R4VN data is used.
#' @param mu Null median or null median paired difference.
#' @param alternative Alternative hypothesis.
#' @param exact Use an exact p-value when possible. `NULL` lets R decide.
#' @param correct Apply continuity correction for the normal approximation.
#' @param conf.int Report a pseudomedian or paired location-shift estimate and
#'   confidence interval when available.
#' @param level Confidence level.
#' @param digits,p_digits Decimal places for estimates and p-values.
#' @param show Logical; open the formatted result in the Viewer. Default `TRUE`.
#' @param console Logical; also print the traditional result in the Console. Default `FALSE`.
#'
#' @return Invisibly returns an object of class `r4vn_stat`.
#' @export
#'
#' @examples
#' d <- data.frame(
#'   before = c(18, 20, 15, 22, 17, 19, 24, 16),
#'   after  = c(15, 18, 14, 19, 16, 17, 20, 15)
#' )
#' signrank(before, data = d, mu = 18)
#' signrank(before, after, data = d)
#' signrank(before, after, data = d, alternative = "greater")
#'
#' # Extended usage examples
#' \donttest{
#' d <- data.frame(
#'   before = c(18, 20, 15, 22, 17, 19, 24, 16),
#'   after = c(15, 18, 14, 19, 16, 17, 20, 15)
#' )
#'
#' # One-sample signed-rank test against a specified median
#' signrank(before, data = d, mu = 18)
#'
#' # Paired signed-rank test; the analyzed difference is before - after
#' signrank(before, after, data = d)
#' signrank(before, after, data = d, alternative = "greater")
#'
#' # Approximation and confidence-interval controls
#' signrank(before, after, data = d, exact = FALSE, correct = FALSE)
#' signrank(before, after, data = d, conf.int = FALSE)
#'
#' usedf(d)
#' signrank(before, after)
#' }
signrank <- function(x, y = NULL, data = NULL, mu = 0,
                     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) {
  env <- parent.frame()
  data <- .r4vn_stat_data(data)
  alternative <- match.arg(alternative)
  x_expr <- substitute(x)
  y_expr <- substitute(y)
  has_y <- !missing(y) && !identical(y_expr, quote(NULL))
  xv <- .r4vn_eval_var(x_expr, data, env, "x")
  if (!is.numeric(xv)) stop("`x` must be numeric.", call. = FALSE)

  if (has_y) {
    yv <- .r4vn_eval_var(y_expr, data, env, "y")
    if (!is.numeric(yv)) stop("`y` must be numeric.", call. = FALSE)
    ok <- is.finite(xv) & is.finite(yv)
    x1 <- xv[ok]
    x2 <- yv[ok]
    if (!length(x1)) stop("No complete pairs are available.", call. = FALSE)
    diff <- x1 - x2
    run <- .r4vn_np_wilcox(list(
      x = x1, y = x2, paired = TRUE, mu = mu,
      alternative = alternative, exact = exact, correct = correct,
      conf.int = conf.int, conf.level = level
    ))
    desc <- rbind(
      .r4vn_np_desc(x1, .r4vn_name(x_expr, "x"), digits),
      .r4vn_np_desc(x2, .r4vn_name(y_expr, "y"), digits),
      .r4vn_np_desc(diff, "Difference (x - y)", digits)
    )
    title <- "Paired Wilcoxon signed-rank test"
    raw <- list(x = x1, y = x2, difference = diff, mu = mu, test = run$fit)
  } else {
    x1 <- xv[is.finite(xv)]
    if (!length(x1)) stop("`x` has no nonmissing observations.", call. = FALSE)
    run <- .r4vn_np_wilcox(list(
      x = x1, mu = mu, paired = FALSE,
      alternative = alternative, exact = exact, correct = correct,
      conf.int = conf.int, conf.level = level
    ))
    desc <- .r4vn_np_desc(x1, .r4vn_name(x_expr, "x"), digits)
    title <- "One-sample Wilcoxon signed-rank test"
    raw <- list(x = x1, mu = mu, test = run$fit)
  }
  test_table <- .r4vn_np_wilcox_table(run$fit, digits, p_digits)
  test_table$Null <- .r4vn_num(mu, digits)
  result <- .r4vn_result(
    title,
    list("Descriptive statistics" = desc, "Test" = test_table),
    notes = if (length(run$notes)) run$notes else NULL,
    raw = raw,
    call = match.call()
  )
  .r4vn_show(result, show)
}

#' Kruskal-Wallis Test
#'
#' Performs the Kruskal-Wallis rank-sum test for a numeric outcome across two or
#' more independent groups. Optional Dunn or pairwise Wilcoxon post-hoc tests
#' and rank-based effect sizes can be requested. With `by = vars(region, sex,
#' treatment)`, region and sex are nested strata and treatment is the innermost
#' Kruskal-Wallis factor.
#'
#' @usage
#' kwallis(
#'   x, by, data = NULL, posthoc = c("none", "dunn", "wilcoxon"),
#'   adjust = "holm", effect = FALSE, digits = 3, p_digits = 3, show = TRUE,
#'   console = FALSE
#' )
#'
#' @param x Numeric outcome.
#' @param by Grouping variable with at least two observed groups, or hierarchical `vars(...)` specification.
#' @param data Data frame. If `NULL`, active R4VN data is used.
#' @param posthoc Post-hoc method: `"none"`, `"dunn"`, or `"wilcoxon"`.
#' @param adjust Multiplicity adjustment accepted by [stats::p.adjust()].
#' @param effect Logical; add epsilon-squared and eta-squared(H).
#' @param digits,p_digits Decimal places for estimates and p-values.
#' @param show Logical; open the formatted result in the Viewer. Default `TRUE`.
#' @param console Logical; also print the traditional result in the Console. Default `FALSE`.
#'
#' @return Invisibly returns an object of class `r4vn_stat`.
#' @export
#'
#' @examples
#' d <- data.frame(
#'   score = c(10, 12, 11, 18, 17, 20, 25, 24, 27),
#'   group = factor(rep(c("A", "B", "C"), each = 3))
#' )
#' kwallis(score, by = group, data = d, posthoc = "dunn", effect = TRUE)
#'
#' # Extended usage examples
#' \donttest{
#' d <- data.frame(
#'   score = c(10, 12, 11, 18, 17, 20, 25, 24, 27),
#'   treatment = factor(rep(c("A", "B", "C"), each = 3))
#' )
#' kwallis(score, by = treatment, data = d)
#' usedf(d)
#' result <- kwallis(score, by = treatment, show = FALSE)
#' result$raw$test
#' }
kwallis <- function(x, by, data = NULL, digits = 3, p_digits = 3, show = TRUE, console = FALSE) {
  env <- parent.frame()
  data <- .r4vn_stat_data(data)
  x_expr <- substitute(x)
  by_expr <- substitute(by)
  xv <- .r4vn_eval_var(x_expr, data, env, "x")
  gv <- .r4vn_eval_var(by_expr, data, env, "by")
  if (!is.numeric(xv)) stop("`x` must be numeric.", call. = FALSE)
  ok <- is.finite(xv) & !is.na(gv)
  xv <- xv[ok]
  gf <- droplevels(factor(gv[ok]))
  if (nlevels(gf) < 2L) stop("`by` must have at least two observed groups.", call. = FALSE)
  groups <- split(xv, gf, drop = TRUE)
  desc <- do.call(rbind, Map(.r4vn_np_desc, groups, names(groups), MoreArgs = list(digits = digits)))
  fit <- stats::kruskal.test(xv, gf)
  test <- data.frame(
    Chi.square = .r4vn_num(unname(fit$statistic), digits),
    df = unname(fit$parameter),
    p = .r4vn_p(fit$p.value, p_digits),
    stringsAsFactors = FALSE,
    check.names = FALSE
  )
  result <- .r4vn_result(
    "Kruskal-Wallis rank-sum test",
    list("Descriptive statistics" = desc, "Test" = test),
    raw = list(x = xv, group = gf, test = fit),
    call = match.call()
  )
  .r4vn_show(result, show)
}

#' Friedman Test for Repeated or Matched Measurements
#'
#' Performs the Friedman rank-sum test for three or more repeated or matched
#' measurements. Data may be supplied in wide form as several numeric variables,
#' or in long form as one outcome, one occasion/treatment variable, and one
#' subject/block identifier.
#'
#' @param ... In wide form, two or more numeric repeated-measure variables. In
#'   long form, exactly one numeric outcome variable.
#' @param by Long-form occasion or treatment variable.
#' @param id Long-form subject or block identifier.
#' @param data Data frame. If `NULL`, active R4VN data is used.
#' @param digits,p_digits Decimal places for estimates and p-values.
#' @param show Logical; open the formatted result in the Viewer. Default `TRUE`.
#' @param console Logical; also print the traditional result in the Console. Default `FALSE`.
#'
#' @details
#' Wide form uses only rows complete across every repeated measurement. Long
#' form must contain one observation for every subject-by-occasion combination;
#' incomplete blocks are rejected by the underlying Friedman test.
#'
#' @return Invisibly returns an object of class `r4vn_stat`.
#' @export
#'
#' @examples
#' wide <- data.frame(
#'   baseline = c(22, 20, 19, 25, 23, 21),
#'   month1   = c(20, 18, 18, 22, 21, 20),
#'   month3   = c(18, 17, 16, 20, 19, 18)
#' )
#' friedman(baseline, month1, month3, data = wide)
#'
#' long <- data.frame(
#'   id = rep(1:6, each = 3),
#'   time = factor(rep(c("Baseline", "Month 1", "Month 3"), 6),
#'                 levels = c("Baseline", "Month 1", "Month 3")),
#'   score = as.vector(t(as.matrix(wide)))
#' )
#' friedman(score, by = time, id = id, data = long)
#'
#' # Extended usage examples
#' \donttest{
#' wide <- data.frame(
#'   baseline = c(22, 20, 19, 25, 23, 21),
#'   month1 = c(20, 18, 18, 22, 21, 20),
#'   month3 = c(18, 17, 16, 20, 19, 18)
#' )
#'
#' # Wide form: each row is a subject and each variable is an occasion
#' friedman(baseline, month1, month3, data = wide)
#'
#' # Long form: outcome, occasion, and subject identifier
#' long <- data.frame(
#'   id = rep(1:6, each = 3),
#'   time = factor(rep(c("Baseline", "Month 1", "Month 3"), 6),
#'                 levels = c("Baseline", "Month 1", "Month 3")),
#'   score = as.vector(t(as.matrix(wide)))
#' )
#' friedman(score, by = time, id = id, data = long)
#'
#' usedf(wide)
#' friedman(baseline, month1, month3)
#' }
friedman <- function(..., by = NULL, id = NULL, data = NULL,
                     digits = 3, p_digits = 3, show = TRUE, console = FALSE) {
  env <- parent.frame()
  data <- .r4vn_stat_data(data)
  exprs <- as.list(substitute(list(...)))[-1L]
  by_expr <- substitute(by)
  id_expr <- substitute(id)
  long_mode <- (!missing(by) && !identical(by_expr, quote(NULL))) ||
    (!missing(id) && !identical(id_expr, quote(NULL)))

  if (long_mode) {
    if (missing(by) || missing(id) || identical(by_expr, quote(NULL)) || identical(id_expr, quote(NULL))) {
      stop("Long-form syntax requires both `by` and `id`.", call. = FALSE)
    }
    if (length(exprs) != 1L) stop("Long-form syntax requires exactly one outcome in `...`.", call. = FALSE)
    y <- .r4vn_eval_var(exprs[[1L]], data, env, "outcome")
    g <- .r4vn_eval_var(by_expr, data, env, "by")
    b <- .r4vn_eval_var(id_expr, data, env, "id")
    if (!is.numeric(y)) stop("The outcome must be numeric.", call. = FALSE)
    ok <- is.finite(y) & !is.na(g) & !is.na(b)
    y <- y[ok]
    g <- droplevels(factor(g[ok]))
    b <- droplevels(factor(b[ok]))
    if (nlevels(g) < 2L) stop("`by` must have at least two observed occasions or treatments.", call. = FALSE)
    key <- interaction(b, g, drop = TRUE)
    if (anyDuplicated(key)) stop("Long-form data contain duplicate subject-by-occasion observations.", call. = FALSE)
    fit <- tryCatch(stats::friedman.test(y, g, b), error = function(e) {
      stop("Friedman test requires a complete unreplicated block design: ", conditionMessage(e), call. = FALSE)
    })
    groups <- split(y, g, drop = TRUE)
    desc <- do.call(rbind, Map(.r4vn_np_desc, groups, names(groups), MoreArgs = list(digits = digits)))
    raw <- list(y = y, treatment = g, block = b, test = fit)
  } else {
    if (length(exprs) < 2L) stop("Supply at least two repeated-measure variables.", call. = FALSE)
    values <- lapply(seq_along(exprs), function(i) {
      z <- .r4vn_eval_var(exprs[[i]], data, env, paste0("measurement ", i))
      if (!is.numeric(z)) stop("Every repeated measurement must be numeric.", call. = FALSE)
      z
    })
    m <- do.call(cbind, values)
    names0 <- vapply(exprs, .r4vn_name, character(1), fallback = "Measurement")
    colnames(m) <- make.unique(names0)
    m <- m[stats::complete.cases(m), , drop = FALSE]
    if (!nrow(m)) stop("No rows are complete across all repeated measurements.", call. = FALSE)
    fit <- stats::friedman.test(m)
    desc <- do.call(rbind, lapply(seq_len(ncol(m)), function(j) .r4vn_np_desc(m[, j], colnames(m)[j], digits)))
    raw <- list(matrix = m, test = fit)
  }

  test <- data.frame(
    Chi.square = .r4vn_num(unname(fit$statistic), digits),
    df = unname(fit$parameter),
    p = .r4vn_p(fit$p.value, p_digits),
    stringsAsFactors = FALSE,
    check.names = FALSE
  )
  result <- .r4vn_result(
    "Friedman rank-sum test",
    list("Descriptive statistics" = desc, "Test" = test),
    raw = raw,
    call = match.call()
  )
  .r4vn_show(result, show)
}

# Shapiro-Wilk Normality Test
#
# Performs the Shapiro-Wilk normality test overall or separately within groups.
# Each tested sample must contain between 3 and 5000 complete observations.
#
# @param x Numeric variable.
# @param by Optional grouping variable.
# @param data Data frame. If `NULL`, active R4VN data is used.
# @param digits,p_digits Decimal places for the W statistic and p-value.
# @param show Logical; open the formatted result in the Viewer. Default `TRUE`.
# @param console Logical; also print the traditional result in the Console. Default `FALSE`.
#
# @return Invisibly returns an object of class `r4vn_stat`.
# @export
#
# @examples
# d <- data.frame(
#   score = c(10, 11, 12, 13, 14, 16, 18, 20, 22, 25),
#   group = rep(c("A", "B"), each = 5)
# )
# swilk(score, data = d)
# swilk(score, by = group, data = d)
#
# # Extended usage examples
# \donttest{
# d <- data.frame(
#   score = c(10, 11, 12, 13, 14, 16, 18, 20, 22, 25),
#   group = factor(rep(c("A", "B"), each = 5))
# )
#
# # Overall and group-specific Shapiro-Wilk tests
# swilk(score, data = d)
# swilk(score, by = group, data = d)
#
# usedf(d)
# result <- swilk(score, by = group, show = FALSE)
# result$sections$Test
# }
.r4vn_swilk_legacy_source <- function(x, by = NULL, data = NULL, digits = 4, p_digits = 3, show = TRUE, console = FALSE) {
  env <- parent.frame()
  data <- .r4vn_stat_data(data)
  x_expr <- substitute(x)
  by_expr <- substitute(by)
  xv <- .r4vn_eval_var(x_expr, data, env, "x")
  if (!is.numeric(xv)) stop("`x` must be numeric.", call. = FALSE)
  has_by <- !missing(by) && !identical(by_expr, quote(NULL))

  run_one <- function(z, label) {
    z <- z[is.finite(z)]
    if (length(z) < 3L || length(z) > 5000L) {
      stop("Each Shapiro-Wilk sample must contain 3 to 5000 complete observations; group `", label, "` has ", length(z), ".", call. = FALSE)
    }
    fit <- stats::shapiro.test(z)
    data.frame(Group = label, n = length(z), W = .r4vn_num(unname(fit$statistic), digits),
               p = .r4vn_p(fit$p.value, p_digits), stringsAsFactors = FALSE,
               check.names = FALSE)
  }

  if (has_by) {
    gv <- .r4vn_eval_var(by_expr, data, env, "by")
    ok <- is.finite(xv) & !is.na(gv)
    groups <- split(xv[ok], droplevels(factor(gv[ok])), drop = TRUE)
    if (!length(groups)) stop("No complete grouped observations are available.", call. = FALSE)
    tab <- do.call(rbind, Map(run_one, groups, names(groups)))
    raw <- list(groups = groups)
  } else {
    z <- xv[is.finite(xv)]
    tab <- run_one(z, .r4vn_name(x_expr, "x"))
    raw <- list(x = z, test = stats::shapiro.test(z))
  }
  result <- .r4vn_result(
    "Shapiro-Wilk normality test",
    list("Test" = tab),
    notes = "A nonsignificant p-value does not prove normality; inspect the distribution and study design as well.",
    raw = raw,
    call = match.call()
  )
  .r4vn_show(result, show)
}


# ============================================================================
# Viewer renderers for non-parametric commands
# ============================================================================
.r4vn_nonparametric_viewer <- function(x, subtitle) {
  blocks <- paste0(vapply(names(x$sections), function(nm) {
    .r4vn_view_section(nm, x$sections[[nm]])
  }, character(1)), collapse = "")
  .r4vn_view_document(x$title, blocks, notes = x$notes, subtitle = subtitle,
                      prefix = "r4vn-nonparametric-")
}

.r4vn_viewer_ranksum <- function(x) .r4vn_nonparametric_viewer(x, "Wilcoxon rank-sum / Mann-Whitney analysis")
.r4vn_viewer_signrank <- function(x) .r4vn_nonparametric_viewer(x, "Wilcoxon signed-rank analysis")
.r4vn_viewer_kwallis <- function(x) .r4vn_nonparametric_viewer(x, "Kruskal-Wallis analysis")
.r4vn_viewer_friedman <- function(x) .r4vn_nonparametric_viewer(x, "Friedman repeated-measures analysis")
.r4vn_viewer_swilk <- function(x) .r4vn_nonparametric_viewer(x, "Shapiro-Wilk normality analysis")