R/statistics.R

In R4VN: Health Data Analysis and Publication-Ready Reporting

# R4VN data-based statistical commands
#
# All commands accept `data`. When `data = NULL`, the active data set selected
# by usedf() or opendata(..., active = TRUE) is used.
#
# Display convention used by this file:
#   show = TRUE      open the RStudio Viewer / browser result (default)
#   console = FALSE  do not duplicate the result in Console (default)
# Set console = TRUE when the traditional Console output is also wanted.

.r4vn_statistics_viewer <- function(x, subtitle = "Analysis using individual-level data") {
  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-statistics-")
}

.r4vn_statistics_display <- function(x, call, show = TRUE, console = FALSE) {
  x$call <- call
  .r4vn_show(x, show = show, console = console)
}

# Dedicated Viewer entry points. They intentionally remain separate even when
# their current layouts are similar, so each command can evolve independently.
.r4vn_viewer_ttest <- function(x) .r4vn_statistics_viewer(x, "t test using individual-level data")
.r4vn_viewer_sdtest <- function(x) .r4vn_statistics_viewer(x, "Variance comparison using individual-level data")
.r4vn_viewer_anova <- function(x) .r4vn_statistics_viewer(x, "One-way analysis of variance")
.r4vn_viewer_prop <- function(x) .r4vn_statistics_viewer(x, "Proportion estimate and confidence interval")
.r4vn_viewer_bitest <- function(x) .r4vn_statistics_viewer(x, "Exact binomial test")
.r4vn_viewer_prtest <- function(x) .r4vn_statistics_viewer(x, "Proportion z test")
.r4vn_viewer_ci <- function(x) .r4vn_statistics_viewer(x, "Confidence interval")
.r4vn_viewer_esize <- function(x) .r4vn_statistics_viewer(x, "Standardized effect size")
.r4vn_viewer_ztest <- function(x) .r4vn_statistics_viewer(x, "z test using known standard deviation")

# ============================================================================
# Function source: ttest.R
# ============================================================================
#' @rdname ttest
#' @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,
                  digits = 3, p_digits = 3, show = TRUE, console = FALSE) {
  call <- match.call(); env <- parent.frame(); data <- .r4vn_stat_data(data); alternative <- match.arg(alternative)
  xv <- .r4vn_eval_var(substitute(x), data, env, "x"); if (!is.numeric(xv)) stop("`x` must be numeric.", call. = FALSE)
  if (!missing(by) && !identical(substitute(by), quote(NULL))) {
    gv <- .r4vn_eval_var(substitute(by), data, env, "by"); ok <- !is.na(xv) & !is.na(gv); gv <- droplevels(factor(gv[ok])); xv <- xv[ok]
    if (nlevels(gv) != 2L) stop("`by` must have exactly two observed groups.", call. = FALSE)
    s <- split(xv, gv)
    out <- ttesti(length(s[[1]]), mean(s[[1]]), stats::sd(s[[1]]), length(s[[2]]), mean(s[[2]]), stats::sd(s[[2]]), equal = equal, alternative = alternative, level = level, digits = digits, p_digits = p_digits, show = FALSE, console = FALSE)
    out$sections[["Descriptive statistics"]]$Group <- names(s)
    return(.r4vn_statistics_display(out, call, show, console))
  }
  if (!missing(y) && !identical(substitute(y), quote(NULL))) {
    yv <- .r4vn_eval_var(substitute(y), data, env, "y"); if (!is.numeric(yv)) stop("`y` must be numeric.", call. = FALSE)
    if (paired) {
      ok <- stats::complete.cases(xv, yv); xv <- xv[ok]; yv <- yv[ok]; rr <- stats::cor(xv, yv)
      out <- ttesti(length(xv), mean(xv), stats::sd(xv), length(yv), mean(yv), stats::sd(yv), paired = TRUE, r = rr, alternative = alternative, level = level, digits = digits, p_digits = p_digits, show = FALSE, console = FALSE)
    } else {
      xv <- xv[!is.na(xv)]; yv <- yv[!is.na(yv)]
      out <- ttesti(length(xv), mean(xv), stats::sd(xv), length(yv), mean(yv), stats::sd(yv), equal = equal, alternative = alternative, level = level, digits = digits, p_digits = p_digits, show = FALSE, console = FALSE)
    }
    return(.r4vn_statistics_display(out, call, show, console))
  }
  xv <- xv[!is.na(xv)]
  out <- ttesti(length(xv), mean(xv), stats::sd(xv), mu = mu, alternative = alternative, level = level, digits = digits, p_digits = p_digits, show = FALSE, console = FALSE)
  .r4vn_statistics_display(out, call, show, console)
}

