# PivoGraph statistical validation: reference values computed with R.
# https://pivograph.com/validation/
#
# Run from the folder that holds this script (writes r-reference.json next to it):
#   Rscript r-reference.R
# Needs R >= 4.4 and the CRAN packages loaded below; MASS ships with R.
#
# Data are R's built-in data sets (or documented subsets), so anyone with R can
# reproduce every number. Conventions follow PivoGraph where R allows a choice:
# group differences are first - second (t-tests: B - A), pairs are listed as
# (1,2), (1,3), ..., (2,3), ... in the order of the factor levels.

suppressPackageStartupMessages({
  library(jsonlite)
  library(multcomp)
  library(PMCMRplus)
  library(coin)
  library(DescTools)
  library(effectsize)
  library(car)
  library(onewaytests)
  library(MASS)
  library(survival)
})

script_dir <- local({
  arg <- grep("^--file=", commandArgs(FALSE), value = TRUE)
  if (length(arg) == 1) dirname(normalizePath(sub("^--file=", "", arg))) else getwd()
})

pairs_of <- function(k) {
  out <- list()
  for (i in 1:(k - 1)) for (j in (i + 1):k) out[[length(out) + 1]] <- c(i, j)
  out
}

as_groups <- function(y, g) {
  g <- factor(g)
  list(names = levels(g), values = lapply(levels(g), function(l) unname(y[g == l])))
}

# Lower-triangle matrices from PMCMRplus: entry for levels (i, j), i < j.
pm_entry <- function(m, lv, i, j) unname(m[lv[j], lv[i]])

# ── Descriptive statistics ──────────────────────────────────────────────────

describe <- function(id, x) {
  n <- length(x)
  s <- sd(x)
  q <- unname(quantile(x, c(0.05, 0.10, 0.25, 0.75, 0.90, 0.95), type = 6))
  mci <- MedianCI(x, conf.level = 0.95, method = "exact")
  positive <- all(x > 0)
  list(
    id = id,
    values = x,
    n = n,
    sum = sum(x),
    mean = mean(x),
    sd = s,
    sem = s / sqrt(n),
    meanCi = unname(t.test(x)$conf.int[1:2]),
    min = min(x),
    q1 = q[3],
    median = median(x),
    q3 = q[4],
    max = max(x),
    percentiles = list(p = c(5, 10, 90, 95), value = q[c(1, 2, 5, 6)]),
    medianCi = if (is.finite(mci[["lwr.ci"]])) unname(mci[c("lwr.ci", "upr.ci")]) else NA,
    medianCiLevel = attr(mci, "conf.level"),
    cv = s / mean(x),
    geoMean = if (positive) exp(mean(log(x))) else NA,
    geoSdFactor = if (positive) exp(sd(log(x))) else NA,
    geoMeanCi = if (positive) exp(unname(t.test(log(x))$conf.int[1:2])) else NA,
    skewness = Skew(x, method = 2),
    kurtosis = Kurt(x, method = 2)
  )
}

# ── Multiple-comparison helpers (formulas, applied to R's raw P values) ───

sidak <- function(p) 1 - (1 - p)^length(p)

holm_sidak <- function(p) {
  m <- length(p)
  o <- order(p)
  adj <- 1 - (1 - p[o])^(m - seq_len(m) + 1)
  adj <- cummax(adj)
  out <- numeric(m)
  out[o] <- pmin(1, adj)
  out
}

# Benjamini, Krieger & Yekutieli (2006) two-stage step-up, as defined in the paper.
bky_discoveries <- function(p, q) {
  m <- length(p)
  q1 <- q / (1 + q)
  r1 <- sum(p.adjust(p, "BH") <= q1)
  if (r1 == 0) return(rep(FALSE, m))
  if (r1 == m) return(rep(TRUE, m))
  m0 <- m - r1
  p.adjust(p, "BH") <= q1 * m / m0
}

