R/immediate.R

In R4VN: Health Data Analysis and Publication-Ready Reporting

# R4VN immediate statistical commands
#
# These functions calculate statistics from typed summary data rather than
# individual-level observations. Each public command remains an independent R
# function even though the functions are stored together in this source file.

# ============================================================================
# Function source: tabi.R
# ============================================================================
#' @rdname tabi
#' @export
tabi <- function(..., row.names = NULL, col.names = NULL, percent = c("none", "row", "col", "total"),
                 row = FALSE, col = FALSE, cell = FALSE, total = FALSE,
                 exp = FALSE, chi = TRUE, fisher = FALSE, lr = FALSE, residual = FALSE,
                 adjresidual = FALSE, correct = FALSE, digits = 1, p_digits = 3,
                 workspace = 2e5, show = TRUE, console = FALSE) {
  m <- .r4vn_matrix_input(..., row.names = row.names, col.names = col.names)
  percent <- .r4vn_percent_option(percent, row, col, cell, total)
  z <- .r4vn_tab_calc(m, percent, exp, chi, fisher, lr, residual, adjresidual, correct, digits, p_digits, workspace)
  .r4vn_show(.r4vn_result("Immediate contingency table", z$sections, raw = z$raw, call = match.call()), show)
}

# ============================================================================
# Function source: ttesti.R
# ============================================================================
#' @rdname ttest
#' @export
ttesti <- function(n1, mean1, sd1, n2 = NULL, mean2 = NULL, sd2 = NULL, mu = 0,
                   equal = FALSE, paired = FALSE, r = NULL,
                   alternative = c("two.sided", "less", "greater"), level = 0.95,
                   digits = 3, p_digits = 3, show = TRUE, console = FALSE) {
  alternative <- match.arg(alternative)
  vals <- c(n1, mean1, sd1); if (any(!is.finite(vals)) || n1 < 2 || sd1 < 0) stop("Invalid first-group summary statistics.", call. = FALSE)
  one <- is.null(n2) && is.null(mean2) && is.null(sd2)
  g1 <- .r4vn_t_group(n1, mean1, sd1, level)
  if (one) {
    se <- sd1 / sqrt(n1); df <- n1 - 1; t <- (mean1 - mu) / se
    p <- switch(alternative, two.sided = 2 * stats::pt(abs(t), df, lower.tail = FALSE), less = stats::pt(t, df), greater = stats::pt(t, df, lower.tail = FALSE))
    desc <- data.frame(Group = "x", n = n1, Mean = .r4vn_num(mean1, digits), SE = .r4vn_num(g1["se"], digits), SD = .r4vn_num(sd1, digits), CI = paste0(.r4vn_num(g1["lower"], digits), " to ", .r4vn_num(g1["upper"], digits)), stringsAsFactors = FALSE)
    test <- data.frame(Contrast = paste0("Mean - ", mu), Difference = .r4vn_num(mean1 - mu, digits), t = .r4vn_num(t, digits), df = .r4vn_num(df, 2), p = .r4vn_p(p, p_digits), stringsAsFactors = FALSE)
    raw <- list(n = n1, mean = mean1, sd = sd1, mu = mu, t = t, df = df, p.value = p)
    return(.r4vn_show(.r4vn_result("One-sample t test", list("Descriptive statistics" = desc, "Test" = test), raw = raw, call = match.call()), show))
  }
  if (any(vapply(list(n2, mean2, sd2), is.null, logical(1)))) stop("For two groups, supply `n2`, `mean2`, and `sd2`.", call. = FALSE)
  if (n2 < 2 || sd2 < 0 || any(!is.finite(c(n2, mean2, sd2)))) stop("Invalid second-group summary statistics.", call. = FALSE)
  g2 <- .r4vn_t_group(n2, mean2, sd2, level); diff <- mean1 - mean2
  if (paired) {
    if (n1 != n2) stop("Paired summaries require equal sample sizes.", call. = FALSE)
    if (is.null(r) || !is.finite(r) || abs(r) > 1) stop("Supply correlation `r` between -1 and 1 for a paired test.", call. = FALSE)
    sd_diff <- sqrt(sd1^2 + sd2^2 - 2 * r * sd1 * sd2); se <- sd_diff / sqrt(n1); df <- n1 - 1
  } else if (equal) {
    df <- n1 + n2 - 2; sp2 <- ((n1 - 1) * sd1^2 + (n2 - 1) * sd2^2) / df; se <- sqrt(sp2 * (1 / n1 + 1 / n2))
  } else {
    v1 <- sd1^2 / n1; v2 <- sd2^2 / n2; se <- sqrt(v1 + v2); df <- (v1 + v2)^2 / (v1^2 / (n1 - 1) + v2^2 / (n2 - 1))
  }
  t <- diff / se; crit <- stats::qt(1 - (1 - level) / 2, df); ci <- diff + c(-1, 1) * crit * se
  p <- switch(alternative, two.sided = 2 * stats::pt(abs(t), df, lower.tail = FALSE), less = stats::pt(t, df), greater = stats::pt(t, df, lower.tail = FALSE))
  desc <- rbind(
    data.frame(Group = "Group 1", n = n1, Mean = .r4vn_num(mean1, digits), SE = .r4vn_num(g1["se"], digits), SD = .r4vn_num(sd1, digits), CI = paste0(.r4vn_num(g1["lower"], digits), " to ", .r4vn_num(g1["upper"], digits)), stringsAsFactors = FALSE),
    data.frame(Group = "Group 2", n = n2, Mean = .r4vn_num(mean2, digits), SE = .r4vn_num(g2["se"], digits), SD = .r4vn_num(sd2, digits), CI = paste0(.r4vn_num(g2["lower"], digits), " to ", .r4vn_num(g2["upper"], digits)), stringsAsFactors = FALSE)
  )
  test <- data.frame(Contrast = "Mean 1 - Mean 2", Difference = .r4vn_num(diff, digits), CI = paste0(.r4vn_num(ci[1], digits), " to ", .r4vn_num(ci[2], digits)), t = .r4vn_num(t, digits), df = .r4vn_num(df, 2), p = .r4vn_p(p, p_digits), stringsAsFactors = FALSE)
  raw <- list(group1 = g1, group2 = g2, difference = diff, se = se, t = t, df = df, p.value = p, conf.int = ci)
  .r4vn_show(.r4vn_result(if (paired) "Paired t test from summaries" else if (equal) "Two-sample t test with equal variances" else "Welch two-sample t test", list("Descriptive statistics" = desc, "Test" = test), raw = raw, call = match.call()), show)
}