# ============================================================================
# Function source: sdtest.R
# ============================================================================
#' @rdname sdtest
#' @export
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(); data <- .r4vn_stat_data(data); alternative <- match.arg(alternative)
  xv <- .r4vn_eval_var(substitute(x), data, env, "x"); if (!is.numeric(xv)) stop("`x` must be numeric.", call. = FALSE)
  if (!missing(by) && !identical(substitute(by), quote(NULL))) {
    g <- .r4vn_eval_var(substitute(by), data, env, "by"); ok <- stats::complete.cases(xv, g); s <- split(xv[ok], droplevels(factor(g[ok]))); if (length(s) != 2L) stop("`by` must have two groups.", call. = FALSE)
    out <- sdtesti(length(s[[1]]), stats::sd(s[[1]]), length(s[[2]]), stats::sd(s[[2]]), alternative = alternative, level = level, digits = digits, p_digits = p_digits, show = FALSE, console = FALSE)
    out$sections[["Descriptive statistics"]]$Group <- names(s)
    return(.r4vn_statistics_display(out, call, show, console))
  }
  if (!missing(y) && !identical(substitute(y), quote(NULL))) {
    yv <- .r4vn_eval_var(substitute(y), data, env, "y"); xv <- xv[!is.na(xv)]; yv <- yv[!is.na(yv)]
    out <- sdtesti(length(xv), stats::sd(xv), length(yv), stats::sd(yv), alternative = alternative, level = level, digits = digits, p_digits = p_digits, show = FALSE, console = FALSE)
  } else {
    xv <- xv[!is.na(xv)]
    out <- sdtesti(length(xv), stats::sd(xv), sd0 = sd0, alternative = alternative, level = level, digits = digits, p_digits = p_digits, show = FALSE, console = FALSE)
  }
  .r4vn_statistics_display(out, call, show, console)
}

# ============================================================================
# Function source: anova.R
# ============================================================================
#' @rdname anova_r4vn
#' @export
anova <- function(object, ..., by = NULL, data = NULL, bartlett = TRUE, level = 0.95,
                  digits = 3, p_digits = 3, show = TRUE, console = FALSE) {
  call <- match.call(); env <- parent.frame(); object_expr <- substitute(object)
  obj <- if (is.null(data) && (missing(by) || identical(substitute(by), quote(NULL)))) tryCatch(eval(object_expr, env), error = function(e) NULL) else NULL
  if (!is.null(obj) && (inherits(obj, c("lm", "glm", "aov", "mlm", "lme", "gls")) || length(list(...)) > 0L)) {
    tab <- stats::anova(obj, ...)
    shown <- .r4vn_result("Analysis of variance", list("ANOVA table" = as.data.frame(tab, check.names = FALSE)),
                          raw = list(anova = tab), call = call)
    .r4vn_show(shown, show = show, console = console)
    return(invisible(tab))
  }
  if (missing(by) || identical(substitute(by), quote(NULL))) stop("Supply a grouping variable in `by` for one-way ANOVA.", call. = FALSE)
  out <- .r4vn_anova_data(object_expr, substitute(by), data, env, bartlett, level, digits, p_digits, show = FALSE)
  .r4vn_statistics_display(out, call, show, console)
}

# ============================================================================
# Function source: prop.R
# ============================================================================
#' @rdname proportion_tests
#' @export
prop <- function(x, data = NULL, event = NULL, p0 = NULL, method = c("all", "exact", "wilson", "wald"),
                 alternative = c("two.sided", "less", "greater"), level = 0.95,
                 digits = 3, p_digits = 3, show = TRUE, console = FALSE) {
  call <- match.call(); env <- parent.frame(); data <- .r4vn_stat_data(data); v <- .r4vn_eval_var(substitute(x), data, env, "x"); v <- v[!is.na(v)]; lev <- unique(as.character(v)); ev <- as.character(event %||% tail(lev, 1L)); if (!ev %in% lev) stop("`event` was not found.", call. = FALSE)
  out <- propi(sum(as.character(v) == ev), length(v), p0, method, alternative, level, digits, p_digits, show = FALSE, console = FALSE)
  .r4vn_statistics_display(out, call, show, console)
}

# ============================================================================
# Function source: bitest.R
# ============================================================================
#' @rdname proportion_tests
#' @export
bitest <- function(x, p = 0.5, data = NULL, event = 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(); data <- .r4vn_stat_data(data)
  v <- .r4vn_eval_var(substitute(x), data, env, "x"); v <- v[!is.na(v)]
  lev <- unique(as.character(v)); ev <- as.character(event %||% tail(lev, 1L)); if (!ev %in% lev) stop("`event` was not found.", call. = FALSE)
  out <- bitesti(length(v), sum(as.character(v) == ev), p = p, alternative = alternative, level = level, digits = digits, p_digits = p_digits, show = FALSE, console = FALSE)
  .r4vn_statistics_display(out, call, show, console)
}