# ── Dunnett and studentized maximum modulus distributions ──────────────────
#
# multcomp's P values come from randomized integration (error about 1e-4 beyond two
# comparisons) and qmvt's critical value from a loose root search (PlantGrowth:
# 2.333538 where the exact 95% point is 2.333412). The reference therefore integrates
# the definitions directly with R's QUADPACK routines; multcomp and PMCMRplus are kept
# as secondary checks.

# P(max |T_i| <= d) for comparisons with a control: rho_ij = lambda_i lambda_j,
# lambda_i = sqrt(n_i / (n_i + n_0)), S = sqrt(chi2_df / df).
dunnett_cdf <- function(d, lambda, df) {
  sigma <- sqrt(1 - lambda^2)
  inner <- function(s) vapply(s, function(si) integrate(function(z) {
    out <- dnorm(z)
    for (i in seq_along(lambda)) {
      out <- out * (pnorm((d * si - lambda[i] * z) / sigma[i]) -
                      pnorm((-d * si - lambda[i] * z) / sigma[i]))
    }
    out
  }, -Inf, Inf, rel.tol = 1e-12)$value, numeric(1))
  integrate(function(s) inner(s) * dchisq(s^2 * df, df) * 2 * s * df, 0, Inf,
            rel.tol = 1e-11)$value
}

dunnett_critical <- function(level, lambda, df) {
  uniroot(function(d) dunnett_cdf(d, lambda, df) - level, c(1, 10), tol = 1e-12)$root
}

# P(max |Z_i| / S <= c) for m independent comparisons sharing S (Dunnett T3).
smm_cdf <- function(c, m, df) {
  integrate(function(s) (2 * pnorm(c * s) - 1)^m * dchisq(s^2 * df, df) * 2 * s * df,
            0, Inf, rel.tol = 1e-13)$value
}

# ── Ordinary one-way ANOVA ──────────────────────────────────────────────────

ordinary <- function(id, y, g) {
  d <- data.frame(y = y, g = factor(g))
  lv <- levels(d$g)
  k <- length(lv)
  fit <- aov(y ~ g, data = d)
  tab <- summary(fit)[[1]]
  dfe <- tab[2, "Df"]
  mse <- tab[2, "Mean Sq"]
  means <- tapply(d$y, d$g, mean)
  ns <- tapply(d$y, d$g, length)
  prs <- pairs_of(k)

  lsd <- lapply(prs, function(pr) {
    i <- pr[1]; j <- pr[2]
    diff <- unname(means[i] - means[j])
    se <- unname(sqrt(mse * (1 / ns[i] + 1 / ns[j])))
    half <- qt(0.975, dfe) * se
    list(first = i - 1, second = j - 1, diff = diff, low = diff - half, high = diff + half,
         t = diff / se, p = 2 * pt(-abs(diff / se), dfe))
  })
  p_raw <- sapply(lsd, `[[`, "p")

  tk <- TukeyHSD(fit)$g
  tukey <- lapply(prs, function(pr) {
    row <- tk[paste(lv[pr[2]], lv[pr[1]], sep = "-"), ]
    list(first = pr[1] - 1, second = pr[2] - 1, diff = -row[["diff"]],
         low = -row[["upr"]], high = -row[["lwr"]], p = row[["p adj"]])
  })

  set.seed(20261030)
  ds <- summary(glht(fit, linfct = mcp(g = "Dunnett")))
  lambda <- unname(sqrt(ns[-1] / (ns[-1] + ns[1])))
  calpha <- dunnett_critical(0.95, lambda, dfe)
  dunnett <- lapply(seq_len(k - 1), function(r) {
    j <- r + 1
    diff <- unname(means[1] - means[j])
    se <- unname(sqrt(mse * (1 / ns[1] + 1 / ns[j])))
    list(first = 0, second = r, diff = diff, low = diff - calpha * se, high = diff + calpha * se,
         p = 1 - dunnett_cdf(abs(diff / se), lambda, dfe),
         pMultcomp = unname(ds$test$pvalues[r]))
  })

  lev <- leveneTest(y ~ g, data = d, center = median)
  bt <- bartlett.test(y ~ g, data = d)
  list(
    id = id,
    groups = as_groups(d$y, d$g),
    f = tab[1, "F value"], dfn = tab[1, "Df"], dfd = dfe, p = tab[1, "Pr(>F)"],
    ssTreatment = tab[1, "Sum Sq"], ssResidual = tab[2, "Sum Sq"],
    rSquared = tab[1, "Sum Sq"] / (tab[1, "Sum Sq"] + tab[2, "Sum Sq"]),
    brownForsythe = list(f = lev[1, "F value"], dfn = lev[1, "Df"], dfd = lev[2, "Df"],
                         p = lev[1, "Pr(>F)"]),
    bartlett = list(statistic = unname(bt$statistic), df = unname(bt$parameter), p = bt$p.value),
    lsd = lsd,
    tukey = tukey,
    dunnett = dunnett,
    dunnettCritical = calpha,
    multcompError = attr(ds$test$pvalues, "error"),
    bonferroni = p.adjust(p_raw, "bonferroni"),
    sidak = sidak(p_raw),
    holmSidak = holm_sidak(p_raw),
    bh = p.adjust(p_raw, "BH"),
    bkyDiscoveries05 = bky_discoveries(p_raw, 0.05)
  )
}

