R/normality.R

In R4VN: Health Data Analysis and Publication-Ready Reporting

# ============================================================================
# R4VN normality/distribution tests
# ============================================================================

.r4vn_skew_kurt <- function(x) {
  x <- x[is.finite(x)]
  n <- length(x)
  if (n < 2L) return(c(skewness = NA_real_, kurtosis = NA_real_))
  m <- mean(x)
  m2 <- mean((x - m)^2)
  if (!is.finite(m2) || m2 <= 0) return(c(skewness = 0, kurtosis = NA_real_))
  m3 <- mean((x - m)^3)
  m4 <- mean((x - m)^4)
  c(skewness = m3 / m2^(3/2), kurtosis = m4 / m2^2)
}

.r4vn_jarque_bera <- function(x) {
  x <- x[is.finite(x)]
  n <- length(x)
  sk <- .r4vn_skew_kurt(x)
  stat <- n / 6 * (sk[["skewness"]]^2 + (sk[["kurtosis"]] - 3)^2 / 4)
  list(statistic = stat, p.value = stats::pchisq(stat, df = 2, lower.tail = FALSE), df = 2)
}

.r4vn_normality_one <- function(x, method) {
  x <- x[is.finite(x)]
  n <- length(x)
  if (n < 3L) {
    return(list(statistic = NA_real_, p.value = NA_real_, note = "At least 3 observations are required."))
  }
  if (length(unique(x)) < 2L) {
    return(list(statistic = NA_real_, p.value = NA_real_, note = "The sample has no variability."))
  }

  if (method == "shapiro") {
    if (n > 5000L) return(list(statistic = NA_real_, p.value = NA_real_, note = "Shapiro-Wilk is limited to 5000 observations."))
    fit <- stats::shapiro.test(x)
    return(list(statistic = unname(fit$statistic), p.value = fit$p.value, note = NULL))
  }

  if (method == "ks") {
    sx <- stats::sd(x)
    if (!is.finite(sx) || sx <= 0) return(list(statistic = NA_real_, p.value = NA_real_, note = "Standard deviation is zero."))
    fit <- suppressWarnings(stats::ks.test(x, "pnorm", mean(x), sx, exact = FALSE))
    return(list(
      statistic = unname(fit$statistic), p.value = fit$p.value,
      note = "Mean and SD are estimated from the sample; the ordinary KS p-value is approximate. Lilliefors is preferred when available."
    ))
  }

  if (method == "jarque.bera") {
    fit <- .r4vn_jarque_bera(x)
    return(list(statistic = fit$statistic, p.value = fit$p.value, note = NULL))
  }

  nortest_map <- c(
    lilliefors = "lillie.test",
    anderson = "ad.test",
    cramer.von.mises = "cvm.test",
    shapiro.francia = "sf.test",
    pearson = "pearson.test"
  )
  if (method %in% names(nortest_map)) {
    if (!requireNamespace("nortest", quietly = TRUE)) {
      return(list(
        statistic = NA_real_, p.value = NA_real_,
        note = "Install the optional `nortest` package to run this method."
      ))
    }
    fun <- get(nortest_map[[method]], envir = asNamespace("nortest"))
    fit <- tryCatch(fun(x), error = function(e) e)
    if (inherits(fit, "error")) {
      return(list(statistic = NA_real_, p.value = NA_real_, note = conditionMessage(fit)))
    }
    return(list(statistic = unname(fit$statistic[1L]), p.value = fit$p.value, note = NULL))
  }

  stop("Unknown normality-test method.", call. = FALSE)
}

.r4vn_normality_method_label <- function(x) {
  switch(x,
    shapiro = "Shapiro-Wilk",
    ks = "Kolmogorov-Smirnov (fitted normal)",
    lilliefors = "Lilliefors",
    anderson = "Anderson-Darling",
    cramer.von.mises = "Cramer-von Mises",
    shapiro.francia = "Shapiro-Francia",
    pearson = "Pearson chi-square",
    jarque.bera = "Jarque-Bera",
    x
  )
}