# ============================================================================
# Function source: sdtesti.R
# ============================================================================
#' @rdname sdtest
#' @export
sdtesti <- function(n1, sd1, n2 = NULL, sd2 = NULL, sd0 = NULL,
                    alternative = c("two.sided", "less", "greater"), level = 0.95,
                    digits = 3, p_digits = 3, show = TRUE, console = FALSE) {
  alternative <- match.arg(alternative)
  if (n1 < 2 || sd1 < 0) stop("Invalid first sample.", call. = FALSE)
  alpha <- 1 - level
  if (is.null(n2) && is.null(sd2)) {
    if (is.null(sd0) || sd0 <= 0) stop("Supply positive `sd0` for a one-sample variance test.", call. = FALSE)
    df <- n1 - 1; stat <- df * sd1^2 / sd0^2
    p <- switch(alternative, two.sided = min(1, 2 * min(stats::pchisq(stat, df), stats::pchisq(stat, df, lower.tail = FALSE))), less = stats::pchisq(stat, df), greater = stats::pchisq(stat, df, lower.tail = FALSE))
    ci_var <- c(df * sd1^2 / stats::qchisq(1 - alpha / 2, df), df * sd1^2 / stats::qchisq(alpha / 2, df))
    tab <- data.frame(n = n1, SD = .r4vn_num(sd1, digits), Variance = .r4vn_num(sd1^2, digits), Null.SD = .r4vn_num(sd0, digits), Chi.square = .r4vn_num(stat, digits), df = df, p = .r4vn_p(p, p_digits), SD.CI = paste0(.r4vn_num(sqrt(ci_var[1]), digits), " to ", .r4vn_num(sqrt(ci_var[2]), digits)), stringsAsFactors = FALSE)
    return(.r4vn_show(.r4vn_result("One-sample variance test", list("Result" = tab), raw = list(statistic = stat, df = df, p.value = p, conf.int.variance = ci_var), call = match.call()), show))
  }
  if (is.null(n2) || is.null(sd2) || n2 < 2 || sd2 < 0) stop("Supply valid `n2` and `sd2`.", call. = FALSE)
  f <- sd1^2 / sd2^2; df1 <- n1 - 1; df2 <- n2 - 1
  p <- switch(alternative, two.sided = min(1, 2 * min(stats::pf(f, df1, df2), stats::pf(f, df1, df2, lower.tail = FALSE))), less = stats::pf(f, df1, df2), greater = stats::pf(f, df1, df2, lower.tail = FALSE))
  ci <- c(f / stats::qf(1 - alpha / 2, df1, df2), f / stats::qf(alpha / 2, df1, df2))
  desc <- data.frame(Group = c("Group 1", "Group 2"), n = c(n1, n2), SD = .r4vn_num(c(sd1, sd2), digits), Variance = .r4vn_num(c(sd1^2, sd2^2), digits), stringsAsFactors = FALSE)
  tst <- data.frame(Variance.ratio = .r4vn_num(f, digits), CI = paste0(.r4vn_num(ci[1], digits), " to ", .r4vn_num(ci[2], digits)), F = .r4vn_num(f, digits), df1 = df1, df2 = df2, p = .r4vn_p(p, p_digits), stringsAsFactors = FALSE)
  .r4vn_show(.r4vn_result("Two-sample variance ratio test", list("Descriptive statistics" = desc, "Test" = tst), raw = list(estimate = f, conf.int = ci, statistic = f, df = c(df1, df2), p.value = p), call = match.call()), show)
}