# ── Welch and Brown-Forsythe ANOVA ──────────────────────────────────────────

welch <- function(id, y, g) {
  d <- data.frame(y = y, g = factor(g))
  lv <- levels(d$g)
  ow <- oneway.test(y ~ g, data = d, var.equal = FALSE)
  bf <- bf.test(y ~ g, data = d, verbose = FALSE)
  gh <- gamesHowellTest(y ~ g, data = d)
  t3 <- dunnettT3Test(y ~ g, data = d)
  prs <- pairs_of(length(lv))
  rows <- lapply(prs, function(pr) {
    a <- d$y[d$g == lv[pr[1]]]
    b <- d$y[d$g == lv[pr[2]]]
    tt <- t.test(a, b, var.equal = FALSE)
    list(first = pr[1] - 1, second = pr[2] - 1, t = unname(tt$statistic),
         df = unname(tt$parameter), p = tt$p.value,
         gamesHowellP = pm_entry(gh$p.value, lv, pr[1], pr[2]),
         dunnettT3P = 1 - smm_cdf(abs(unname(tt$statistic)), length(prs), unname(tt$parameter)),
         dunnettT3PMCMRplus = pm_entry(t3$p.value, lv, pr[1], pr[2]))
  })
  list(
    id = id,
    groups = as_groups(d$y, d$g),
    f = unname(ow$statistic), dfn = unname(ow$parameter[1]), dfd = unname(ow$parameter[2]),
    p = ow$p.value,
    brownForsythe = list(f = bf$statistic, dfn = bf$parameter[1], dfd = bf$parameter[2],
                         p = bf$p.value),
    comparisons = rows
  )
}

# ── Kruskal-Wallis and Dunn ─────────────────────────────────────────────────

# Exact permutation P conditional on the observed (mid)ranks: share of all ways to
# split the ranks into groups of the observed sizes with sum(R^2 / n) >= observed.
exact_kw <- function(y, g) {
  g <- factor(g)
  r <- rank(y)
  sizes <- as.vector(table(g))
  k <- length(sizes)
  observed <- sum(tapply(r, g, sum)^2 / sizes)
  tol <- 1e-9 * max(1, abs(observed))
  hits <- 0
  leaves <- 0
  walk <- function(remaining, gi, acc) {
    if (gi == k - 1) {
      cm <- combn(length(remaining), sizes[gi])
      s <- colSums(matrix(remaining[cm], nrow = sizes[gi]))
      last <- sum(remaining) - s
      v <- acc + s^2 / sizes[gi] + last^2 / sizes[k]
      hits <<- hits + sum(v >= observed - tol)
      leaves <<- leaves + length(v)
      return(invisible())
    }
    cm <- combn(length(remaining), sizes[gi])
    for (c in seq_len(ncol(cm))) {
      pick <- cm[, c]
      walk(remaining[-pick], gi + 1, acc + sum(remaining[pick])^2 / sizes[gi])
    }
  }
  walk(r, 1, 0)
  hits / leaves
}

