R/T_test.R

In autotestR: Automated Functions for Basic Statistical Tests

#' Student's t-test (with automatic diagnostics)
#'
#' Performs Student's t-test to compare means between two independent groups,
#' with automatic checks for normality and homogeneity of variances.
#' If assumptions are violated, the Mann-Whitney test is automatically applied
#' (without generating a plot).
#'
#' @param ... Two numeric vectors or a data frame with exactly two columns.
#' @param title Logical. If true, return a plot entitled.
#' @param title_text Plot title (string). Default: "t-test".
#' @param xlab X-axis label in the plot (string). Default: "Group".
#' @param ylab Y-axis label in the plot (string). Default: "Value".
#' @param style Plot aesthetic generated by the function.
#' @param help Logical. If TRUE, shows a detailed explanation of the function.
#'        Default: FALSE.
#' @param verbose Logical. If TRUE, prints detailed messages. Default: TRUE.
#' @return Invisible list with summary, test test_result, method and (optionally) plot.
#' @export
#'
#' @examples
#' set.seed(123)
#' df <- data.frame(
#'   control = rnorm(30, 10),
#'   treatment = rnorm(30, 15)
#' )
#' test.t(df)

test.t <- function(...,
                   title = TRUE,
                   title_text = "t-test" ,
                   xlab = "Group",
                   ylab = "Value",
                   style = c("boxplot", "violin", "mono", "halfeye"),
                   help = FALSE,
                   verbose = TRUE) {

  args <- list(...)
  style <- match.arg(style)

  # ============================
  # Quick help
  # ============================
  if (help) {
    if (verbose) {
      message("
Function test.t()

Description:
  Performs Student's t-test to compare the means of two independent groups,
  with automatic checks for normality (Shapiro-Wilk) and homogeneity of variances
  (Levene's test).

If assumptions are violated, the Mann-Whitney test is automatically applied
and no plot is generated.

Accepted inputs:
  - Two numeric vectors (e.g., group1, group2)
  - A data frame with exactly two numeric columns

Example:
  set.seed(123)
  df <- data.frame(
    control = rnorm(30, 10),
    treatment = rnorm(30, 15)
  )
  test.t(df)
")
    }
    return(invisible(NULL))
  }

  # ============================
  # Detect data frame input
  # ============================
  if (length(args) == 1 && is.data.frame(args[[1]])) {
    df <- args[[1]]
    if (ncol(df) != 2) {
      stop("The data frame must contain exactly two numeric columns.")
    }
    groups <- as.list(df)
    group_names <- colnames(df)
  } else {
    groups <- args
    if (length(groups) != 2) {
      stop("Provide exactly two groups or a data frame with two columns.")
    }
    group_names <- sapply(substitute(list(...))[-1], deparse)
  }

  # ============================
  # Validation
  # ============================
  if (!all(sapply(groups, is.numeric))) {
    stop("Both groups must be numeric.")
  }
  if (any(sapply(groups, function(g) sd(g, na.rm = TRUE) == 0))) {
    stop("One of the groups has zero variance (constant values).")
  }

  # ============================
  # Required packages
  # ============================
  required_packages <- c("ggplot2", "dplyr", "car", "ggdist")
  for (pkg in required_packages) {
    if (!requireNamespace(pkg, quietly = TRUE)) {
      stop(
        paste0(
          "Package '", pkg,
          "' is not installed. Install it with install.packages('", pkg, "')"
        )
      )
    }
  }

  # ============================
  # Build long-format data frame
  # ============================
  values <- unlist(groups)
  group <- factor(
    rep(group_names, times = sapply(groups, length)),
    levels = group_names
  )
  data <- data.frame(value = values, group = group)

  # ============================
  # Normality test (Shapiro-Wilk)
  # ============================
  normality <- sapply(groups, function(g) {
    if (length(g) < 3) return(NA)
    stats::shapiro.test(g[!is.na(g)])$p.value
  })
  is_normal <- all(normality > 0.05, na.rm = TRUE)

  # ============================
  # Homogeneity test (Levene)
  # ============================
  homogeneity <- tryCatch(
    car::leveneTest(value ~ group, data = data, na.action = na.omit),
    error = function(e) data.frame(`Pr(>F)` = NA)
  )
  is_homogeneous <- !is.na(homogeneity$`Pr(>F)`[1]) &&
    homogeneity$`Pr(>F)`[1] > 0.05

  # --- Assumption checks ---
  violated_normality   <- !is_normal
  violated_homogeneity <- !is_homogeneous

  assumption_warning <- violated_normality || violated_homogeneity

  # ============================
  # Student t-test
  # ============================
  test_result <- stats::t.test(groups[[1]], groups[[2]])

  method <- "Student's t-test"

  ci_low  <- test_result$conf.int[1]
  ci_high <- test_result$conf.int[2]

  p_value <- test_result$p.value
  p_label <- .format_pval(p_value)

  # ============================
  # Effect size: Cohen's d
  # ============================
  x <- groups[[1]]
  y <- groups[[2]]

  nx <- sum(!is.na(x))
  ny <- sum(!is.na(y))

  mean_diff <- mean(x, na.rm = TRUE) - mean(y, na.rm = TRUE)

  sd_pooled <- sqrt(
    ((nx - 1) * sd(x, na.rm = TRUE)^2 +
       (ny - 1) * sd(y, na.rm = TRUE)^2) /
      (nx + ny - 2)
  )

  cohen_d <- mean_diff / sd_pooled

  # ============================
  # Descriptive summary
  # ============================
  summary_table <- data |>
    dplyr::group_by(group) |>
    dplyr::summarise(
      Mean = round(mean(value, na.rm = TRUE), 2),
      SD = round(sd(value, na.rm = TRUE), 2),
      .groups = "drop"
    )

  if (verbose && assumption_warning) {

    .print_block("Assumption check", function() {

      if (violated_normality) {
        cat("Warning: Normality assumption violated (Shapiro-Wilk p < 0.05)\n")
      }

      if (violated_homogeneity) {
        cat("Warning: Variance homogeneity violated (Levene p < 0.05)\n")
      }

      cat("Interpret results with caution.\n")

    })
  }

  # ============================
  # Labels and colors
  # ============================
    p_label <- signif(p_value, 3)

    p_label_sub <- .format_p(p_value)

    signif_label <- if (p_value < 0.001) {
      "***"
    } else if (p_value < 0.01) {
      "**"
    } else if (p_value < 0.05) {
      "*"
    } else {
      ""
    }

    # --- Subtitle ---
    subtitle_text <- paste0(
      "diff = ", round(mean_diff, 2),
      " | p ", ifelse(p_value < 0.001, "", "= "), p_label_sub
    )

    # --- Colors ---
    y_pos <- max(values, na.rm = TRUE) +
      0.1 * diff(range(values, na.rm = TRUE))

    vivid_colors <- scales::hue_pal()(length(unique(data$group)))
    mono_colors <- c("grey75", "grey25")

  # ============================
  # STYLE 1 (Boxplot + jitter)
  # ============================
    if (style == "boxplot") {
      g <- ggplot2::ggplot(data, ggplot2::aes(x = group, y = value, fill = group)) +
        ggplot2::geom_boxplot(alpha = 0.75, outlier.shape = NA, width = 0.7, linewidth = 0.7) +
        ggplot2::geom_jitter(width = 0.1, alpha = 0.2, color = "grey25", size = 1.8) +
        ggplot2::annotate(
          "text", x = mean(1:2), y = y_pos,
          label = signif_label, size = 6,
          col = "grey25"
        ) +
        ggplot2::theme_minimal(base_size = 12) +
        ggplot2::scale_fill_manual(values = vivid_colors) +
        ggplot2::labs(
          title = if (title) title_text else NULL,
          subtitle = subtitle_text,
          x = "",
          y = ylab
        ) +
        ggplot2::theme(
          legend.position = "none",
          plot.margin = ggplot2::margin(5.5, 5.5, 10, 5.5),
          axis.text.x = ggplot2::element_text(
            angle = 45, hjust = 1, size = 12
          )
        )
    }

  # ============================
  # STYLE 2 (Clean violin)
  # ============================
    if (style == "violin") {
      g <- ggplot2::ggplot(data, ggplot2::aes(x = group, y = value, fill = group)) +
        ggplot2::geom_violin(
          trim = FALSE, alpha = .6,
          color = NA, adjust = .6
        ) +
        ggplot2::geom_boxplot(
          width = .18, outlier.shape = NA,
          color = "gray20", linewidth = .4
        ) +
        ggplot2::annotate(
          "text", x = mean(1:2), y = y_pos,
          label = signif_label, size = 6,
          col = "grey25"
        ) +
        ggplot2::geom_point(
          position = ggplot2::position_jitter(width = .1),
          alpha = .2, size = 1.8, color = "gray25"
        ) +
        ggplot2::scale_fill_manual(values = vivid_colors) +
        ggplot2::theme_minimal(base_size = 12) +
        ggplot2::labs(
          title = if (title) title_text else NULL,
          subtitle = subtitle_text,
          x = "",
          y = ylab
        ) +
        ggplot2::theme(
          legend.position = "none",
          axis.text.x = ggplot2::element_text(
            angle = 45, hjust = 1, size = 12
          )
        )
    }

  # ============================
  # STYLE 3 (Premium monochrome)
  # ============================
  if (style == "mono") {
    g <- ggplot2::ggplot(data, ggplot2::aes(x = group, y = value, fill = group)) +
      ggplot2::geom_boxplot(alpha = 0.75, outlier.shape = NA, width = 0.7, linewidth = 0.7) +
      ggplot2::geom_jitter(width = 0.1, alpha = 0.2, color = "grey25", size = 1.8) +
      ggplot2::annotate(
        "text", x = mean(1:2), y = y_pos,
        label = signif_label, size = 6,
        col = "grey25"
      ) +
      ggplot2::theme_minimal(base_size = 12) +
      ggplot2::scale_fill_manual(values = mono_colors) +
      ggplot2::labs(
        title = if (title) title_text else NULL,
        subtitle = subtitle_text,
        x = "",
        y = ylab
      ) +
      ggplot2::theme(
        legend.position = "none",
        plot.margin = ggplot2::margin(5.5, 5.5, 10, 5.5),
        axis.text.x = ggplot2::element_text(
          angle = 45, hjust = 1, size = 12
        )
      )
  }

  # ============================
  # STYLE 4 (ggdist half-eye + median)
  # ============================
    if (style == "halfeye") {

      if (!requireNamespace("ggdist", quietly = TRUE))
        stop("Package 'ggdist' is required for style = halfeye'.")

      g <- ggplot2::ggplot(
        data,
        ggplot2::aes(x = group, y = value, fill = group)
      ) +
        ggdist::stat_halfeye(
          alpha = 0.6,
          trim = FALSE,
          adjust = 0.6,
          width = 0.6,
          .width = c(0.5, 0.8, 0.95),
          justification = -0.2,
          slab_color = "gray20",
          interval_color = "gray20"
        ) +
        ggplot2::annotate(
          "text", x = mean(1:2), y = y_pos,
          label = signif_label, size = 6,
          col = "grey25"
        ) +
        ggplot2::geom_point(
          position = ggplot2::position_nudge(x = 0.15),
          size = 1.1, alpha = 0.4, color = "black"
        ) +
        ggdist::stat_pointinterval(
          position = ggplot2::position_nudge(x = 0.2),
          point_color = "black",
          interval_color = "black",
          .width = 0.95
        ) +
        ggplot2::theme_minimal(base_size = 12) +
        ggplot2::scale_fill_manual(values = vivid_colors) +
        ggplot2::labs(
          title = if (title) title_text else NULL,
          subtitle = subtitle_text,
          x = "",
          y = ylab
        ) +
        ggplot2::theme(
          legend.position = "none",
          axis.text.x = ggplot2::element_text(
            angle = 45, hjust = 1, size = 12
          )
        )
    }

    print(g)

  # ============================
  # Output
  # ============================
  obj <- list(
    summary = summary_table,
    result = test_result,
    method = method,
    normality = normality,
    homogeneity = homogeneity,
    plot = g
  )

  # ============================
  # Return
  # ============================
  if (verbose) {

    .print_header(method)

    .print_block("Summary", function() {
      print(summary_table, row.names = FALSE)
    })

    .print_block("Statistics", function() {

      cat("t statistic:          ", round(test_result$statistic, 3), "\n",
          "Degrees of freedon:   ", round(test_result$parameter, 1), "\n",
          "p-value:              ", p_label, "\n", sep = "")
      cat("Mean difference:      ", round(mean_diff, 2), " [", round(ci_low, 2), ", ", round(ci_high, 2), "]\n", sep = "")
      cat("Cohen's d:            ", round(cohen_d, 3), "\n", sep = "")

    })
  }

  return(invisible(list(result = obj)))

}