# ============================================================================
# Function source: anovai.R
# ============================================================================
#' @rdname anova_r4vn
#' @export
anovai <- function(..., group.names = NULL, bartlett = TRUE, level = 0.95,
                   digits = 3, p_digits = 3, show = TRUE, console = FALSE) {
  groups <- lapply(list(...), .r4vn_summary_group); mat <- do.call(rbind, groups); colnames(mat) <- c("n", "mean", "sd")
  if (nrow(mat) < 2L || any(mat[, "n"] < 2) || any(mat[, "sd"] < 0)) stop("At least two valid groups are required.", call. = FALSE)
  k <- nrow(mat); N <- sum(mat[, "n"]); grand <- sum(mat[, "n"] * mat[, "mean"]) / N
  ssb <- sum(mat[, "n"] * (mat[, "mean"] - grand)^2); ssw <- sum((mat[, "n"] - 1) * mat[, "sd"]^2)
  dfb <- k - 1; dfw <- N - k; msb <- ssb / dfb; msw <- ssw / dfw; f <- msb / msw; p <- stats::pf(f, dfb, dfw, lower.tail = FALSE)
  eta2 <- ssb / (ssb + ssw); omega2 <- (ssb - dfb * msw) / (ssb + ssw + msw)
  gn <- group.names %||% paste0("Group ", seq_len(k)); crits <- stats::qt(1 - (1 - level) / 2, mat[, "n"] - 1); ses <- mat[, "sd"] / sqrt(mat[, "n"])
  desc <- data.frame(Group = gn, n = mat[, "n"], Mean = .r4vn_num(mat[, "mean"], digits), SD = .r4vn_num(mat[, "sd"], digits), SE = .r4vn_num(ses, digits), CI = paste0(.r4vn_num(mat[, "mean"] - crits * ses, digits), " to ", .r4vn_num(mat[, "mean"] + crits * ses, digits)), stringsAsFactors = FALSE)
  av <- data.frame(Source = c("Between", "Within", "Total"), SS = .r4vn_num(c(ssb, ssw, ssb + ssw), digits), df = c(dfb, dfw, N - 1), MS = c(.r4vn_num(msb, digits), .r4vn_num(msw, digits), ""), F = c(.r4vn_num(f, digits), "", ""), p = c(.r4vn_p(p, p_digits), "", ""), stringsAsFactors = FALSE)
  eff <- data.frame(Measure = c("Eta-squared", "Omega-squared"), Estimate = .r4vn_num(c(eta2, omega2), digits), stringsAsFactors = FALSE)
  sections <- list("Group statistics" = desc, "ANOVA" = av, "Effect size" = eff); braw <- NULL
  if (isTRUE(bartlett)) {
    vi <- mat[, "sd"]^2; ni <- mat[, "n"]; sp2 <- sum((ni - 1) * vi) / sum(ni - 1)
    if (all(vi > 0) && sp2 > 0) {
      num <- sum(ni - 1) * log(sp2) - sum((ni - 1) * log(vi)); corr <- 1 + (sum(1 / (ni - 1)) - 1 / sum(ni - 1)) / (3 * (k - 1)); bk <- num / corr; bp <- stats::pchisq(bk, k - 1, lower.tail = FALSE)
    } else { bk <- bp <- NA_real_ }
    sections[["Bartlett test of equal variances"]] <- data.frame(Chi.square = .r4vn_num(bk, digits), df = k - 1, p = .r4vn_p(bp, p_digits), stringsAsFactors = FALSE); braw <- list(statistic = bk, df = k - 1, p.value = bp)
  }
  raw <- list(groups = mat, grand.mean = grand, anova = list(ssb = ssb, ssw = ssw, df = c(dfb, dfw), F = f, p.value = p, eta2 = eta2, omega2 = omega2), bartlett = braw)
  .r4vn_show(.r4vn_result("One-way ANOVA from summary statistics", sections, raw = raw, call = match.call()), show)
}

# ============================================================================
# Function source: propi.R
# ============================================================================
#' @rdname proportion_tests
#' @export
propi <- function(events, total, 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) {
  method <- match.arg(method); alternative <- match.arg(alternative); .r4vn_check_counts(c(events, total)); if (total <= 0 || events > total) stop("Require 0 <= events <= total and total > 0.", call. = FALSE)
  methods <- if (method == "all") c("exact", "wilson", "wald") else method; cis <- lapply(methods, function(m) .r4vn_prop_ci(events, total, level, m))
  tab <- do.call(rbind, lapply(seq_along(methods), function(i) data.frame(Method = tools::toTitleCase(methods[i]), Proportion = .r4vn_num(events / total, digits), Lower = .r4vn_num(cis[[i]][1], digits), Upper = .r4vn_num(cis[[i]][2], digits), stringsAsFactors = FALSE)))
  sections <- list("Estimate" = data.frame(Events = events, Total = total, Proportion = .r4vn_num(events / total, digits), Percent = .r4vn_num(100 * events / total, 1), stringsAsFactors = FALSE), "Confidence intervals" = tab)
  bt <- NULL
  if (!is.null(p0)) { bt <- stats::binom.test(events, total, p = p0, alternative = alternative, conf.level = level); sections[["Exact binomial test"]] <- data.frame(Null.proportion = .r4vn_num(p0, digits), p = .r4vn_p(bt$p.value, p_digits), stringsAsFactors = FALSE) }
  .r4vn_show(.r4vn_result("Proportion from summary counts", sections, raw = list(events = events, total = total, proportion = events / total, conf.int = cis, test = bt), call = match.call()), show)
}

# ============================================================================
# Function source: bitesti.R
# ============================================================================
#' @rdname proportion_tests
#' @export
bitesti <- function(total, events, p = 0.5, alternative = c("two.sided", "less", "greater"), level = 0.95,
                    digits = 3, p_digits = 3, show = TRUE, console = FALSE) {
  call <- match.call(); alternative <- match.arg(alternative)
  out <- propi(events, total, p0 = p, method = "exact", alternative = alternative,
               level = level, digits = digits, p_digits = p_digits,
               show = FALSE, console = FALSE)
  out$call <- call
  .r4vn_show(out, show = show, console = console)
}