kruskal <- function(id, y, g, exact = FALSE) {
  d <- data.frame(y = y, g = factor(g))
  lv <- levels(d$g)
  kt <- kruskal.test(y ~ g, data = d)
  dn <- suppressWarnings(kwAllPairsDunnTest(y ~ g, data = d, p.adjust.method = "none"))
  dnb <- suppressWarnings(kwAllPairsDunnTest(y ~ g, data = d, p.adjust.method = "bonferroni"))
  rows <- lapply(pairs_of(length(lv)), function(pr) {
    list(first = pr[1] - 1, second = pr[2] - 1,
         z = abs(pm_entry(dn$statistic, lv, pr[1], pr[2])),
         p = pm_entry(dn$p.value, lv, pr[1], pr[2]),
         pBonferroni = pm_entry(dnb$p.value, lv, pr[1], pr[2]))
  })
  list(
    id = id,
    groups = as_groups(d$y, d$g),
    h = unname(kt$statistic), df = unname(kt$parameter), pChiSquare = kt$p.value,
    pExact = if (exact) exact_kw(d$y, d$g) else NA,
    dunn = rows
  )
}

# ── Friedman, repeated-measures ANOVA ───────────────────────────────────────

permutations <- function(k) {
  if (k == 1) return(matrix(1L, 1, 1))
  sub <- permutations(k - 1)
  do.call(rbind, lapply(seq_len(k), function(first) {
    rest <- setdiff(seq_len(k), first)
    cbind(first, matrix(rest[sub], nrow = nrow(sub)))
  }))
}

# Exact P conditional on the within-row ranks (ties kept): every distinct
# arrangement of each row's ranks is equally likely.
exact_friedman <- function(m) {
  k <- ncol(m)
  ranks <- t(apply(m, 1, rank))
  observed <- sum(colSums(ranks)^2)
  tol <- 1e-9 * max(1, observed)
  idx <- permutations(k)
  sums <- matrix(0, 1, k)
  for (r in seq_len(nrow(ranks))) {
    arr <- unique(matrix(ranks[r, idx], ncol = k))
    sums <- sums[rep(seq_len(nrow(sums)), each = nrow(arr)), , drop = FALSE] +
      arr[rep(seq_len(nrow(arr)), times = nrow(sums)), , drop = FALSE]
  }
  mean(rowSums(sums^2) >= observed - tol)
}

friedman <- function(id, m, exact = FALSE) {
  lv <- colnames(m)
  ft <- friedman.test(m)
  sg <- frdAllPairsSiegelTest(m, p.adjust.method = "none")
  rows <- lapply(pairs_of(ncol(m)), function(pr) {
    list(first = pr[1] - 1, second = pr[2] - 1,
         z = abs(pm_entry(sg$statistic, lv, pr[1], pr[2])),
         p = pm_entry(sg$p.value, lv, pr[1], pr[2]))
  })
  list(
    id = id,
    names = lv,
    rows = unname(lapply(seq_len(nrow(m)), function(r) unname(m[r, ]))),
    statistic = unname(ft$statistic), df = unname(ft$parameter), pChiSquare = ft$p.value,
    pExact = if (exact) exact_friedman(m) else NA,
    dunn = rows
  )
}

repeated <- function(id, m) {
  lv <- colnames(m)
  k <- ncol(m)
  n <- nrow(m)
  fit <- lm(m ~ 1)
  av <- Anova(fit, idata = data.frame(cond = factor(lv, levels = lv)), idesign = ~cond, type = "III")
  sm <- summary(av, multivariate = FALSE)
  uni <- sm$univariate.tests["cond", ]
  gg <- sm$pval.adjustments["cond", ]
  # Pairwise paired comparisons (paired SE per pair, df = n - 1).
  rows <- lapply(pairs_of(k), function(pr) {
    diffs <- m[, pr[1]] - m[, pr[2]]
    t <- mean(diffs) / (sd(diffs) / sqrt(n))
    list(first = pr[1] - 1, second = pr[2] - 1, diff = mean(diffs), t = t,
         p = 2 * pt(-abs(t), n - 1),
         pTukey = ptukey(sqrt(2) * abs(t), k, n - 1, lower.tail = FALSE))
  })
  list(
    id = id,
    names = lv,
    rows = unname(lapply(seq_len(n), function(r) unname(m[r, ]))),
    f = uni[["F value"]], dfn = uni[["num Df"]], dfd = uni[["den Df"]], p = uni[["Pr(>F)"]],
    epsilon = gg[["GG eps"]], pGeisserGreenhouse = gg[["Pr(>F[GG])"]],
    comparisons = rows
  )
}