# ============================================================================
# Function source: prtest.R
# ============================================================================
#' @rdname proportion_tests
#' @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(); data <- .r4vn_stat_data(data); v <- .r4vn_eval_var(substitute(x), data, env, "x"); lev <- unique(as.character(v[!is.na(v)])); ev <- as.character(event %||% tail(lev, 1L))
  if (missing(by) || identical(substitute(by), quote(NULL))) {
    ok <- !is.na(v); out <- prtesti(sum(as.character(v[ok]) == ev), sum(ok), p0 = p0, alternative = alternative, level = level, digits = digits, p_digits = p_digits, show = FALSE, console = FALSE)
  } else {
    g <- .r4vn_eval_var(substitute(by), data, env, "by"); ok <- stats::complete.cases(v, g); g <- droplevels(factor(g[ok])); v <- v[ok]; if (nlevels(g) != 2L) stop("`by` must have two groups.", call. = FALSE); s <- split(as.character(v) == ev, g)
    out <- prtesti(sum(s[[1]]), length(s[[1]]), sum(s[[2]]), length(s[[2]]), alternative = alternative, level = level, digits = digits, p_digits = p_digits, show = FALSE, console = FALSE)
    out$sections[["Group proportions"]]$Group <- names(s)
  }
  .r4vn_statistics_display(out, call, show, console)
}

# ============================================================================
# Function source: ci.R
# ============================================================================
#' @rdname ci_r4vn
#' @export
ci <- function(x, data = NULL, type = c("auto", "mean", "proportion", "variance"), event = NULL,
               method = c("exact", "wilson", "wald"), level = 0.95, digits = 3,
               show = TRUE, console = FALSE) {
  call <- match.call(); env <- parent.frame(); data <- .r4vn_stat_data(data); v <- .r4vn_eval_var(substitute(x), data, env, "x"); v <- v[!is.na(v)]; type <- match.arg(type)
  if (type == "auto") type <- if (is.numeric(v) && length(unique(v)) > 2L) "mean" else "proportion"
  if (type == "mean") out <- cii(length(v), mean(v), stats::sd(v), type = "mean", level = level, digits = digits, show = FALSE, console = FALSE)
  else if (type == "variance") out <- cii(length(v), variance = stats::var(v), type = "variance", level = level, digits = digits, show = FALSE, console = FALSE)
  else {
    lev <- unique(as.character(v)); ev <- as.character(event %||% tail(lev, 1L))
    out <- cii(length(v), events = sum(as.character(v) == ev), type = "proportion", method = method, level = level, digits = digits, show = FALSE, console = FALSE)
  }
  .r4vn_statistics_display(out, call, show, console)
}

# ============================================================================
# Function source: esize.R
# ============================================================================
#' @rdname esize
#' @export
esize <- function(x, by, data = NULL, level = 0.95, digits = 3,
                  show = TRUE, console = FALSE) {
  call <- match.call(); env <- parent.frame(); data <- .r4vn_stat_data(data); xv <- .r4vn_eval_var(substitute(x), data, env, "x"); g <- .r4vn_eval_var(substitute(by), data, env, "by"); ok <- stats::complete.cases(xv, g); s <- split(xv[ok], droplevels(factor(g[ok]))); if (length(s) != 2L) stop("`by` must have two groups.", call. = FALSE)
  out <- esizei(length(s[[1]]), mean(s[[1]]), stats::sd(s[[1]]), length(s[[2]]), mean(s[[2]]), stats::sd(s[[2]]), level = level, digits = digits, show = FALSE, console = FALSE)
  .r4vn_statistics_display(out, call, show, console)
}

# ============================================================================
# Function source: ztest.R
# ============================================================================
#' @rdname ztest
#' @export
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(); data <- .r4vn_stat_data(data); xv <- .r4vn_eval_var(substitute(x), data, env, "x")
  if (!missing(by) && !identical(substitute(by), quote(NULL))) {
    g <- .r4vn_eval_var(substitute(by), data, env, "by"); ok <- stats::complete.cases(xv, g); s <- split(xv[ok], droplevels(factor(g[ok]))); if (length(s) != 2L || is.null(sigma2)) stop("Two groups and `sigma2` are required.", call. = FALSE)
    out <- ztesti(length(s[[1]]), mean(s[[1]]), sigma1, length(s[[2]]), mean(s[[2]]), sigma2, alternative = alternative, level = level, digits = digits, p_digits = p_digits, show = FALSE, console = FALSE)
  } else if (!missing(y) && !identical(substitute(y), quote(NULL))) {
    yv <- .r4vn_eval_var(substitute(y), data, env, "y"); xv <- xv[!is.na(xv)]; yv <- yv[!is.na(yv)]
    out <- ztesti(length(xv), mean(xv), sigma1, length(yv), mean(yv), sigma2, alternative = alternative, level = level, digits = digits, p_digits = p_digits, show = FALSE, console = FALSE)
  } else {
    xv <- xv[!is.na(xv)]
    out <- ztesti(length(xv), mean(xv), sigma1, mu = mu, alternative = alternative, level = level, digits = digits, p_digits = p_digits, show = FALSE, console = FALSE)
  }
  .r4vn_statistics_display(out, call, show, console)
}