# ============================================================================
# Function source: prtesti.R
# ============================================================================
#' @rdname proportion_tests
#' @export
prtesti <- function(events1, n1, events2 = NULL, n2 = NULL, p0 = NULL,
                    alternative = c("two.sided", "less", "greater"), level = 0.95,
                    digits = 3, p_digits = 3, show = TRUE, console = FALSE) {
  alternative <- match.arg(alternative); .r4vn_check_counts(c(events1, n1)); if (n1 <= 0 || events1 > n1) stop("Invalid first proportion.", call. = FALSE); p1 <- events1 / n1; zcrit <- stats::qnorm(1 - (1 - level) / 2)
  if (is.null(events2) && is.null(n2)) {
    if (is.null(p0) || p0 <= 0 || p0 >= 1) stop("Supply `p0` between 0 and 1.", call. = FALSE)
    se0 <- sqrt(p0 * (1 - p0) / n1); z <- (p1 - p0) / se0; p <- switch(alternative, two.sided = 2 * stats::pnorm(abs(z), lower.tail = FALSE), less = stats::pnorm(z), greater = stats::pnorm(z, lower.tail = FALSE)); ci <- p1 + c(-1, 1) * zcrit * sqrt(p1 * (1 - p1) / n1)
    tab <- data.frame(Proportion = .r4vn_num(p1, digits), Null = .r4vn_num(p0, digits), Difference = .r4vn_num(p1 - p0, digits), CI = paste0(.r4vn_num(max(0, ci[1]), digits), " to ", .r4vn_num(min(1, ci[2]), digits)), z = .r4vn_num(z, digits), p = .r4vn_p(p, p_digits), stringsAsFactors = FALSE)
    return(.r4vn_show(.r4vn_result("One-sample proportion z test", list("Test" = tab), raw = list(estimate = p1, difference = p1 - p0, z = z, p.value = p, conf.int = ci), call = match.call()), show))
  }
  if (xor(is.null(events2), is.null(n2))) stop("Supply both `events2` and `n2` for a two-sample test.", call. = FALSE)
  .r4vn_check_counts(c(events2, n2)); if (n2 <= 0 || events2 > n2) stop("Invalid second proportion.", call. = FALSE); p2 <- events2 / n2; diff <- p1 - p2; pooled <- (events1 + events2) / (n1 + n2); se0 <- sqrt(pooled * (1 - pooled) * (1 / n1 + 1 / n2)); z <- diff / se0
  p <- switch(alternative, two.sided = 2 * stats::pnorm(abs(z), lower.tail = FALSE), less = stats::pnorm(z), greater = stats::pnorm(z, lower.tail = FALSE)); seci <- sqrt(p1 * (1 - p1) / n1 + p2 * (1 - p2) / n2); ci <- diff + c(-1, 1) * zcrit * seci
  desc <- data.frame(Group = c("Group 1", "Group 2"), Events = c(events1, events2), Total = c(n1, n2), Proportion = .r4vn_num(c(p1, p2), digits), stringsAsFactors = FALSE)
  tst <- data.frame(Difference = .r4vn_num(diff, digits), CI = paste0(.r4vn_num(ci[1], digits), " to ", .r4vn_num(ci[2], digits)), z = .r4vn_num(z, digits), p = .r4vn_p(p, p_digits), stringsAsFactors = FALSE)
  .r4vn_show(.r4vn_result("Two-sample proportion z test", list("Group proportions" = desc, "Test" = tst), raw = list(proportions = c(p1, p2), difference = diff, z = z, p.value = p, conf.int = ci), call = match.call()), show)
}

# ============================================================================
# Function source: cii.R
# ============================================================================
#' @rdname ci_r4vn
#' @export
cii <- function(n, mean = NULL, sd = NULL, events = NULL, variance = NULL,
                type = c("auto", "mean", "proportion", "variance"), method = c("exact", "wilson", "wald"),
                level = 0.95, digits = 3, show = TRUE, console = FALSE) {
  type <- match.arg(type); method <- match.arg(method)
  if (type == "auto") type <- if (!is.null(events)) "proportion" else if (!is.null(variance)) "variance" else "mean"
  if (type == "proportion") {
    call <- match.call()
    out <- propi(events, n, method = method, level = level, digits = digits, show = FALSE, console = FALSE)
    out$call <- call
    return(.r4vn_show(out, show = show, console = console))
  }
  if (type == "mean") {
    if (is.null(mean) || is.null(sd) || n < 2 || sd < 0) stop("Supply valid `n`, `mean`, and `sd`.", call. = FALSE); se <- sd / sqrt(n); crit <- stats::qt(1 - (1 - level) / 2, n - 1); ci <- mean + c(-1, 1) * crit * se
    tab <- data.frame(n = n, Mean = .r4vn_num(mean, digits), SE = .r4vn_num(se, digits), SD = .r4vn_num(sd, digits), Lower = .r4vn_num(ci[1], digits), Upper = .r4vn_num(ci[2], digits), stringsAsFactors = FALSE)
    return(.r4vn_show(.r4vn_result("Confidence interval for a mean", list("Estimate" = tab), raw = list(estimate = mean, se = se, conf.int = ci), call = match.call()), show))
  }
  if (is.null(variance) || variance < 0 || n < 2) stop("Supply valid `n` and `variance`.", call. = FALSE); df <- n - 1; alpha <- 1 - level; ci <- c(df * variance / stats::qchisq(1 - alpha / 2, df), df * variance / stats::qchisq(alpha / 2, df))
  tab <- data.frame(n = n, Variance = .r4vn_num(variance, digits), SD = .r4vn_num(sqrt(variance), digits), Variance.CI = paste0(.r4vn_num(ci[1], digits), " to ", .r4vn_num(ci[2], digits)), SD.CI = paste0(.r4vn_num(sqrt(ci[1]), digits), " to ", .r4vn_num(sqrt(ci[2]), digits)), stringsAsFactors = FALSE)
  .r4vn_show(.r4vn_result("Confidence interval for variance", list("Estimate" = tab), raw = list(estimate = variance, conf.int = ci), call = match.call()), show)
}