# ── Two-sample tests (differences are B - A) ────────────────────────────────

ttest_unpaired <- function(id, a, b) {
  w <- t.test(b, a, var.equal = FALSE)
  s <- t.test(b, a, var.equal = TRUE)
  big <- if (var(a) >= var(b)) a else b
  small <- if (var(a) >= var(b)) b else a
  vt <- var.test(big, small)
  tt <- function(x) list(t = abs(unname(x$statistic)), df = unname(x$parameter), p = x$p.value,
                         diff = unname(x$estimate[1] - x$estimate[2]),
                         low = x$conf.int[1], high = x$conf.int[2])
  list(
    id = id, a = a, b = b,
    welch = tt(w), student = tt(s),
    rSquaredWelch = unname(w$statistic^2 / (w$statistic^2 + w$parameter)),
    varianceF = list(f = unname(vt$statistic), dfn = unname(vt$parameter[1]),
                     dfd = unname(vt$parameter[2]), p = vt$p.value),
    cohensD = cohens_d(b, a, pooled_sd = TRUE)$Cohens_d,
    hedgesG = hedges_g(b, a, pooled_sd = TRUE)$Hedges_g
  )
}

ttest_paired <- function(id, a, b) {
  p <- t.test(b, a, paired = TRUE)
  ct <- cor.test(a, b, alternative = "greater")
  list(
    id = id, a = a, b = b,
    t = abs(unname(p$statistic)), df = unname(p$parameter), p = p$p.value,
    diff = unname(p$estimate), low = p$conf.int[1], high = p$conf.int[2],
    rSquared = unname(p$statistic^2 / (p$statistic^2 + p$parameter)),
    dz = repeated_measures_d(b, a, method = "z", adjust = FALSE)[[1]],
    pairingR = unname(ct$estimate), pairingP = ct$p.value
  )
}

# R's exact interval for wilcox.test(conf.int = TRUE): achieved level of
# [d(qu), d(n1 n2 - qu + 1)] with qu = qwilcox(alpha / 2, n1, n2).
wilcox_achieved_level <- function(n1, n2, level = 0.95) {
  qu <- qwilcox((1 - level) / 2, n1, n2)
  if (qu == 0) qu <- 1
  1 - 2 * pwilcox(trunc(qu) - 1, n1, n2)
}

mann_whitney <- function(id, a, b, kind) {
  if (kind == "exact") {
    w <- wilcox.test(b, a, exact = TRUE, conf.int = TRUE)
    p_exact <- w$p.value
    level <- wilcox_achieved_level(length(b), length(a))
  } else {
    w <- wilcox.test(b, a, exact = FALSE, correct = TRUE, conf.int = TRUE)
    p_exact <- if (kind == "ties") {
      d <- data.frame(y = c(a, b), g = factor(rep(c("A", "B"), c(length(a), length(b)))))
      pvalue(wilcox_test(y ~ g, data = d, distribution = "exact"))
    } else NA
    level <- NA
  }
  wb <- unname(w$statistic) # U of B
  list(
    id = id, a = a, b = b, kind = kind,
    u = min(wb, length(a) * length(b) - wb),
    p = if (kind == "approximate") w$p.value else p_exact,
    pApproximate = if (kind == "exact") NA else w$p.value,
    hodgesLehmann = unname(w$estimate),
    ci = unname(w$conf.int[1:2]),
    ciLevel = level
  )
}

