Normality & Outlier Screening — Is Your Data Ready?

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))

u <- (n * G^2) / ((n - 1)^2)
  if (u >= 1) {
    p <- 0
  } else {
    t2 <- u * (n - 2) / (1 - u)
    tval <- sqrt(t2)
    p <- n * 2 * pt(-tval, df = n - 2)
  }
  p <- min(1, max(0, p))

tcrit <- qt(alpha / (2 * n), df = n - 2, lower.tail = FALSE)
  crit <- ((n - 1) / sqrt(n)) * sqrt(tcrit^2 / (n - 2 + tcrit^2))

  list(G = G, crit = as.numeric(crit), p = p,
       value = as.numeric(x[i]), flagged = isTRUE(G > crit))
}

compute_shared <- function(df, params, col_map = list()) {
  # === SHARED EXPORTS ===
  #   initial_rows/final_rows/rows_removed  $ row accounting
  #   col_names          $ named character — semantic -> humanized names
  #   used_cols          $ character — semantic names screened (incl. constant)
  #   dropped_cols       $ character — excluded (non-numeric / empty)
  #   stats              $ per-column list of computed stats (aligned to used_cols)
  #   normality_df       $ column, n, skew, kurtosis, shapiro_p, verdict
  #   outliers_df        $ column, mild_outliers, extreme_outliers,
  #                        most_extreme_value, grubbs_p
  #   skew_by_column_df  $ column, abs_skew  (bar; constant/NA -> 0)
  #   n_used/n_constant  $ counts
  #   n_normal/n_nonnormal/n_nonparam  $ verdict + recommendation counts
  #   total_mild/total_extreme         $ outlier totals across columns
  #   n_sampled          $ how many columns were Shapiro-sampled at 5,000
  #   most_name/most_p   $ single most non-normal column + its Shapiro p
  #   metrics / json_output
  # === /SHARED EXPORTS ===

initial_rows <- nrow(df)
  col_cols <- grep("^column_[0-9]+$", names(df), value = TRUE)
  col_cols <- col_cols[order(as.integer(sub("^column_", "", col_cols)))]
  if (length(col_cols) < 1) {
    stop("column_mapping must map at least one numeric column(column_1).")
  }
  col_names <- setNames(humanize_semantic(col_cols, col_map), col_cols)

dropped_cols <- character(0)
  for (cc in col_cols) {
    v <- df[[cc]]
    if (!is.numeric(v)) {
      conv <- suppressWarnings(as.numeric(as.character(v)))
      n_orig <- sum(!is.na(v) & as.character(v) != "")
      if (n_orig > 0 && sum(!is.na(conv)) >= 0.95 * n_orig) {
        df[[cc]] <- conv
      } else {
        dropped_cols <- c(dropped_cols, cc); next
      }
    }
    v <- suppressWarnings(as.numeric(df[[cc]]))
    if (all(is.na(v))) { dropped_cols <- c(dropped_cols, cc); next }
    df[[cc]] <- v
  }
  used_cols <- setdiff(col_cols, dropped_cols)
  if (length(used_cols) < 1) {
    stop(paste0("Screening needs at least one usable numeric column; none of ",
                "the mapped columns(",
                paste(unname(col_names[col_cols]), collapse = ", "),
                ") were numeric."))
  }
  final_rows <- initial_rows
  rows_removed <- 0L
  if (final_rows < 10) stop(sprintf("Only %d rows — screening needs at least 10.", final_rows))

stats <- list()
  for (cc in used_cols) {
    hn <- col_names[[cc]]
    x_all <- suppressWarnings(as.numeric(df[[cc]]))
    x <- x_all[is.finite(x_all)]
    n <- length(x)
    n_missing <- sum(!is.finite(x_all))

    mean_v <- if (n > 0) mean(x) else NA_real_
    median_v <- if (n > 0) median(x) else NA_real_
    sd_v <- if (n > 1) sd(x) else NA_real_
    constant <- (n < 2) || !is.finite(sd_v) || sd_v == 0

    sk <- if (constant) NA_real_ else skewness(x)
    ku <- if (constant) NA_real_ else excess_kurtosis(x)

shapiro_p <- NA_real_
    sampled <- FALSE
    shap_n <- n
    if (!constant && n >= 3) {
      xs <- x
      if (n > 5000) {
        set.seed(42)
        xs <- x[sample(n, 5000)]
        sampled <- TRUE
        shap_n <- 5000L
      }
      sp <- tryCatch(shapiro.test(xs)$p.value, error = function(e) NA_real_)
      shapiro_p <- as.numeric(sp)
    }

mild <- 0L; extreme <- 0L
    if (n > 0 && !constant) {
      q <- quantile(x, c(0.25, 0.75), names = FALSE, type = 7)
      iqr <- q[2] - q[1]
      if (is.finite(iqr) && iqr > 0) {
        mild <- as.integer(sum(x < q[1] - 1.5 * iqr | x > q[2] + 1.5 * iqr))
        extreme <- as.integer(sum(x < q[1] - 3 * iqr | x > q[2] + 3 * iqr))
      }
    }

gr <- grubbs_one(x)

rz_flagged <- 0L
    if (n > 0 && !constant) {
      med <- median(x)
      mad_s <- 1.4826 * median(abs(x - med))
      if (!is.finite(mad_s) || mad_s == 0) mad_s <- sd_v
      if (is.finite(mad_s) && mad_s > 0) {
        rz_flagged <- as.integer(sum(abs((x - med) / mad_s) > 3.5))
      }
    }

    verdict <- if (constant) {
      "constant(no variation)"
    } else if (!is.na(sk) && sk > 1) {
      "right-skewed"
    } else if (!is.na(sk) && sk < -1) {
      "left-skewed"
    } else if (!is.na(ku) && ku > 2) {
      "heavy-tailed"
    } else if (!is.na(shapiro_p) && shapiro_p >= 0.05) {
      "approximately normal"
    } else if (!is.na(shapiro_p) && shapiro_p < 0.05) {
      "mildly non-normal"
    } else {
      "approximately normal"
    }

    nonparam <- (!is.na(shapiro_p) && shapiro_p < 0.05) ||
                (!is.na(sk) && abs(sk) > 1)
    recommendation <- if (constant) {
      "constant — excluded from testing"
    } else if (nonparam) {
      "consider a nonparametric method or transformation"
    } else {
      "parametric methods are appropriate"
    }

    stats[[cc]] <- list(
      name = hn, n = n, n_missing = n_missing,
      mean = mean_v, median = median_v, sd = sd_v, constant = constant,
      skew = sk, kurtosis = ku, abs_skew = if (is.na(sk)) 0 else abs(sk),
      shapiro_p = shapiro_p, sampled = sampled, shap_n = shap_n,
      mild = mild, extreme = extreme,
      grubbs_G = gr$G, grubbs_crit = gr$crit, grubbs_p = gr$p,
      grubbs_value = gr$value, grubbs_flagged = gr$flagged,
      most_extreme_value = gr$value, robust_z_flagged = rz_flagged,
      verdict = verdict, recommendation = recommendation, nonparam = nonparam
    )
  }

constant_vec <- vapply(used_cols, function(cc) isTRUE(stats[[cc]]$constant), logical(1))
  verdict_vec <- vapply(used_cols, function(cc) stats[[cc]]$verdict, character(1))
  shapiro_vec <- vapply(used_cols, function(cc) stats[[cc]]$shapiro_p, numeric(1))
  absskew_vec <- vapply(used_cols, function(cc) stats[[cc]]$abs_skew, numeric(1))
  nonparam_vec <- vapply(used_cols, function(cc) isTRUE(stats[[cc]]$nonparam), logical(1))

  n_used <- length(used_cols)
  n_constant <- sum(constant_vec)
  n_testable <- n_used - n_constant
  n_normal <- sum(verdict_vec == "approximately normal")
  n_nonnormal <- max(0L, n_testable - n_normal)
  n_nonparam <- sum(nonparam_vec)
  total_mild <- sum(vapply(used_cols, function(cc) stats[[cc]]$mild, integer(1)))
  total_extreme <- sum(vapply(used_cols, function(cc) stats[[cc]]$extreme, integer(1)))
  n_sampled <- sum(vapply(used_cols, function(cc) isTRUE(stats[[cc]]$sampled), logical(1)))

testable_idx <- which(!constant_vec)
  most_name <- "none"; most_p <- NA_real_
  if (length(testable_idx) > 0) {
    sp_ok <- testable_idx[is.finite(shapiro_vec[testable_idx])]
    if (length(sp_ok) > 0) {
      mi <- sp_ok[which.min(shapiro_vec[sp_ok])]
    } else {
      as_ok <- testable_idx[is.finite(absskew_vec[testable_idx])]
      mi <- if (length(as_ok) > 0) as_ok[which.max(absskew_vec[as_ok])] else testable_idx[1]
    }
    most_name <- hn_vec[mi]
    most_p <- shapiro_vec[mi]
  }

  metrics <- list(
    `Columns Screened`         = as.integer(n_used),
    `Approximately Normal`     = as.integer(n_normal),
    `Non-Normal or Skewed`     = as.integer(n_nonnormal),
    `Extreme Outliers Flagged` = as.integer(total_extreme),
    `Nonparametric Recommended`= as.integer(n_nonparam),
    `Most Non-Normal Column`   = most_name
  )

  json_output <- list(
    answer = paste0(
      "Screened ", n_used, " numeric column", .s(n_used), " across ",
      format(final_rows, big.mark = ","), " row", .s(final_rows), ": ",
      n_normal, " look", if (n_normal == 1) "s" else "",
      " approximately normal and ", n_nonnormal, " ",
      if (n_nonnormal == 1) "is" else "are", " skewed or non-normal. ",
      total_extreme, " extreme outlier", .s(total_extreme),
      " were found across the columns. ",
      if (most_name != "none") paste0(
        "The single most non-normal column is ", most_name,
        " (Shapiro-Wilk ", fmt_pp(most_p), "). ") else "",
      n_nonparam, " column", .s(n_nonparam), " ",
      if (n_nonparam == 1) "is" else "are",
      " better served by a nonparametric method or a transformation."
    ),
    cards = lapply(
      c("tldr", "overview", "preprocessing", "normality_table",
        "outlier_summary", "recommendation"),
      function(cid) list(id = cid, metrics = metrics)
    )
  )

  list(
    initial_rows = initial_rows, final_rows = final_rows,
    rows_removed = rows_removed,
    col_names = col_names, used_cols = used_cols, dropped_cols = dropped_cols,
    stats = stats,
    normality_df = normality_df, outliers_df = outliers_df,
    skew_by_column_df = skew_by_column_df,
    n_used = n_used, n_constant = n_constant, n_testable = n_testable,
    n_normal = n_normal, n_nonnormal = n_nonnormal, n_nonparam = n_nonparam,
    total_mild = total_mild, total_extreme = total_extreme,
    n_sampled = n_sampled, most_name = most_name, most_p = most_p,
    metrics = metrics, json_output = json_output
  )
}