# ============================================================================
# Function source: esizei.R
# ============================================================================
#' @rdname esize
#' @export
esizei <- function(n1, mean1, sd1, n2, mean2, sd2, level = 0.95, digits = 3, show = TRUE, console = FALSE) {
  if (any(c(n1, n2) < 2) || any(c(sd1, sd2) < 0)) stop("Invalid group summaries.", call. = FALSE); df <- n1 + n2 - 2; sp <- sqrt(((n1 - 1) * sd1^2 + (n2 - 1) * sd2^2) / df); d <- (mean1 - mean2) / sp; J <- 1 - 3 / (4 * df - 1); g <- J * d; glass1 <- (mean1 - mean2) / sd1; glass2 <- (mean1 - mean2) / sd2
  se_d <- sqrt((n1 + n2) / (n1 * n2) + d^2 / (2 * df)); z <- stats::qnorm(1 - (1 - level) / 2); ci_d <- d + c(-1, 1) * z * se_d; ci_g <- J * ci_d
  tab <- data.frame(Measure = c("Cohen d", "Hedges g", "Glass delta (SD1)", "Glass delta (SD2)"), Estimate = .r4vn_num(c(d, g, glass1, glass2), digits), Lower = c(.r4vn_num(ci_d[1], digits), .r4vn_num(ci_g[1], digits), "", ""), Upper = c(.r4vn_num(ci_d[2], digits), .r4vn_num(ci_g[2], digits), "", ""), stringsAsFactors = FALSE)
  .r4vn_show(.r4vn_result("Effect sizes from two group summaries", list("Standardized effects" = tab), raw = list(cohen_d = d, hedges_g = g, glass = c(glass1, glass2), conf.int.d = ci_d, conf.int.g = ci_g), call = match.call()), show)
}

# ============================================================================
# Function source: ztesti.R
# ============================================================================
#' @rdname ztest
#' @export
ztesti <- function(n1, mean1, sd1, n2 = NULL, mean2 = NULL, sd2 = NULL, mu = 0,
                   alternative = c("two.sided", "less", "greater"), level = 0.95,
                   digits = 3, p_digits = 3, show = TRUE, console = FALSE) {
  alternative <- match.arg(alternative)
  second_supplied <- c(!is.null(n2), !is.null(mean2), !is.null(sd2))
  if (any(second_supplied) && !all(second_supplied)) stop("Supply `n2`, `mean2`, and `sd2` together.", call. = FALSE)
  one <- !any(second_supplied)
  if (!is.numeric(n1) || length(n1) != 1L || !is.finite(n1) || n1 <= 0 || !is.numeric(sd1) || length(sd1) != 1L || !is.finite(sd1) || sd1 <= 0) stop("Supply valid positive `n1` and `sd1`.", call. = FALSE)
  if (!one && (!is.numeric(n2) || length(n2) != 1L || !is.finite(n2) || n2 <= 0 || !is.numeric(sd2) || length(sd2) != 1L || !is.finite(sd2) || sd2 <= 0)) stop("Supply valid positive `n2` and `sd2`.", call. = FALSE)
  diff <- if (one) mean1 - mu else mean1 - mean2; se <- if (one) sd1 / sqrt(n1) else sqrt(sd1^2 / n1 + sd2^2 / n2); z <- diff / se
  p <- switch(alternative, two.sided = 2 * stats::pnorm(abs(z), lower.tail = FALSE), less = stats::pnorm(z), greater = stats::pnorm(z, lower.tail = FALSE)); crit <- stats::qnorm(1 - (1 - level) / 2); ci <- diff + c(-1, 1) * crit * se
  tab <- data.frame(Difference = .r4vn_num(diff, digits), SE = .r4vn_num(se, digits), CI = paste0(.r4vn_num(ci[1], digits), " to ", .r4vn_num(ci[2], digits)), z = .r4vn_num(z, digits), p = .r4vn_p(p, p_digits), stringsAsFactors = FALSE)
  .r4vn_show(.r4vn_result(if (one) "One-sample z test" else "Two-sample z test", list("Test" = tab), raw = list(difference = diff, se = se, z = z, p.value = p, conf.int = ci), call = match.call()), show)
}