# Differences are rounded to 12 significant digits before ranking, as in PivoGraph:
# floating-point subtraction otherwise splits ties (anorexia: 64 distinct |d| instead of 61).
wilcoxon <- function(id, a, b, kind) {
  d <- signif(b - a, 12)
  nz <- d[d != 0]
  ranks <- rank(abs(nz))
  sum_pos <- sum(ranks[nz > 0])
  sum_neg <- sum(ranks[nz < 0])
  p <- switch(kind,
    exact = wilcox.test(d, exact = TRUE)$p.value,
    ties = pvalue(wilcoxsign_test(d ~ rep(0, length(d)), zero.method = "Wilcoxon",
                                  distribution = "exact")),
    approximate = wilcox.test(d, exact = FALSE, correct = TRUE)$p.value
  )
  list(
    id = id, a = a, b = b, kind = kind,
    pairs = length(nz), zeros = sum(d == 0),
    sumPositive = sum_pos, sumNegative = sum_neg,
    p = unname(p),
    medianDifference = median(d),
    spearman = cor(a, b, method = "spearman")
  )
}

# ── XY: area under the curve, correlation ───────────────────────────────────

trapezoid <- function(id, x, y) {
  o <- order(x)
  x <- x[o]; y <- y[o]
  list(id = id, x = x, y = y, area = sum(diff(x) * (head(y, -1) + tail(y, -1)) / 2))
}

correlation <- function(id, x, y) {
  pe <- cor.test(x, y, method = "pearson")
  sp <- suppressWarnings(cor.test(x, y, method = "spearman", exact = FALSE))
  list(id = id, x = x, y = y,
       pearson = list(r = unname(pe$estimate), p = pe$p.value, ci = unname(pe$conf.int[1:2])),
       spearman = list(r = unname(sp$estimate), p = sp$p.value))
}

# ── Survival: Kaplan–Meier and curve comparison ─────────────────────────────

# Gehan–Breslow–Wilcoxon from its definition: the log-rank score and its
# hypergeometric covariance, each event time weighted by the number at risk.
gehan_breslow_wilcoxon <- function(time, status, group) {
  lv <- levels(group); k <- length(lv)
  u <- numeric(k); v <- matrix(0, k, k)
  for (t in sort(unique(time[status == 1]))) {
    nk <- sapply(lv, function(l) sum(time[group == l] >= t))
    dk <- sapply(lv, function(l) sum(time[group == l] == t & status[group == l] == 1))
    n <- sum(nk); d <- sum(dk); p <- nk / n
    u <- u + n * (dk - d * p)
    if (n > 1) v <- v + n^2 * d * (n - d) / (n - 1) * (diag(p) - outer(p, p))
  }
  m <- k - 1
  chisq <- as.numeric(t(u[1:m]) %*% solve(v[1:m, 1:m, drop = FALSE]) %*% u[1:m])
  list(chiSquare = chisq, df = m, p = pchisq(chisq, m, lower.tail = FALSE))
}

