# Thesis-It statistics families, contract 1.
# Copyright (C) 2026 Thesis-It.
#
# This program is free software: you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by the Free
# Software Foundation, either version 3 of the License, or (at your option)
# any later version. It is distributed WITHOUT ANY WARRANTY. See
# <https://www.gnu.org/licenses/gpl-3.0.html>.
#
# Each family takes a data frame and column positions and returns a flat
# numeric vector, read by its fixed layout. They run inside webR on the
# student's device.

.thesisit <- new.env()
local({
  num <- function(d, k) suppressWarnings(as.numeric(d[[k]]))
  cat <- function(d, k) { x <- as.character(d[[k]]); x[!is.na(x) & trimws(x) == ""] <- NA; x }
  has <- function(x) if (length(x) == 0) NA_real_ else x

  assign("descriptives", function(d, cols) {
    out <- vapply(cols, function(k) {
      x <- num(d, k); missing <- sum(is.na(x)); x <- x[!is.na(x)]; n <- length(x)
      m <- if (n > 0) mean(x) else NA_real_
      s <- if (n > 1) sd(x) else NA_real_
      g1 <- if (n > 2 && isTRUE(s > 0)) n / ((n - 1) * (n - 2)) * sum(((x - m) / s)^3) else NA_real_
      g2 <- if (n > 3 && isTRUE(s > 0)) n * (n + 1) / ((n - 1) * (n - 2) * (n - 3)) * sum(((x - m) / s)^4) - 3 * (n - 1)^2 / ((n - 2) * (n - 3)) else NA_real_
      c(n, missing, m, s, if (n > 0) median(x) else NA_real_, if (n > 0) min(x) else NA_real_, if (n > 0) max(x) else NA_real_, g1, g2)
    }, numeric(9))
    as.vector(out)
  }, envir = .thesisit)

  assign("reliability", function(d, cols, level) {
    items <- vapply(cols, function(k) num(d, k), numeric(nrow(d)))
    items <- items[stats::complete.cases(items), , drop = FALSE]
    n <- nrow(items); k <- ncol(items)
    if (n < 2 || k < 2) return(c(k, n, rep(NA_real_, 4)))
    alpha <- k / (k - 1) * (1 - sum(apply(items, 2, stats::var)) / stats::var(rowSums(items)))
    a <- 1 - level
    lower <- 1 - (1 - alpha) * stats::qf(1 - a / 2, n - 1, (n - 1) * (k - 1))
    upper <- 1 - (1 - alpha) * stats::qf(a / 2, n - 1, (n - 1) * (k - 1))
    r <- stats::cor(items); rbar <- mean(r[upper.tri(r)])
    standardized <- k * rbar / (1 + (k - 1) * rbar)
    c(k, n, alpha, lower, upper, standardized)
  }, envir = .thesisit)

  assign("correlation", function(d, cols, method) {
    x <- vapply(cols, function(k) num(d, k), numeric(nrow(d)))
    k <- length(cols)
    r <- p <- n <- matrix(NA_real_, k, k)
    for (i in seq_len(k)) for (j in seq_len(k)) {
      ok <- !is.na(x[, i]) & !is.na(x[, j]); n[i, j] <- sum(ok)
      if (i != j && n[i, j] > 2 && stats::sd(x[ok, i]) > 0 && stats::sd(x[ok, j]) > 0) {
        test <- suppressWarnings(stats::cor.test(x[ok, i], x[ok, j], method = method, exact = FALSE))
        r[i, j] <- unname(test$estimate); p[i, j] <- test$p.value
      }
    }
    means <- vapply(seq_len(k), function(i) mean(x[, i], na.rm = TRUE), numeric(1))
    sds <- vapply(seq_len(k), function(i) stats::sd(x[, i], na.rm = TRUE), numeric(1))
    ns <- vapply(seq_len(k), function(i) sum(!is.na(x[, i])), numeric(1))
    listwise <- sum(stats::complete.cases(x))
    c(as.vector(r), as.vector(p), as.vector(n), means, sds, ns, listwise)
  }, envir = .thesisit)

  assign("ttest", function(d, outcome, group, level) {
    y <- num(d, outcome); g <- cat(d, group)
    ok <- !is.na(y) & !is.na(g); y <- y[ok]; g <- factor(g[ok])
    if (nlevels(g) != 2) stop("two groups")
    lv <- levels(g); a <- g == lv[1]; b <- g == lv[2]
    student <- stats::t.test(y[a], y[b], var.equal = TRUE, conf.level = level)
    welch <- stats::t.test(y[a], y[b], var.equal = FALSE, conf.level = level)
    deviation <- abs(y - stats::ave(y, g))
    levene <- stats::anova(stats::lm(deviation ~ g))
    n1 <- sum(a); n2 <- sum(b); m1 <- mean(y[a]); m2 <- mean(y[b]); s1 <- stats::sd(y[a]); s2 <- stats::sd(y[b])
    pooled <- sqrt(((n1 - 1) * s1^2 + (n2 - 1) * s2^2) / (n1 + n2 - 2))
    cohen <- (m1 - m2) / pooled
    hedges <- cohen * (1 - 3 / (4 * (n1 + n2) - 9))
    row <- function(t) c(unname(t$statistic), unname(t$parameter), t$p.value, m1 - m2, unname(t$stderr), t$conf.int[1], t$conf.int[2])
    c(levene[["F value"]][1], levene[["Pr(>F)"]][1], row(student), row(welch), cohen, hedges, n1, n2, m1, m2, s1, s2)
  }, envir = .thesisit)

  assign("chisquare", function(d, first, second) {
    x <- cat(d, first); y <- cat(d, second)
    ok <- !is.na(x) & !is.na(y)
    table <- base::table(factor(x[ok]), factor(y[ok]))
    if (min(dim(table)) < 2) stop("two categories each")
    test <- suppressWarnings(stats::chisq.test(table, correct = FALSE))
    expected <- test$expected
    n <- sum(table); k <- min(dim(table))
    value <- unname(test$statistic)
    v <- sqrt(value / (n * (k - 1)))
    square <- all(dim(table) == 2)
    phi <- if (square) (table[1, 1] * table[2, 2] - table[1, 2] * table[2, 1]) / sqrt(prod(rowSums(table)) * prod(colSums(table))) else NA_real_
    fisher <- if (square) stats::fisher.test(table)$p.value else NA_real_
    c(value, unname(test$parameter), test$p.value, sum(expected < 5), 100 * mean(expected < 5), min(expected), n, v, phi, fisher, nrow(table), ncol(table))
  }, envir = .thesisit)

  assign("regression", function(d, outcome, predictors, level) {
    frame <- data.frame(y = num(d, outcome))
    for (i in seq_along(predictors)) frame[[paste0("x", i)]] <- num(d, predictors[i])
    frame <- frame[stats::complete.cases(frame), , drop = FALSE]
    fit <- stats::lm(y ~ ., data = frame)
    if (any(is.na(stats::coef(fit)))) stop("a predictor is a copy of another")
    s <- summary(fit); ci <- stats::confint(fit, level = level); co <- s$coefficients
    sdy <- stats::sd(frame$y)
    beta <- c(NA_real_, vapply(seq_along(predictors), function(i) co[i + 1, 1] * stats::sd(frame[[paste0("x", i)]]) / sdy, numeric(1)))
    tolerance <- c(NA_real_, vapply(seq_along(predictors), function(i) {
      if (length(predictors) < 2) return(1)
      others <- setdiff(paste0("x", seq_along(predictors)), paste0("x", i))
      1 - summary(stats::lm(stats::reformulate(others, paste0("x", i)), data = frame))$r.squared
    }, numeric(1)))
    rows <- as.vector(t(cbind(co[, 1], co[, 2], beta, co[, 3], co[, 4], ci[, 1], ci[, 2], tolerance, 1 / tolerance)))
    f <- s$fstatistic
    residuals <- stats::residuals(fit)
    dw <- sum(diff(residuals)^2) / sum(residuals^2)
    model <- c(sqrt(s$r.squared), s$r.squared, s$adj.r.squared, s$sigma, unname(f[1]), unname(f[2]), unname(f[3]), stats::pf(f[1], f[2], f[3], lower.tail = FALSE), dw, nrow(frame))
    c(rows, model)
  }, envir = .thesisit)

  assign("anova", function(d, outcome, factor) {
    y <- num(d, outcome); g <- cat(d, factor)
    ok <- !is.na(y) & !is.na(g); y <- y[ok]; g <- base::factor(g[ok])
    if (nlevels(g) < 2) stop("two groups")
    tab <- summary(stats::aov(y ~ g))[[1]]
    ss <- tab[["Sum Sq"]]; df <- tab[["Df"]]; ms <- tab[["Mean Sq"]]
    total <- sum(ss)
    eta <- ss[1] / total
    omega <- (ss[1] - df[1] * ms[2]) / (total + ms[2])
    c(ss[1], df[1], ms[1], tab[["F value"]][1], tab[["Pr(>F)"]][1], eta, omega, ss[2], df[2], ms[2], total, sum(df), length(y))
  }, envir = .thesisit)
})
invisible(NULL)