# ============================================================================
# Function source: epii.R
# ============================================================================
#' @rdname epi
#' @export
epii <- function(a = NULL, b = NULL, c = NULL, d = NULL, by = NULL, level = 0.95,
                 correction = 0.5, digits = 3, p_digits = 3, show = TRUE, console = FALSE) {
  if (!is.null(by)) {
    strata <- .r4vn_strata(by); if (length(strata) < 2L) stop("Stratified analysis requires at least two strata.", call. = FALSE); m <- Reduce(`+`, strata)
    if (!all(vapply(list(a,b,c,d), is.null, logical(1)))) { supplied <- .r4vn_as_2x2(c(a,b,c,d)); if (!isTRUE(all.equal(unname(supplied), unname(m)))) warning("Supplied crude table differs from the sum of strata; the summed table is used.", call. = FALSE) }
  } else {
    strata <- NULL
    if (is.matrix(a) || is.table(a) || (is.numeric(a) && length(a) == 4L && all(vapply(list(b,c,d), is.null, logical(1))))) m <- .r4vn_as_2x2(a) else m <- .r4vn_as_2x2(c(a,b,c,d))
  }
  s <- .r4vn_epi_single(m, level, correction)
  display <- cbind(m, Total = rowSums(m)); display <- rbind(display, Total = c(colSums(m), sum(m)))
  measures <- data.frame(Measure = c("Risk/prevalence exposed", "Risk/prevalence unexposed", "Odds ratio", "Risk ratio", "Prevalence ratio", "Risk difference", if (s$RR >= 1) "Attributable fraction exposed" else "Prevented fraction exposed", if (s$RD > 0) "Number needed to harm" else if (s$RD < 0) "Number needed to treat" else "Number needed", if (s$PAR >= 0) "Population attributable risk" else "Population prevented risk", if (s$PAF >= 0) "Population attributable fraction" else "Population prevented fraction"), Estimate = .r4vn_num(c(s$risk.exposed,s$risk.unexposed,s$OR,s$RR,s$PR,s$RD,s$AF.exposed,s$NNT,abs(s$PAR),abs(s$PAF)), digits), Lower = c("","",.r4vn_num(s$OR.ci[1],digits),.r4vn_num(s$RR.ci[1],digits),.r4vn_num(s$PR.ci[1],digits),.r4vn_num(s$RD.ci[1],digits),"","","",""), Upper = c("","",.r4vn_num(s$OR.ci[2],digits),.r4vn_num(s$RR.ci[2],digits),.r4vn_num(s$PR.ci[2],digits),.r4vn_num(s$RD.ci[2],digits),"","","",""), stringsAsFactors = FALSE)
  tests <- data.frame(Test = c("Pearson chi-square", "Fisher exact"), Statistic = c(.r4vn_num(s$chi.square, digits), ""), p = c(.r4vn_p(s$chi.p,p_digits), .r4vn_p(s$fisher.p,p_digits)), stringsAsFactors = FALSE)
  sections <- list("2 x 2 table" = display, "Epidemiological measures" = measures, "Association tests" = tests); strat_raw <- NULL
  if (length(strata)) {
    sr <- lapply(strata, .r4vn_epi_single, level = level, correction = correction)
    strata_tab <- do.call(rbind, lapply(seq_along(sr), function(i) data.frame(Stratum = names(sr)[i], OR = .r4vn_ci(sr[[i]]$OR,sr[[i]]$OR.ci[1],sr[[i]]$OR.ci[2],digits), RR = .r4vn_ci(sr[[i]]$RR,sr[[i]]$RR.ci[1],sr[[i]]$RR.ci[2],digits), RD = .r4vn_ci(sr[[i]]$RD,sr[[i]]$RD.ci[1],sr[[i]]$RD.ci[2],digits), stringsAsFactors = FALSE)))
    arr <- array(0, dim = c(2,2,length(strata)), dimnames = list(c("Exposed","Unexposed"),c("Case","Noncase"),names(strata))); for (i in seq_along(strata)) arr[,,i] <- strata[[i]]
    mh <- tryCatch(stats::mantelhaen.test(arr, correct = FALSE, exact = FALSE), error = function(e) NULL)
    ai <- vapply(strata,function(x)x[1,1],numeric(1)); bi <- vapply(strata,function(x)x[1,2],numeric(1)); ci0 <- vapply(strata,function(x)x[2,1],numeric(1)); di <- vapply(strata,function(x)x[2,2],numeric(1)); n1i <- ai+bi; n0i <- ci0+di; ni <- n1i+n0i
    rr_mh <- sum(ai*n0i/ni)/sum(ci0*n1i/ni); varlog <- 1/pmax(ai,correction)-1/pmax(n1i,correction)+1/pmax(ci0,correction)-1/pmax(n0i,correction); w <- 1/varlog; se_pool <- sqrt(1/sum(w)); z <- stats::qnorm(1-(1-level)/2); rr_ci <- exp(log(rr_mh)+c(-1,1)*z*se_pool)
    logor <- log(vapply(sr,function(x) if (is.finite(x$OR)&&x$OR>0) x$OR else NA_real_,numeric(1))); vor <- vapply(strata,function(x){q<-if(any(x==0))x+correction else x;sum(1/q)},numeric(1)); good <- is.finite(logor)&is.finite(vor)&vor>0; q <- if(sum(good)>1){ww<-1/vor[good];sum(ww*(logor[good]-sum(ww*logor[good])/sum(ww))^2)}else NA_real_; qp <- if(is.na(q))NA_real_ else stats::pchisq(q,sum(good)-1,lower.tail=FALSE)
    gd <- do.call(rbind,lapply(seq_along(strata),function(i){x<-strata[[i]];data.frame(case=c(x[1,1],x[2,1]),noncase=c(x[1,2],x[2,2]),exposure=c(1,0),stratum=factor(names(strata)[i],levels=names(strata)))})); fit0 <- tryCatch(stats::glm(cbind(case,noncase)~exposure+stratum,family=stats::binomial(),data=gd),error=function(e)NULL); fit1 <- tryCatch(stats::glm(cbind(case,noncase)~exposure*stratum,family=stats::binomial(),data=gd),error=function(e)NULL); ip <- if(is.null(fit0)||is.null(fit1))NA_real_ else suppressWarnings(stats::anova(fit0,fit1,test="Chisq")$`Pr(>Chi)`[2])
    mh_or <- if (is.null(mh)) NA_real_ else unname(mh$estimate)
    or_change_mh <- if (is.finite(s$OR) && is.finite(mh_or) && mh_or != 0) 100 * (s$OR - mh_or) / mh_or else NA_real_
    or_change_crude <- if (is.finite(s$OR) && s$OR != 0 && is.finite(mh_or)) 100 * (s$OR - mh_or) / s$OR else NA_real_
    common <- data.frame(
      Measure = c("Crude OR", "Mantel-Haenszel common OR",
                  "% difference (crude - MH) / MH", "% difference (crude - MH) / crude",
                  "Mantel-Haenszel pooled RR (approx. CI)", "Mantel-Haenszel association test",
                  "Woolf homogeneity test", "Exposure-by-stratum interaction (LR)"),
      Estimate = c(.r4vn_num(s$OR,digits), if(is.null(mh))"" else .r4vn_num(mh_or,digits),
                   if(is.finite(or_change_mh)) paste0(.r4vn_num(or_change_mh,digits), "%") else "",
                   if(is.finite(or_change_crude)) paste0(.r4vn_num(or_change_crude,digits), "%") else "",
                   .r4vn_num(rr_mh,digits), if(is.null(mh))"" else .r4vn_num(unname(mh$statistic),digits),
                   .r4vn_num(q,digits), ""),
      Lower = c(.r4vn_num(s$OR.ci[1],digits), if(is.null(mh))"" else .r4vn_num(mh$conf.int[1],digits), "", "", .r4vn_num(rr_ci[1],digits), "", "", ""),
      Upper = c(.r4vn_num(s$OR.ci[2],digits), if(is.null(mh))"" else .r4vn_num(mh$conf.int[2],digits), "", "", .r4vn_num(rr_ci[2],digits), "", "", ""),
      p = c("", "", "", "", "", if(is.null(mh))"" else .r4vn_p(mh$p.value,p_digits), .r4vn_p(qp,p_digits), .r4vn_p(ip,p_digits)),
      stringsAsFactors = FALSE
    )
    sections[["Stratum-specific estimates"]] <- strata_tab; sections[["Adjusted and interaction analysis"]] <- common; strat_raw <- list(strata=sr,mantel.haenszel=mh,rr.mh=rr_mh,rr.mh.ci=rr_ci,or.crude=s$OR,or.mh=mh_or,or.percent.diff.mh=or_change_mh,or.percent.diff.crude=or_change_crude,homogeneity=list(statistic=q,p.value=qp),interaction.p=ip)
  }
  notes <- if (s$zero.correction) paste0("A continuity correction of ", correction, " was used only for log-scale confidence intervals because at least one cell was zero.") else NULL
  .r4vn_show(.r4vn_result(if (length(strata)) "Stratified epidemiological 2 x 2 analysis" else "Epidemiological 2 x 2 analysis", sections, notes, raw = list(crude=s,stratified=strat_raw,table=m), call = match.call()), show)
}

