Group Comparison — t-test, ANOVA & Nonparametric

suppressPackageStartupMessages(library(DT))
suppressPackageStartupMessages(library(htmlwidgets))
suppressPackageStartupMessages(library(arrow))
suppressPackageStartupMessages(library(knitr))
suppressPackageStartupMessages(library(rmarkdown))
suppressPackageStartupMessages(library(dplyr))
suppressPackageStartupMessages(library(tidyr))
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(stringr))
suppressPackageStartupMessages(library(lubridate))
suppressPackageStartupMessages(library(broom))
suppressPackageStartupMessages(library(Matrix))
suppressPackageStartupMessages(library(cluster))
suppressPackageStartupMessages(library(data.table))

compute_shared <- function(df, params, col_map = list()) {
  # === SHARED EXPORTS ===
  #   initial_rows/final_rows/rows_removed  $ row accounting
  #   outcome_h / group_h  $ humanized user names for the two mapped columns
  #   k / group_levels     $ number of groups after cleaning + their names
  #   n_na_outcome         $ rows dropped for missing/non-numeric outcome
  #   dropped_groups_df    $ data.frame(group, n) — groups dropped (n < 3)
  #   lumped_levels        $ character — levels folded into "Other"
  #   group_summary_df     $ group, n, mean, sd, median, ci_low, ci_high
  #   test_results_df      $ test, statistic, p_value, effect_size, interpretation
  #   pairwise_df          $ comparison, difference, ci_low, ci_high, adj_p, significant
  #   boxplot_df           $ group_name, outcome_value (<= 2000 sampled rows)
  #   shapiro_p / out_frac / nonparam_preferred  $ assumption diagnostics
  #   primary_test / primary_p / primary_sig     $ the headline result
  #   effect_value / effect_text                 $ headline effect size
  #   top_sig_pair         $ 1-row slice of pairwise_df or NULL
  #   metrics / json_output
  # === /SHARED EXPORTS ===

initial_rows <- nrow(df)
  outcome_h <- humanize_semantic("outcome", col_map)
  group_h   <- humanize_semantic("group", col_map)
  if (!("outcome" %in% names(df)) || !("group" %in% names(df))) {
    stop(sprintf("Group comparison needs both &#x27;%s' (the numeric outcome) and '%s' (the groups) mapped.",
                 outcome_h, group_h))
  }

v <- df$outcome
  if (!is.numeric(v)) {
    ch <- as.character(v)
    non_blank <- !is.na(ch) & trimws(ch) != ""
    conv <- suppressWarnings(as.numeric(ch))
    if (sum(non_blank) == 0 ||
        sum(!is.na(conv[non_blank])) < 0.95 * sum(non_blank)) {
      stop(sprintf("The outcome column &#x27;%s' does not look numeric — fewer than 95%% of its values parse as numbers. Pick a numeric column to compare.",
                   outcome_h))
    }
    v <- conv
  }
  df$outcome <- v

  g <- as.character(df$group)
  g[is.na(g) | trimws(g) == ""] <- "Missing"

  keep <- !is.na(df$outcome)
  n_na_outcome <- sum(!keep)
  df <- df[keep, , drop = FALSE]
  g  <- g[keep]
  if (nrow(df) == 0) {
    stop(sprintf("No rows with a usable numeric value in &#x27;%s' remained after cleaning.", outcome_h))
  }

tab <- table(g)
  small <- names(tab)[tab < 3]
  dropped_groups_df <- data.frame(group = character(0), n = integer(0),
                                  stringsAsFactors = FALSE)
  if (length(small) > 0) {
    dropped_groups_df <- data.frame(group = small, n = as.integer(tab[small]),
                                    stringsAsFactors = FALSE)
    sel <- !(g %in% small)
    df <- df[sel, , drop = FALSE]
    g  <- g[sel]
  }

  lumped_levels <- character(0)
  tab <- sort(table(g), decreasing = TRUE)
  if (length(tab) > 8) {
    keep_lv <- names(tab)[1:8]
    lumped_levels <- setdiff(names(tab), keep_lv)
    g[g %in% lumped_levels] <- "Other"
  }

tab <- table(g)
  small2 <- names(tab)[tab < 3]
  if (length(small2) > 0) {
    dropped_groups_df <- rbind(dropped_groups_df,
                               data.frame(group = small2, n = as.integer(tab[small2]),
                                          stringsAsFactors = FALSE))
    sel <- !(g %in% small2)
    df <- df[sel, , drop = FALSE]
    g  <- g[sel]
  }

  gf <- factor(g)
  k  <- nlevels(gf)
  if (k < 2) {
    stop(sprintf("Group comparison needs at least 2 groups in &#x27;%s' with 3 or more rows each; only %d usable group(s) remained after cleaning. Check that '%s' really splits the data into groups.",
                 group_h, k, group_h))
  }
  y <- df$outcome
  final_rows <- length(y)
  rows_removed <- initial_rows - final_rows
  if (final_rows < 10) {
    stop(sprintf("Only %d usable rows remained — at least 10 are needed to compare groups.", final_rows))
  }
  if (isTRUE(stats::var(y) == 0)) {
    stop(sprintf("The outcome &#x27;%s' has no variation at all (every value is identical) — there is nothing to compare.", outcome_h))
  }
  group_levels <- levels(gf)

group_summary_df <- do.call(rbind, lapply(group_levels, function(l) {
    x <- y[gf == l]
    n <- length(x)
    m <- mean(x)
    s <- stats::sd(x)
    half <- if (!is.na(s) && s > 0 && n > 1) stats::qt(0.975, n - 1) * s / sqrt(n) else 0
    data.frame(group = l, n = n,
               mean = round(m, 3),
               sd = round(if (is.na(s)) 0 else s, 3),
               median = round(stats::median(x), 3),
               ci_low = round(m - half, 3),
               ci_high = round(m + half, 3),
               stringsAsFactors = FALSE)
  }))
  rownames(group_summary_df) <- NULL

dat <- data.frame(y = y, gf = gf)
  fit <- stats::aov(y ~ gf, data = dat)
  res <- stats::residuals(fit)
  shapiro_p <- NA_real_
  if (final_rows >= 3 && final_rows <= 5000) {
    shapiro_p <- tryCatch(stats::shapiro.test(res)$p.value,
                          error = function(e) NA_real_)
  }
  iqr <- stats::IQR(res)
  qs  <- stats::quantile(res, c(0.25, 0.75), names = FALSE)
  out_frac <- if (is.na(iqr) || iqr == 0) 0 else
    mean(res < qs[1] - 3 * iqr | res > qs[2] + 3 * iqr)
  heavy_outliers <- out_frac > 0.01
  nonparam_preferred <- (!is.na(shapiro_p) && shapiro_p < 0.01) || heavy_outliers

  sig_phrase <- function(p) {
    if (is.na(p)) "could not be computed"
    else if (p < 0.001) "highly significant(p < 0.001)"
    else if (p < 0.05) sprintf("significant(p = %.3g)", p)
    else sprintf("not significant(p = %.3g)", p)
  }
  d_word <- function(d) {
    a <- abs(d)
    if (is.na(a)) "unknown"
    else if (a < 0.2) "negligible" else if (a < 0.5) "small"
    else if (a < 0.8) "medium" else "large"
  }
  eta_word <- function(e) {
    if (is.na(e)) "unknown"
    else if (e < 0.01) "negligible" else if (e < 0.06) "small"
    else if (e < 0.14) "medium" else "large"
  }

test_rows <- list()
  add_test <- function(test, statistic, p_value, effect_size, interpretation) {
    test_rows[[length(test_rows) + 1]] <<- data.frame(
      test = test,
      statistic = round(statistic, 3),
      p_value = signif(p_value, 3),
      effect_size = if (is.na(effect_size)) NA_real_ else round(effect_size, 3),
      interpretation = interpretation,
      stringsAsFactors = FALSE
    )
  }

  pairwise_df <- NULL
  effect_value <- NA_real_
  effect_text <- ""
  primary_test <- NA_character_
  primary_p <- NA_real_

  if (k == 2) {

l1 <- group_levels[1]; l2 <- group_levels[2]
    x1 <- y[gf == l1]; x2 <- y[gf == l2]
    n1 <- length(x1); n2 <- length(x2)
    s1 <- stats::sd(x1); s2 <- stats::sd(x2)
    sp <- sqrt(((n1 - 1) * s1^2 + (n2 - 1) * s2^2) / (n1 + n2 - 2))
    cohens_d <- if (!is.na(sp) && sp > 0) (mean(x1) - mean(x2)) / sp else NA_real_
    higher <- if (mean(x1) >= mean(x2)) l1 else l2

    welch <- tryCatch(stats::t.test(x1, x2), error = function(e) NULL)
    student <- tryCatch(stats::t.test(x1, x2, var.equal = TRUE),
                        error = function(e) NULL)
    mw <- tryCatch(suppressWarnings(stats::wilcox.test(x1, x2)),
                   error = function(e) NULL)
    rank_biserial <- if (!is.null(mw)) 2 * unname(mw$statistic) / (n1 * n2) - 1 else NA_real_

    if (!is.null(welch)) {
      add_test("Welch t-test(unequal variances)", unname(welch$statistic),
               welch$p.value, cohens_d,
               sprintf("Difference in mean %s is %s; %s effect(Cohen&#x27;s d), %s mean is higher.",
                       outcome_h, sig_phrase(welch$p.value), d_word(cohens_d), higher))
    }
    if (!is.null(student)) {
      add_test("Student t-test(equal variances)", unname(student$statistic),
               student$p.value, cohens_d,
               sprintf("Classic t-test cross-check: %s.", sig_phrase(student$p.value)))
    }
    if (!is.null(mw)) {
      add_test("Mann-Whitney U", unname(mw$statistic), mw$p.value, rank_biserial,
               sprintf("Rank-based(no normality assumption): %s; rank-biserial r = %s.",
                       sig_phrase(mw$p.value),
                       ifelse(is.na(rank_biserial), "NA", round(rank_biserial, 3))))
    }
    if (!is.na(shapiro_p)) {
      add_test("Shapiro-Wilk normality(residuals)", NA_real_, shapiro_p, NA_real_,
               if (shapiro_p < 0.01)
                 "Normality clearly violated — favor the rank-based result."
               else "No strong evidence against normality.")
    }

dmean <- mean(x1) - mean(x2)
    if (!is.null(welch)) {
      ci <- as.numeric(welch$conf.int)
      pw_p <- welch$p.value
    } else if (!is.null(student)) {
      ci <- as.numeric(student$conf.int)
      pw_p <- student$p.value
    } else {
      ci <- c(dmean, dmean)
      pw_p <- NA_real_
    }
    pairwise_df <- data.frame(
      comparison = paste0(l1, " - ", l2),
      difference = round(dmean, 3),
      ci_low = round(ci[1], 3),
      ci_high = round(ci[2], 3),
      adj_p = signif(pw_p, 3),
      significant = ifelse(is.na(pw_p), "", ifelse(pw_p < 0.05, "yes", "no")),
      stringsAsFactors = FALSE
    )

    effect_value <- cohens_d
    effect_text <- sprintf("Cohen&#x27;s d = %s (%s effect)",
                           ifelse(is.na(cohens_d), "NA", round(cohens_d, 2)),
                           d_word(cohens_d))
    if (nonparam_preferred && !is.null(mw)) {
      primary_test <- "Mann-Whitney U"; primary_p <- mw$p.value
    } else if (!is.null(welch)) {
      primary_test <- "Welch t-test"; primary_p <- welch$p.value
    } else if (!is.null(student)) {
      primary_test <- "Student t-test"; primary_p <- student$p.value
    } else if (!is.null(mw)) {
      primary_test <- "Mann-Whitney U"; primary_p <- mw$p.value
    } else {
      stop(sprintf("No comparison test could be computed for &#x27;%s' between the two %s groups (the data may be essentially constant).",
                   outcome_h, group_h))
    }
  } else {

a_tab <- summary(fit)[[1]]
    ss_b <- a_tab[["Sum Sq"]][1]
    ss_w <- a_tab[["Sum Sq"]][2]
    eta2 <- if (!is.na(ss_b) && !is.na(ss_w) && (ss_b + ss_w) > 0)
      ss_b / (ss_b + ss_w) else NA_real_
    f_stat <- a_tab[["F value"]][1]
    anova_p <- a_tab[["Pr(>F)"]][1]

    welch_a <- tryCatch(stats::oneway.test(y ~ gf, data = dat),
                        error = function(e) NULL)
    kw <- tryCatch(stats::kruskal.test(y ~ gf, data = dat),
                   error = function(e) NULL)

    if (!is.na(anova_p)) {
      add_test("One-way ANOVA", f_stat, anova_p, eta2,
               sprintf("Mean %s differs across %s groups: %s; eta-squared = %s(%s effect).",
                       outcome_h, group_h, sig_phrase(anova_p),
                       ifelse(is.na(eta2), "NA", round(eta2, 3)), eta_word(eta2)))
    }
    if (!is.null(welch_a)) {
      add_test("Welch ANOVA(unequal variances)", unname(welch_a$statistic),
               welch_a$p.value, eta2,
               sprintf("Variance-robust cross-check: %s.", sig_phrase(welch_a$p.value)))
    }
    if (!is.null(kw)) {
      add_test("Kruskal-Wallis", unname(kw$statistic), kw$p.value, NA_real_,
               sprintf("Rank-based(no normality assumption): %s.", sig_phrase(kw$p.value)))
    }
    if (!is.na(shapiro_p)) {
      add_test("Shapiro-Wilk normality(residuals)", NA_real_, shapiro_p, NA_real_,
               if (shapiro_p < 0.01)
                 "Normality clearly violated — favor the rank-based result."
               else "No strong evidence against normality.")
    }

    tk <- tryCatch(stats::TukeyHSD(fit)$gf, error = function(e) NULL)
    if (!is.null(tk) && nrow(tk) > 0) {
      pairwise_df <- data.frame(
        comparison = rownames(tk),
        difference = round(as.numeric(tk[, "diff"]), 3),
        ci_low = round(as.numeric(tk[, "lwr"]), 3),
        ci_high = round(as.numeric(tk[, "upr"]), 3),
        adj_p = signif(as.numeric(tk[, "p adj"]), 3),
        stringsAsFactors = FALSE
      )
      pairwise_df$significant <- ifelse(is.na(pairwise_df$adj_p), "",
                                        ifelse(pairwise_df$adj_p < 0.05, "yes", "no"))
      ok <- !is.na(pairwise_df$difference)
      pairwise_df <- pairwise_df[ok, , drop = FALSE]
      pairwise_df <- pairwise_df[order(-abs(pairwise_df$difference)), , drop = FALSE]
      pairwise_df <- head(pairwise_df, 15)
      rownames(pairwise_df) <- NULL
    } else {
      pairwise_df <- data.frame(comparison = character(0), difference = numeric(0),
                                ci_low = numeric(0), ci_high = numeric(0),
                                adj_p = numeric(0), significant = character(0),
                                stringsAsFactors = FALSE)
    }

    effect_value <- eta2
    effect_text <- sprintf("eta-squared = %s(%s effect — share of variation in %s explained by %s)",
                           ifelse(is.na(eta2), "NA", round(eta2, 3)),
                           eta_word(eta2), outcome_h, group_h)
    if (nonparam_preferred && !is.null(kw)) {
      primary_test <- "Kruskal-Wallis"; primary_p <- kw$p.value
    } else if (!is.null(welch_a) && !is.na(welch_a$p.value)) {
      primary_test <- "Welch ANOVA"; primary_p <- welch_a$p.value
    } else if (!is.na(anova_p)) {
      primary_test <- "One-way ANOVA"; primary_p <- anova_p
    } else if (!is.null(kw)) {
      primary_test <- "Kruskal-Wallis"; primary_p <- kw$p.value
    } else {
      stop(sprintf("No comparison test could be computed for &#x27;%s' across the %s groups (the data may be essentially constant).",
                   outcome_h, group_h))
    }
  }

  test_results_df <- do.call(rbind, test_rows)
  rownames(test_results_df) <- NULL
  primary_sig <- !is.na(primary_p) && primary_p < 0.05

top_sig_pair <- NULL
  if (nrow(pairwise_df) > 0) {
    sig_rows <- pairwise_df[!is.na(pairwise_df$difference) &
                              pairwise_df$significant == "yes", , drop = FALSE]
    if (nrow(sig_rows) > 0) top_sig_pair <- sig_rows[1, ]  # already |diff|-sorted
  }
  n_sig_pairs <- sum(pairwise_df$significant == "yes", na.rm = TRUE)

set.seed(42)
  sidx <- if (final_rows > 2000) sample(final_rows, 2000) else seq_len(final_rows)
  boxplot_df <- data.frame(
    group_name = as.character(gf[sidx]),
    outcome_value = y[sidx],
    stringsAsFactors = FALSE
  )

  metrics <- list(
    `Observations`      = final_rows,
    `Groups Compared`   = k,
    `Primary Test`      = primary_test,
    `Primary p-value`   = signif(primary_p, 3),
    `Effect Size`       = if (is.na(effect_value)) NA_real_ else round(effect_value, 3),
    `Significant Pairs` = as.integer(n_sig_pairs)
  )

  json_output <- list(
    answer = paste0(
      "Compared ", outcome_h, " across ", k, " ", group_h, " groups(",
      format(final_rows, big.mark = ","), " rows): the ", primary_test,
      " finds the difference ", sig_phrase(primary_p), "; ", effect_text, ". ",
      if (!is.null(top_sig_pair)) paste0(
        "Largest significant gap: ", top_sig_pair$comparison, " (difference = ",
        top_sig_pair$difference, ", adjusted p = ", top_sig_pair$adj_p, "); ")
      else "No individual pair reaches significance; ",
      n_sig_pairs, " of ", nrow(pairwise_df), " pairwise comparison(s) significant at p<0.05."
    ),
    cards = lapply(
      c("tldr", "overview", "preprocessing", "group_boxplot",
        "group_summary", "test_results", "pairwise_differences"),
      function(cid) list(id = cid, metrics = metrics)
    )
  )

  list(
    initial_rows = initial_rows, final_rows = final_rows,
    rows_removed = rows_removed, n_na_outcome = n_na_outcome,
    outcome_h = outcome_h, group_h = group_h,
    k = k, group_levels = group_levels,
    dropped_groups_df = dropped_groups_df, lumped_levels = lumped_levels,
    group_summary_df = group_summary_df,
    test_results_df = test_results_df,
    pairwise_df = pairwise_df,
    boxplot_df = boxplot_df,
    shapiro_p = shapiro_p, out_frac = out_frac,
    nonparam_preferred = nonparam_preferred,
    primary_test = primary_test, primary_p = primary_p,
    primary_sig = primary_sig,
    effect_value = effect_value, effect_text = effect_text,
    top_sig_pair = top_sig_pair, n_sig_pairs = n_sig_pairs,
    metrics = metrics, json_output = json_output
  )
}