#' Normality tests for one or more variables
#'
#' Runs a consistent set of normality/distribution tests overall or within
#' groups. `vars()` may select many numeric variables in one call, and the R4VN
#' hierarchical convention is supported: `by = vars(region, sex)` means test
#' separately within each region and then within sex inside each region.
#'
#' @usage normtest(x = NULL, vars = NULL, by = NULL, data = NULL, method = c("all", "shapiro", "ks", "lilliefors", "anderson", "cramer.von.mises", "shapiro.francia", "pearson", "jarque.bera"), digits = 4, p_digits = 3, show = TRUE, console = FALSE)
#'
#' @param x Optional numeric variable. `x = vars(age, bmi)` is also accepted.
#' @param vars Optional `vars(...)` selector for testing several numeric variables.
#' @param by Optional grouping specification. A single variable gives group-specific tests. With `by = vars(region, sex, outcome)`, `region` and `sex` are nested strata and `outcome` is the innermost grouping variable.
#' @param data Data frame. When omitted, the active R4VN data frame is used.
#' @param method Test method. `"all"` runs Shapiro-Wilk, fitted-normal Kolmogorov-Smirnov, Lilliefors, Anderson-Darling, Cramer-von Mises, Shapiro-Francia, Pearson chi-square, and Jarque-Bera. Methods provided by the optional `nortest` package are reported as unavailable rather than causing `method = "all"` to fail when that package is absent.
#' @param digits,p_digits Decimal places for test statistics and p-values.
#' @param show Logical; open the formatted result in the Viewer.
#' @param console Logical; also print the result in the Console.
#'
#' @details
#' The ordinary one-sample Kolmogorov-Smirnov test assumes a fully specified
#' reference distribution. `normtest(method = "ks")` estimates the mean and
#' standard deviation from the same sample for convenience, so its p-value is
#' only approximate; the Lilliefors method is generally preferable for this
#' use. Normality tests should be interpreted together with histograms, Q-Q
#' plots, sample size, and the intended analysis rather than as a mechanical
#' pass/fail rule.
#'
#' @return Invisibly returns an object of class `r4vn_stat`. The `Test` section contains one row per variable, stratum, group, and method.
#' @export
#'
#' @examples
#' d <- data.frame(
#'   age = c(31, 42, 38, 50, 46, 35, 55, 61, 44, 39),
#'   bmi = c(21, 24, 26, 23, 28, 20, 25, 29, 22, 27),
#'   sex = factor(rep(c("Female", "Male"), 5)),
#'   region = factor(rep(c("North", "South"), each = 5))
#' )
#' normtest(age, data = d, method = "shapiro", show = FALSE)
#' normtest(vars = vars(age, bmi), data = d, method = c("shapiro", "jarque.bera"), show = FALSE)
#' normtest(vars = vars(age, bmi), by = sex, data = d, method = "shapiro", show = FALSE)
#' normtest(vars = vars(age, bmi), by = vars(region, sex), data = d, method = "shapiro", show = FALSE)
#'
#' @seealso [swilk()], [varform()], [ghist()]
normtest <- function(x = NULL, vars = NULL, by = NULL, data = NULL,
                     method = c("all", "shapiro", "ks", "lilliefors", "anderson",
                                "cramer.von.mises", "shapiro.francia", "pearson", "jarque.bera"),
                     digits = 4, p_digits = 3, show = TRUE, console = FALSE) {
  call <- match.call()
  env <- parent.frame()
  data <- .r4vn_stat_data(data)
  x_expr <- substitute(x)
  vars_expr <- substitute(vars)
  by_expr <- substitute(by)

  if (!missing(vars) && !.r4vn_expr_is_null(vars_expr)) {
    var_names <- .r4vn_vars_spec_names(vars_expr, data, env, arg = "vars", numeric_only = TRUE)
  } else if (!missing(x) && !.r4vn_expr_is_null(x_expr)) {
    if (.r4vn_is_vars_call(x_expr)) {
      var_names <- .r4vn_vars_spec_names(x_expr, data, env, arg = "x", numeric_only = TRUE)
    } else if (is.symbol(x_expr) && as.character(x_expr) %in% names(data)) {
      var_names <- as.character(x_expr)
    } else {
      x_value <- tryCatch(eval(x_expr, envir = env), error = function(e) NULL)
      if (is.character(x_value) && length(x_value) && all(x_value %in% names(data))) {
        var_names <- as.character(x_value)
      } else {
      # Preserve support for an arbitrary numeric expression by placing it in a
      # temporary analysis column without mutating the user's data.
      value <- .r4vn_eval_var(x_expr, data, env, "x")
      if (!is.numeric(value)) stop("`x` must be numeric.", call. = FALSE)
      temp_name <- ".r4vn_normtest_x"
      while (temp_name %in% names(data)) temp_name <- paste0(temp_name, "_")
        data[[temp_name]] <- value
        var_names <- temp_name
      }
    }
  } else {
    stop("Supply `x` or `vars = vars(...)`.", call. = FALSE)
  }

  methods <- as.character(method)
  allowed <- c("all", "shapiro", "ks", "lilliefors", "anderson",
               "cramer.von.mises", "shapiro.francia", "pearson", "jarque.bera")
  bad <- setdiff(methods, allowed)
  if (length(bad)) stop("Unknown `method`: ", paste(bad, collapse = ", "), ".", call. = FALSE)
  if ("all" %in% methods) methods <- setdiff(allowed, "all")
  methods <- unique(methods)

  by_spec <- .r4vn_by_spec(by_expr, data, env, allow_null = TRUE)
  strata <- .r4vn_strata_indices(data, by_spec$strata)
  if (!length(strata)) stop("No observations remain after stratification.", call. = FALSE)

  rows <- list()
  raw <- list()
  notes <- character()
  r <- 0L

  for (vn in var_names) {
    if (!is.numeric(data[[vn]])) stop("Variable `", vn, "` must be numeric.", call. = FALSE)
    for (si in seq_along(strata)) {
      idx <- strata[[si]]
      stratum_label <- .r4vn_stratum_label(data, by_spec$strata, idx)
      group_indices <- if (is.null(by_spec$by)) {
        list(Overall = idx)
      } else {
        g <- data[[by_spec$by]][idx]
        good <- !is.na(g)
        split(idx[good], droplevels(factor(g[good])), drop = TRUE)
      }
      if (!length(group_indices)) next

      for (gn in names(group_indices)) {
        z <- data[[vn]][group_indices[[gn]]]
        z <- z[is.finite(z)]
        for (m in methods) {
          fit <- .r4vn_normality_one(z, m)
          r <- r + 1L
          rows[[r]] <- data.frame(
            Variable = if (identical(vn, ".r4vn_normtest_x")) .r4vn_name(x_expr, "x") else vn,
            Stratum = if (length(by_spec$strata)) stratum_label else "",
            Group = if (is.null(by_spec$by)) "" else gn,
            Method = .r4vn_normality_method_label(m),
            n = length(z),
            Statistic = .r4vn_num(fit$statistic, digits),
            p = .r4vn_p(fit$p.value, p_digits),
            stringsAsFactors = FALSE,
            check.names = FALSE
          )
          raw[[paste(vn, stratum_label, gn, m, sep = "|")]] <- fit
          if (!is.null(fit$note) && nzchar(fit$note)) notes <- c(notes, paste0(.r4vn_normality_method_label(m), ": ", fit$note))
        }
      }
    }
  }

  tab <- if (length(rows)) do.call(rbind, rows) else data.frame()
  if (!length(by_spec$strata) && "Stratum" %in% names(tab)) tab$Stratum <- NULL
  if (is.null(by_spec$by) && "Group" %in% names(tab)) tab$Group <- NULL

  result <- .r4vn_result(
    "Normality tests",
    list("Test" = tab),
    notes = unique(c(
      notes,
      "A nonsignificant p-value does not prove normality; inspect the distribution and the planned analysis as well."
    )),
    raw = list(tests = raw, variables = var_names, by = by_spec),
    call = call
  )
  .r4vn_show(result, show = show, console = console)
}