# ============================================================================
# Function source: cci.R
# ============================================================================
#' @rdname epi
#' @export
cci <- function(a, b, c, d, ...) epii(a, b, c, d, ...)

# ============================================================================
# Function source: csi.R
# ============================================================================
#' @rdname epi
#' @export
csi <- function(a, b, c, d, ...) epii(a, b, c, d, ...)

# ============================================================================
# Function source: mcci.R
# ============================================================================
#' @rdname mcc
#' @export
mcci <- function(a, b, c, d, level = 0.95, digits = 3, p_digits = 3, show = TRUE, console = FALSE) {
  .r4vn_check_counts(c(a,b,c,d)); disc <- b+c; or <- b/c; bt <- if(disc>0) stats::binom.test(b,disc,p=.5,conf.level=level) else NULL; pci <- if(is.null(bt)) c(NA,NA) else unname(bt$conf.int); ci <- pci/(1-pci)
  tab <- matrix(c(a,b,c,d),2,2,byrow=TRUE,dimnames=list("Case"=c("Exposed","Unexposed"),"Control"=c("Exposed","Unexposed")))
  res <- data.frame(Discordant.case.exposed=b,Discordant.control.exposed=c,Matched.OR=.r4vn_num(or,digits),Lower=.r4vn_num(ci[1],digits),Upper=.r4vn_num(ci[2],digits),Exact.p=if(is.null(bt))"" else .r4vn_p(bt$p.value,p_digits),stringsAsFactors=FALSE)
  .r4vn_show(.r4vn_result("Matched case-control analysis",list("Matched-pair table"=tab,"Matched estimate"=res),raw=list(OR=or,conf.int=ci,p.value=if(is.null(bt))NA else bt$p.value),call=match.call()),show)
}