# status: 1 = event, 0 = censored. Groups in factor-level order; ratios are first / second.
survival_curves <- function(id, time, status, group, km = FALSE) {
  group <- factor(group); lv <- levels(group); k <- length(lv)
  lr <- survdiff(Surv(time, status) ~ group)
  fit <- survfit(Surv(time, status) ~ group, conf.type = "log-log")
  q <- quantile(fit, 0.5)
  med <- as.numeric(q$quantile)
  out <- list(
    id = id,
    groups = lapply(lv, function(l) list(name = l, time = time[group == l],
                                         status = status[group == l])),
    n = as.numeric(lr$n), observed = as.numeric(lr$obs), expected = as.numeric(lr$exp),
    logRank = list(chiSquare = lr$chisq, df = k - 1,
                   p = pchisq(lr$chisq, k - 1, lower.tail = FALSE)),
    gehanBreslowWilcoxon = gehan_breslow_wilcoxon(time, status, group),
    # Median survival and its 95% CI from the log-log pointwise limits (NA = not reached).
    median = med, medianLower = as.numeric(q$lower), medianUpper = as.numeric(q$upper))
  if (k >= 3) {
    # Log-rank test for trend, groups scored 1, 2, ..., k.
    w <- seq_len(k); u <- lr$obs - lr$exp
    chi <- sum(w * u)^2 / as.numeric(t(w) %*% lr$var %*% w)
    out$trend <- list(chiSquare = chi, df = 1, p = pchisq(chi, 1, lower.tail = FALSE))
  }
  if (k == 2) {
    z <- qnorm(0.975); o <- lr$obs; e <- lr$exp; v <- lr$var[1, 1]
    mh <- (o[1] - e[1]) / v
    lrr <- (o[1] / e[1]) / (o[2] / e[2])
    out$hazardRatio <- list(
      mantelHaenszel = list(ratio = exp(mh), ci = exp(mh + c(-1, 1) * z / sqrt(v))),
      logRank = list(ratio = lrr, ci = exp(log(lrr) + c(-1, 1) * z * sqrt(1 / e[1] + 1 / e[2]))))
    if (!anyNA(med)) {
      # Ratio of medians; 95% CI assuming exponential survival in both groups.
      r <- med[1] / med[2]
      out$medianRatio <- list(ratio = r, ci = r * exp(c(-1, 1) * z * sqrt(1 / o[1] + 1 / o[2])))
    }
  }
  if (km) {
    s <- summary(fit)
    strata <- as.integer(s$strata)
    out$km <- lapply(seq_len(k), function(i) list(
      time = s$time[strata == i], survival = s$surv[strata == i],
      lower = s$lower[strata == i], upper = s$upper[strata == i]))
  }
  out
}

# ── Data ────────────────────────────────────────────────────────────────────

lung_ecog <- subset(lung, ph.ecog %in% 0:2)

rounding <- matrix(c(
  5.40, 5.50, 5.55, 5.85, 5.70, 5.75, 5.20, 5.60, 5.50, 5.55, 5.50, 5.40,
  5.90, 5.85, 5.70, 5.45, 5.55, 5.60, 5.40, 5.40, 5.35, 5.45, 5.50, 5.35,
  5.25, 5.15, 5.00, 5.85, 5.80, 5.70, 5.25, 5.20, 5.10, 5.65, 5.55, 5.45,
  5.60, 5.35, 5.45, 5.05, 5.00, 4.95, 5.50, 5.50, 5.40, 5.45, 5.55, 5.50,
  5.55, 5.55, 5.35, 5.45, 5.50, 5.55, 5.50, 5.45, 5.25, 5.65, 5.60, 5.40,
  5.70, 5.65, 5.55, 6.30, 6.30, 6.25),
  nrow = 22, byrow = TRUE, dimnames = list(1:22, c("Round Out", "Narrow Angle", "Wide Angle")))

# Hollander & Wolfe (1973) examples used in R's help pages.
mucociliary <- list(
  y = c(2.9, 3.0, 2.5, 2.6, 3.2, 3.8, 2.7, 4.0, 2.4, 2.8, 3.4, 3.7, 2.2, 2.0),
  g = factor(rep(c("Normal", "Obstructive", "Asbestosis"), c(5, 4, 5)),
             levels = c("Normal", "Obstructive", "Asbestosis")))
hw_x <- c(0.80, 0.83, 1.89, 1.04, 1.45, 1.38, 1.91, 1.64, 0.73, 1.46)
hw_y <- c(1.15, 0.88, 0.90, 0.74, 1.21)
depression_x <- c(1.83, 0.50, 1.62, 2.48, 1.68, 1.88, 1.55, 3.06, 1.30)
depression_y <- c(0.878, 0.647, 0.598, 2.05, 1.06, 1.29, 1.06, 3.14, 1.29)

insects_small <- droplevels(subset(InsectSprays, spray %in% c("C", "D", "E")))
insects_small <- do.call(rbind, lapply(split(insects_small, insects_small$spray), head, 4))
plant <- PlantGrowth
theoph1 <- subset(Theoph, Subject == "1")