#' Shapiro-Wilk normality test
#'
#' Convenience wrapper for `normtest(method = "shapiro")`. It retains the
#' familiar `swilk(x)` command while adding `vars()` and hierarchical `by`.
#'
#' @usage swilk(x = NULL, vars = NULL, by = NULL, data = NULL, digits = 4, p_digits = 3, show = TRUE, console = FALSE)
#'
#' @param x Optional numeric variable.
#' @param vars Optional `vars(...)` selector for several numeric variables.
#' @param by Optional group or hierarchical `by = vars(...)` specification.
#' @param data Data frame; active data are used when omitted.
#' @param digits,p_digits Decimal places.
#' @param show Logical; show the formatted result.
#' @param console Logical; also print to Console.
#' @return Invisibly returns an `r4vn_stat` object.
#' @export
#' @examples
#' d <- data.frame(x = rnorm(20), y = rnorm(20), g = rep(c("A", "B"), 10))
#' swilk(x, data = d, show = FALSE)
#' swilk(vars = vars(x, y), by = g, data = d, show = FALSE)
swilk <- function(x = NULL, vars = NULL, by = NULL, data = NULL,
                  digits = 4, p_digits = 3, show = TRUE, console = FALSE) {
  env <- parent.frame()
  cl <- match.call()
  requested_show <- show
  requested_console <- console
  cl[[1L]] <- quote(normtest)
  cl$method <- "shapiro"
  cl$show <- FALSE
  cl$console <- FALSE
  out <- eval(cl, envir = env)
  out$title <- "Shapiro-Wilk normality test"
  out$call <- match.call()
  .r4vn_show(out, show = requested_show, console = requested_console)
}