# ============================================================================
# Function source: iri.R
# ============================================================================
#' @rdname ir
#' @export
iri <- function(cases.exposed, cases.unexposed, time.exposed, time.unexposed,
                level = 0.95, digits = 4, p_digits = 3, show = TRUE, console = FALSE) {
  .r4vn_check_counts(c(cases.exposed, cases.unexposed), "case counts")
  if(any(c(time.exposed,time.unexposed)<=0) || any(!is.finite(c(time.exposed,time.unexposed)))) stop("Person-time must be finite and positive.",call.=FALSE)
  r1<-cases.exposed/time.exposed;r0<-cases.unexposed/time.unexposed;rd<-r1-r0;rr<-r1/r0;pt<-tryCatch(stats::poisson.test(c(cases.exposed,cases.unexposed),T=c(time.exposed,time.unexposed),conf.level=level),error=function(e)NULL);ci<-if(is.null(pt))c(NA,NA) else unname(pt$conf.int);z<-stats::qnorm(1-(1-level)/2);se_rd<-sqrt(cases.exposed/time.exposed^2+cases.unexposed/time.unexposed^2);ci_rd<-rd+c(-1,1)*z*se_rd;afe<-if(is.finite(rr)&&rr>=1)(rr-1)/rr else if(is.finite(rr))1-rr else NA;afp<-if((cases.exposed+cases.unexposed)>0)(cases.exposed/(cases.exposed+cases.unexposed))*afe else NA
  desc<-data.frame(Group=c("Exposed","Unexposed","Total"),Cases=c(cases.exposed,cases.unexposed,cases.exposed+cases.unexposed),Person.time=c(time.exposed,time.unexposed,time.exposed+time.unexposed),Rate=.r4vn_num(c(r1,r0,(cases.exposed+cases.unexposed)/(time.exposed+time.unexposed)),digits),stringsAsFactors=FALSE)
  meas<-data.frame(Measure=c("Incidence-rate difference","Incidence-rate ratio",if(rr>=1)"Attributable fraction exposed" else "Prevented fraction exposed","Population attributable fraction"),Estimate=.r4vn_num(c(rd,rr,afe,afp),digits),Lower=c(.r4vn_num(ci_rd[1],digits),.r4vn_num(ci[1],digits),"",""),Upper=c(.r4vn_num(ci_rd[2],digits),.r4vn_num(ci[2],digits),"",""),p=c("",if(is.null(pt))"" else .r4vn_p(pt$p.value,p_digits),"",""),stringsAsFactors=FALSE)
  .r4vn_show(.r4vn_result("Incidence-rate comparison",list("Rates"=desc,"Rate measures"=meas),raw=list(rate=c(r1,r0),difference=rd,ratio=rr,ratio.ci=ci,p.value=if(is.null(pt))NA else pt$p.value),call=match.call()),show)
}



# ============================================================================
# Viewer renderers for immediate commands
# Each public command has its own entry point so its layout can be customized
# later without changing the computational code or another command's Viewer.
# ============================================================================
.r4vn_immediate_viewer <- function(x, subtitle = "Immediate analysis from summary data") {
  blocks <- paste0(vapply(names(x$sections), function(nm) {
    obj <- x$sections[[nm]]
    .r4vn_view_section(nm, obj)
  }, character(1)), collapse = "")
  .r4vn_view_document(x$title, blocks, notes = x$notes, subtitle = subtitle, prefix = "r4vn-immediate-")
}

.r4vn_viewer_tabi <- function(x) .r4vn_immediate_viewer(x, "Contingency table from entered counts")
.r4vn_viewer_ttesti <- function(x) .r4vn_immediate_viewer(x, "t test from entered summary statistics")
.r4vn_viewer_sdtesti <- function(x) .r4vn_immediate_viewer(x, "Variance test from entered summary statistics")
.r4vn_viewer_anovai <- function(x) .r4vn_immediate_viewer(x, "One-way ANOVA from entered summary statistics")
.r4vn_viewer_propi <- function(x) .r4vn_immediate_viewer(x, "Proportion analysis from entered counts")
.r4vn_viewer_bitesti <- function(x) .r4vn_immediate_viewer(x, "Exact binomial test from entered counts")
.r4vn_viewer_prtesti <- function(x) .r4vn_immediate_viewer(x, "Proportion z test from entered counts")
.r4vn_viewer_cii <- function(x) .r4vn_immediate_viewer(x, "Confidence interval from summary statistics")
.r4vn_viewer_esizei <- function(x) .r4vn_immediate_viewer(x, "Effect size from entered group summaries")
.r4vn_viewer_ztesti <- function(x) .r4vn_immediate_viewer(x, "z test from entered summary statistics")
.r4vn_viewer_epii <- function(x) .r4vn_immediate_viewer(x, "Epidemiological analysis from a 2 x 2 table")
.r4vn_viewer_mcci <- function(x) .r4vn_immediate_viewer(x, "Matched case-control analysis from entered counts")
.r4vn_viewer_iri <- function(x) .r4vn_immediate_viewer(x, "Incidence-rate analysis from entered counts and person-time")