ref <- list(
  meta = list(
    r = R.version.string,
    packages = lapply(c(multcomp = "multcomp", PMCMRplus = "PMCMRplus", coin = "coin",
                        DescTools = "DescTools", effectsize = "effectsize", car = "car",
                        onewaytests = "onewaytests", survival = "survival"),
                      function(p) as.character(packageVersion(p)))
  ),
  descriptive = list(
    describe("PlantGrowth: ctrl", plant$weight[plant$group == "ctrl"]),
    describe("chickwts: casein", chickwts$weight[chickwts$feed == "casein"]),
    describe("InsectSprays: C", InsectSprays$count[InsectSprays$spray == "C"]),
    describe("precip", unname(precip))
  ),
  ordinaryAnova = list(
    ordinary("PlantGrowth", plant$weight, plant$group),
    ordinary("chickwts", chickwts$weight, chickwts$feed),
    ordinary("InsectSprays", InsectSprays$count, InsectSprays$spray)
  ),
  welchAnova = list(
    welch("PlantGrowth", plant$weight, plant$group),
    welch("InsectSprays", InsectSprays$count, InsectSprays$spray)
  ),
  kruskalWallis = list(
    kruskal("Hollander & Wolfe (1973): mucociliary efficiency", mucociliary$y, mucociliary$g, exact = TRUE),
    kruskal("InsectSprays: C, D, E [1:4]", insects_small$count, insects_small$spray, exact = TRUE),
    kruskal("InsectSprays", InsectSprays$count, InsectSprays$spray)
  ),
  friedman = list(
    friedman("RoundingTimes", rounding),
    friedman("RoundingTimes [1:7, ]", rounding[1:7, ], exact = TRUE)
  ),
  repeatedMeasures = list(repeated("RoundingTimes", rounding)),
  tTestUnpaired = list(
    ttest_unpaired("PlantGrowth: ctrl, trt1", plant$weight[plant$group == "ctrl"],
                   plant$weight[plant$group == "trt1"]),
    ttest_unpaired("sleep: group 1, 2", sleep$extra[1:10], sleep$extra[11:20])
  ),
  tTestPaired = list(
    ttest_paired("sleep: group 1, 2", sleep$extra[1:10], sleep$extra[11:20]),
    ttest_paired("anorexia: Prewt, Postwt", anorexia$Prewt, anorexia$Postwt)
  ),
  mannWhitney = list(
    mann_whitney("Hollander & Wolfe (1973): x, y", hw_x, hw_y, "exact"),
    mann_whitney("InsectSprays: A, B", InsectSprays$count[InsectSprays$spray == "A"],
                 InsectSprays$count[InsectSprays$spray == "B"], "ties"),
    mann_whitney("faithful: waiting[1:60], waiting[61:130]", faithful$waiting[1:60],
                 faithful$waiting[61:130], "approximate")
  ),
  wilcoxon = list(
    wilcoxon("Hollander & Wolfe (1973): depression scale", depression_x, depression_y, "exact"),
    wilcoxon("sleep: group 1, 2", sleep$extra[1:10], sleep$extra[11:20], "ties"),
    wilcoxon("anorexia: Prewt, Postwt", anorexia$Prewt, anorexia$Postwt, "approximate")
  ),
  areaUnderCurve = list(trapezoid("Theoph: Subject 1", theoph1$Time, theoph1$conc)),
  correlation = list(correlation("cars: speed, dist", cars$speed, cars$dist)),
  survival = list(
    survival_curves("gehan: 6-MP, control", gehan$time, gehan$cens, gehan$treat, km = TRUE),
    survival_curves("aml: Maintained, Nonmaintained", aml$time, aml$status, aml$x),
    survival_curves("lung: ph.ecog 0, 1, 2", lung_ecog$time, lung_ecog$status - 1,
                    paste("ph.ecog", lung_ecog$ph.ecog))
  )
)

out <- file.path(script_dir, "r-reference.json")
write_json(ref, out, auto_unbox = TRUE, digits = NA, pretty = TRUE, na = "null")
cat("wrote", out, "\n")