From 9d12e1fcfd531e2adf0746512773b979d045acc5 Mon Sep 17 00:00:00 2001 From: Codex Date: Sun, 9 Aug 2026 10:15:02 +0000 Subject: [PATCH] Add proportion-standardized Pillai estimates --- NEWS.md | 14 + R/pillai_bhatt.R | 213 ++++++++++- README.md | 29 ++ man/pillai_overlap.Rd | 68 +++- tests/testthat/test-pillai-standardized.R | 431 ++++++++++++++++++++++ 5 files changed, 745 insertions(+), 10 deletions(-) create mode 100644 tests/testthat/test-pillai-standardized.R diff --git a/NEWS.md b/NEWS.md index fef501b..0e0dbe5 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,19 @@ # phontrast 2.3.1 +## Proportion-standardized Pillai estimates + +- Added the opt-in `proportion_standardized = TRUE` mode to + `pillai_overlap()`. It appends the plug-in and unbiased squared Mahalanobis + separation estimates and Becker's proportion-standardized Pillai score while + leaving the default two-field return and raw `pillai` value unchanged. +- Negative unbiased separation estimates follow an explicitly labelled + fallback path: `pillai_eq` is `NA`, `pillai_eq_fallback` is `TRUE`, and the + multiplicatively corrected first term appears only as `d2_fallback` beside + its closed-form upward-bias term. +- Added definedness and fragility diagnostics for the two-category, + no-covariate estimator chain, along with regression vectors, base-R MANOVA + comparisons, closed-form anchors, and null Monte-Carlo acceptance tests. + ## CRAN resubmission fixes - Replaced the non-running `extract_mfcc()` example with a fast, executable diff --git a/R/pillai_bhatt.R b/R/pillai_bhatt.R index c506e3e..d05ddd8 100644 --- a/R/pillai_bhatt.R +++ b/R/pillai_bhatt.R @@ -1,21 +1,96 @@ #' Pillai trace for multivariate overlap #' #' Computes the Pillai-Bartlett trace from a MANOVA of features ~ category. -#' This is a convenience wrapper for comparison with JSD. +#' Optionally adds separation estimates and a proportion-standardized Pillai +#' score for an exactly two-category, no-covariate design. #' #' @param data Data frame. #' @param features Character vector of numeric feature columns. -#' @param category_col String; column giving categories (>= 2). +#' @param category_col String; column giving exactly two categories. +#' @param proportion_standardized Logical; if `TRUE`, append the plug-in and +#' unbiased squared Mahalanobis separation estimates and the +#' proportion-standardized Pillai score. The default, `FALSE`, preserves the +#' original two-field return value. #' -#' @return A list with elements `pillai` and `p_value`. +#' @details +#' With `proportion_standardized = FALSE`, the function returns the ordinary +#' Pillai trace and its p-value exactly as before. Raw Pillai describes the +#' categories in the realized data; the proportion-standardized score targets +#' the balanced-design score for the same estimated underlying separation. +#' +#' For two categories with realized counts \eqn{n_1} and \eqn{n_2}, total +#' \eqn{N}, \eqn{p} features, error degrees of freedom \eqn{\nu_e=N-2}, and +#' \eqn{H=2n_1n_2/N}, the optional estimator chain is +#' \deqn{\hat D^2 = 2\nu_e[V/(1-V)]/H,} +#' \deqn{\tilde D^2 = [(\nu_e-p-1)/\nu_e]\hat D^2 - 2p/H,} +#' followed, when \eqn{\tilde D^2 >= 0}, by +#' \deqn{V_{eq}=\tilde D^2/(4+\tilde D^2).} +#' The unbiased estimator is as given by Lachenbruch and Mickey (1968), and +#' \eqn{V_{eq}} implements Becker's (1986) correction and his two-group +#' multivariate generalization. +#' +#' If \eqn{\tilde D^2 < 0}, `pillai_eq` is `NA`, +#' `pillai_eq_fallback` is `TRUE`, and `d2_fallback` contains only the +#' multiplicatively corrected first term. This common near-merger outcome is +#' labelled separately because the fallback retains the split-dependent bias +#' `bias_2p_over_H` and is not a balanced-design equivalent. +#' +#' The estimator chain assumes multivariate normality within each category and +#' a common within-category covariance. It supports no covariates. Both +#' categories must contain at least two complete tokens, and the within-class +#' error SSCP must be nonsingular. When \eqn{\nu_e-p-1 <= 0}, all optional +#' fields are returned as typed `NA` values with a warning. A minority category +#' with fewer than \eqn{p+1} tokens remains computable but is flagged by +#' `fragile_minority`. +#' +#' Because \eqn{x/(4+x)} is strictly concave, `pillai_eq` is slightly downward +#' biased for the balanced-design target. This bias grows with the sampling +#' variance of \eqn{\tilde D^2} and is largest for small, imbalanced, +#' near-merged samples. +#' +#' @return A list with elements `pillai` and `p_value`. When +#' `proportion_standardized = TRUE`, the list additionally contains `n1`, +#' `n2`, `H`, `d2_plugin`, `d2_unbiased`, `pillai_eq`, +#' `pillai_eq_fallback`, `d2_fallback`, `bias_2p_over_H`, and +#' `fragile_minority`. Counts follow the category factor-level order used by +#' the model. On the fallback path, `pillai_eq` is `NA` and +#' `d2_fallback` is labelled separately; off that path, `d2_fallback` is +#' `NA`. +#' +#' @references +#' Becker, G. (1986). Correcting the point-biserial correlation for attenuation +#' owing to unequal sample size. *Journal of Experimental Education*, 55(1), +#' 5-8. +#' +#' Lachenbruch, P. A., & Mickey, M. R. (1968). Estimation of error rates in +#' discriminant analysis. *Technometrics*, 10(1), 1-11. +#' +#' Berry, G. M. (2026). Beyond the null: Calibration, balance, and the +#' interpretation of Pillai scores. Preprint. #' @export #' @importFrom stats manova cov -pillai_overlap <- function(data, features, category_col) { +pillai_overlap <- function(data, + features, + category_col, + proportion_standardized = FALSE) { + .check_bool(proportion_standardized, "proportion_standardized") .check_columns(data, c(category_col, features)) data <- .metric_data(data, c(category_col, features)) .check_numeric_features(data, features) .two_levels(data[[category_col]], "category_col") + standardized_counts <- NULL + if (proportion_standardized) { + standardized_counts <- .pillai_category_counts(data[[category_col]]) + if (any(standardized_counts < 2L)) { + stop( + "Proportion-standardized Pillai estimates require at least two ", + "complete tokens in each category.", + call. = FALSE + ) + } + } + if (length(features) == 1L) { model_df <- data.frame( y = data[[features]], @@ -31,12 +106,27 @@ pillai_overlap <- function(data, features, category_col) { ) } ) + if (proportion_standardized) { + .check_pillai_error_sscp(fit, length(features)) + } aov_tab <- stats::anova(fit) ss_effect <- aov_tab[1, "Sum Sq"] ss_resid <- aov_tab[nrow(aov_tab), "Sum Sq"] - return(list( + result <- list( pillai = as.numeric(ss_effect / (ss_effect + ss_resid)), p_value = as.numeric(aov_tab[1, "Pr(>F)"]) + ) + if (!proportion_standardized) { + return(result) + } + return(c( + result, + .pillai_standardized_fields( + result$pillai, + standardized_counts[[1L]], + standardized_counts[[2L]], + length(features) + ) )) } @@ -58,6 +148,9 @@ pillai_overlap <- function(data, features, category_col) { ) } ) + if (proportion_standardized) { + .check_pillai_error_sscp(m, length(features)) + } s <- tryCatch( summary(m, test = "Pillai"), error = function(e) { @@ -69,10 +162,118 @@ pillai_overlap <- function(data, features, category_col) { ) } ) - list( + result <- list( pillai = s$stats[1, "Pillai"], p_value = s$stats[1, "Pr(>F)"] ) + if (!proportion_standardized) { + return(result) + } + c( + result, + .pillai_standardized_fields( + result$pillai, + standardized_counts[[1L]], + standardized_counts[[2L]], + length(features) + ) + ) +} + +.pillai_category_counts <- function(category) { + category <- category[!is.na(category)] + if (is.factor(category)) { + category <- droplevels(category) + } else { + category <- factor(category) + } + as.integer(table(category)) +} + +.check_pillai_error_sscp <- function(fit, p) { + residual_matrix <- as.matrix(stats::residuals(fit)) + error_sscp <- crossprod(residual_matrix) + factorization <- tryCatch( + chol(error_sscp), + error = function(e) NULL + ) + if (nrow(error_sscp) != p || ncol(error_sscp) != p || + any(!is.finite(error_sscp)) || is.null(factorization)) { + stop( + "Proportion-standardized Pillai estimates require a nonsingular ", + "within-class error SSCP matrix; no ridge correction is applied.", + call. = FALSE + ) + } + invisible(TRUE) +} + +.pillai_standardized_na_fields <- function() { + list( + n1 = NA_integer_, + n2 = NA_integer_, + H = NA_real_, + d2_plugin = NA_real_, + d2_unbiased = NA_real_, + pillai_eq = NA_real_, + pillai_eq_fallback = NA, + d2_fallback = NA_real_, + bias_2p_over_H = NA_real_, + fragile_minority = NA + ) +} + +.pillai_standardized_fields <- function(V, n1, n2, p) { + if (!is.numeric(V) || length(V) != 1L || !is.finite(V) || V < 0 || V >= 1) { + stop("`V` must be a single finite Pillai score in [0, 1).", call. = FALSE) + } + .check_positive_count(n1, "n1") + .check_positive_count(n2, "n2") + .check_positive_count(p, "p") + if (n1 < 2L || n2 < 2L) { + stop("`n1` and `n2` must each be at least 2.", call. = FALSE) + } + + n1 <- as.integer(n1) + n2 <- as.integer(n2) + p <- as.integer(p) + N <- n1 + n2 + nu_e <- N - 2 + if (nu_e - p - 1 <= 0) { + warning( + "Proportion-standardized Pillai estimates are undefined because ", + "nu_e - p - 1 <= 0 (N = ", N, ", p = ", p, "); all optional ", + "fields were returned as NA.", + call. = FALSE + ) + return(.pillai_standardized_na_fields()) + } + + H <- 2 * n1 * n2 / N + pi1 <- n1 / N + theta_hat <- V / (1 - V) + d2_plugin <- nu_e * theta_hat / (N * pi1 * (1 - pi1)) + first_term <- (nu_e - p - 1) / nu_e * d2_plugin + bias_2p_over_H <- 2 * p / H + d2_unbiased <- first_term - bias_2p_over_H + pillai_eq_fallback <- d2_unbiased < 0 + + list( + n1 = n1, + n2 = n2, + H = H, + d2_plugin = d2_plugin, + d2_unbiased = d2_unbiased, + pillai_eq = if (pillai_eq_fallback) { + NA_real_ + } else { + d2_unbiased / (4 + d2_unbiased) + }, + pillai_eq_fallback = pillai_eq_fallback, + d2_fallback = if (pillai_eq_fallback) first_term else NA_real_, + bias_2p_over_H = bias_2p_over_H, + fragile_minority = min(n1, n2) < p + 1L + ) } #' Group-level Pillai scores diff --git a/README.md b/README.md index 6b5db48..68a7218 100644 --- a/README.md +++ b/README.md @@ -217,6 +217,35 @@ Use lower-level helpers such as `jsd_kde_nd()`, `percent_overlap_kde()`, `pillai_overlap()`, and `bhattacharyya_mvnorm()` when you are validating a method, debugging one contrast, or need direct control over one metric. +### Proportion-standardized Pillai estimates + +For an exactly two-category design, `pillai_overlap()` can optionally estimate +the squared Mahalanobis separation and the Pillai score the same estimated +system would have under a balanced category split: + +```r +pillai_overlap( + data = vowels, + features = c("f1", "f2"), + category_col = "vowel", + proportion_standardized = TRUE +) +``` + +The default remains `FALSE`, so existing calls still return only `pillai` and +`p_value`. Raw Pillai describes the categories in the realized speech data; +the opt-in `pillai_eq` estimate characterizes the corresponding balanced-design +system. The extended result also reports `d2_plugin`, `d2_unbiased`, the class +counts and harmonic mean, and the finite-sample bias term. + +When the unbiased separation estimate is negative, `pillai_eq` is `NA` and the +result sets `pillai_eq_fallback = TRUE`. The nonnegative fallback is reported +only as `d2_fallback`, alongside its known upward bias; it is never silently +mapped onto the corrected-score scale. This is a common outcome near merger, +not an exceptional error case. These estimates assume multivariate normality +within each category and a common within-category covariance, and no covariates +or multicategory extension is performed. + > `compare_overlap_metrics()` from phonJSD still works but is **deprecated**: it > now calls `phontrast()` with `output = "wide"`. Switch calls to `phontrast()`. diff --git a/man/pillai_overlap.Rd b/man/pillai_overlap.Rd index c6c46fc..48639f9 100644 --- a/man/pillai_overlap.Rd +++ b/man/pillai_overlap.Rd @@ -4,19 +4,79 @@ \alias{pillai_overlap} \title{Pillai trace for multivariate overlap} \usage{ -pillai_overlap(data, features, category_col) +pillai_overlap(data, features, category_col, proportion_standardized = FALSE) } \arguments{ \item{data}{Data frame.} \item{features}{Character vector of numeric feature columns.} -\item{category_col}{String; column giving categories (>= 2).} +\item{category_col}{String; column giving exactly two categories.} + +\item{proportion_standardized}{Logical; if \code{TRUE}, append the plug-in and +unbiased squared Mahalanobis separation estimates and the +proportion-standardized Pillai score. The default, \code{FALSE}, preserves the +original two-field return value.} } \value{ -A list with elements \code{pillai} and \code{p_value}. +A list with elements \code{pillai} and \code{p_value}. When +\code{proportion_standardized = TRUE}, the list additionally contains \code{n1}, +\code{n2}, \code{H}, \code{d2_plugin}, \code{d2_unbiased}, \code{pillai_eq}, +\code{pillai_eq_fallback}, \code{d2_fallback}, \code{bias_2p_over_H}, and +\code{fragile_minority}. Counts follow the category factor-level order used by +the model. On the fallback path, \code{pillai_eq} is \code{NA} and +\code{d2_fallback} is labelled separately; off that path, \code{d2_fallback} is +\code{NA}. } \description{ Computes the Pillai-Bartlett trace from a MANOVA of features ~ category. -This is a convenience wrapper for comparison with JSD. +Optionally adds separation estimates and a proportion-standardized Pillai +score for an exactly two-category, no-covariate design. +} +\details{ +With \code{proportion_standardized = FALSE}, the function returns the ordinary +Pillai trace and its p-value exactly as before. Raw Pillai describes the +categories in the realized data; the proportion-standardized score targets +the balanced-design score for the same estimated underlying separation. + +For two categories with realized counts \eqn{n_1} and \eqn{n_2}, total +\eqn{N}, \eqn{p} features, error degrees of freedom \eqn{\nu_e=N-2}, and +\eqn{H=2n_1n_2/N}, the optional estimator chain is +\deqn{\hat D^2 = 2\nu_e[V/(1-V)]/H,} +\deqn{\tilde D^2 = [(\nu_e-p-1)/\nu_e]\hat D^2 - 2p/H,} +followed, when \eqn{\tilde D^2 >= 0}, by +\deqn{V_{eq}=\tilde D^2/(4+\tilde D^2).} +The unbiased estimator is as given by Lachenbruch and Mickey (1968), and +\eqn{V_{eq}} implements Becker's (1986) correction and his two-group +multivariate generalization. + +If \eqn{\tilde D^2 < 0}, \code{pillai_eq} is \code{NA}, +\code{pillai_eq_fallback} is \code{TRUE}, and \code{d2_fallback} contains only the +multiplicatively corrected first term. This common near-merger outcome is +labelled separately because the fallback retains the split-dependent bias +\code{bias_2p_over_H} and is not a balanced-design equivalent. + +The estimator chain assumes multivariate normality within each category and +a common within-category covariance. It supports no covariates. Both +categories must contain at least two complete tokens, and the within-class +error SSCP must be nonsingular. When \eqn{\nu_e-p-1 <= 0}, all optional +fields are returned as typed \code{NA} values with a warning. A minority category +with fewer than \eqn{p+1} tokens remains computable but is flagged by +\code{fragile_minority}. + +Because \eqn{x/(4+x)} is strictly concave, \code{pillai_eq} is slightly downward +biased for the balanced-design target. This bias grows with the sampling +variance of \eqn{\tilde D^2} and is largest for small, imbalanced, +near-merged samples. +} +\references{ +Becker, G. (1986). Correcting the point-biserial correlation for attenuation +owing to unequal sample size. \emph{Journal of Experimental Education}, 55(1), +5-8. + +Lachenbruch, P. A., & Mickey, M. R. (1968). Estimation of error rates in +discriminant analysis. \emph{Technometrics}, 10(1), 1-11. + +Berry, G. M. (2026). Beyond the null: Calibration, balance, and the +interpretation of Pillai scores. Preprint. } diff --git a/tests/testthat/test-pillai-standardized.R b/tests/testthat/test-pillai-standardized.R new file mode 100644 index 0000000..17d361a --- /dev/null +++ b/tests/testthat/test-pillai-standardized.R @@ -0,0 +1,431 @@ +expect_relative_equal <- function(actual, expected, tolerance = 1e-9) { + relative_error <- abs(actual - expected) / + max(abs(expected), .Machine$double.eps) + expect_true( + relative_error <= tolerance, + info = paste0( + "relative error ", format(relative_error, digits = 16), + " exceeds ", format(tolerance, scientific = TRUE) + ) + ) +} + +test_that("default Pillai return and base-R value remain unchanged", { + set.seed(101) + data <- data.frame( + category = rep(c("a", "b"), each = 30), + f1 = c(rnorm(30), rnorm(30, 0.8)), + f2 = c(rnorm(30), rnorm(30, 0.3)) + ) + + default <- pillai_overlap(data, c("f1", "f2"), "category") + explicit_default <- pillai_overlap( + data, + c("f1", "f2"), + "category", + proportion_standardized = FALSE + ) + base_v <- summary( + stats::manova(cbind(f1, f2) ~ category, data = data), + test = "Pillai" + )$stats[1, "Pillai"] + + expect_named(default, c("pillai", "p_value")) + expect_identical(default, explicit_default) + expect_equal(default$pillai, base_v, tolerance = 1e-6) + + one_feature <- pillai_overlap(data, "f1", "category") + one_fit <- stats::lm(f1 ~ category, data = data) + one_aov <- stats::anova(one_fit) + expected_one <- one_aov[1, "Sum Sq"] / sum(one_aov[, "Sum Sq"]) + expect_named(one_feature, c("pillai", "p_value")) + expect_equal(one_feature$pillai, expected_one, tolerance = 1e-12) + + one_standardized <- pillai_overlap( + data, + "f1", + "category", + proportion_standardized = TRUE + ) + one_expected <- getFromNamespace(".pillai_standardized_fields", "phontrast")( + one_feature$pillai, 30, 30, 1 + ) + expect_identical(one_standardized[c("pillai", "p_value")], one_feature) + expect_equal( + one_standardized[setdiff(names(one_standardized), c("pillai", "p_value"))], + one_expected, + tolerance = 1e-12 + ) +}) + +test_that("worked estimator-chain vectors are reproduced", { + chain <- getFromNamespace(".pillai_standardized_fields", "phontrast") + vectors <- data.frame( + id = paste0("TV", 1:8), + n1 = c(50, 50, 20, 10, 5, 20, 2, 30), + n2 = c(50, 50, 80, 90, 95, 80, 58, 30), + p = c(2, 2, 2, 2, 2, 3, 2, 2), + V = c(0.15, 0.02, 0.12, 0.02, 0.15, 0.12, 0.25, 0.035088), + H = c(50, 50, 32, 18, 9.5, 32, 3.8666666666666667, 30), + bias = c( + 0.08, 0.08, 0.125, 0.2222222222222222, + 0.42105263157894735, 0.1875, 1.0344827586206897, + 0.13333333333333333 + ), + d2_plugin = c( + 0.6917647058823531, 0.08, 0.8352272727272727, + 0.2222222222222222, 3.640866873065016, 0.8352272727272727, + 10, 0.14060722635846584 + ), + d2_unbiased = c( + 0.5905882352941177, -0.0024489795918367363, + 0.6846590909090908, -0.006802721088435382, + 3.108359133126936, 0.6136363636363635, + 8.448275862068964, 0.0000011054548670341724 + ), + fallback = c(FALSE, TRUE, FALSE, TRUE, FALSE, FALSE, FALSE, FALSE), + pillai_eq = c( + 0.12865197334700154, NA, 0.14614918132201332, NA, + 0.43728222996515687, 0.13300492610837436, + 0.6786703601108033, 2.763636403816603e-7 + ), + d2_fallback = c( + NA, 0.07755102040816327, NA, 0.21541950113378683, + NA, NA, NA, NA + ), + fragile = c(FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, FALSE) + ) + + expected_names <- c( + "n1", "n2", "H", "d2_plugin", "d2_unbiased", "pillai_eq", + "pillai_eq_fallback", "d2_fallback", "bias_2p_over_H", + "fragile_minority" + ) + + for (i in seq_len(nrow(vectors))) { + v <- vectors[i, ] + out <- chain(v$V, v$n1, v$n2, v$p) + + expect_identical(names(out), expected_names, info = v$id) + expect_identical(out$n1, as.integer(v$n1), info = v$id) + expect_identical(out$n2, as.integer(v$n2), info = v$id) + expect_relative_equal(out$H, v$H) + expect_relative_equal(out$bias_2p_over_H, v$bias) + expect_relative_equal(out$d2_plugin, v$d2_plugin) + expect_relative_equal(out$d2_unbiased, v$d2_unbiased) + expect_identical(out$pillai_eq_fallback, v$fallback, info = v$id) + expect_identical(out$fragile_minority, v$fragile, info = v$id) + + if (is.na(v$pillai_eq)) { + expect_true(is.na(out$pillai_eq), info = v$id) + } else { + expect_relative_equal(out$pillai_eq, v$pillai_eq) + } + if (is.na(v$d2_fallback)) { + expect_true(is.na(out$d2_fallback), info = v$id) + } else { + expect_relative_equal(out$d2_fallback, v$d2_fallback) + } + } + + boundary <- chain(0.035088, 30, 30, 2) + expect_gt(boundary$d2_unbiased, 0) + expect_false(boundary$pillai_eq_fallback) +}) + +test_that("public opt-in appends fields and labels fallback results", { + points <- data.frame( + f1 = c(0, 1, 0, 1), + f2 = c(0, 0, 1, 1) + ) + data <- rbind( + data.frame(category = "a", points), + data.frame(category = "b", points) + ) + + raw <- pillai_overlap(data, c("f1", "f2"), "category") + standardized <- pillai_overlap( + data, + c("f1", "f2"), + "category", + proportion_standardized = TRUE + ) + + expect_identical(standardized[c("pillai", "p_value")], raw) + expect_named( + standardized, + c( + "pillai", "p_value", "n1", "n2", "H", "d2_plugin", + "d2_unbiased", "pillai_eq", "pillai_eq_fallback", + "d2_fallback", "bias_2p_over_H", "fragile_minority" + ) + ) + expect_true(standardized$pillai_eq_fallback) + expect_true(is.na(standardized$pillai_eq)) + expect_equal(standardized$d2_fallback, 0, tolerance = 1e-15) + expect_lt(standardized$d2_unbiased, 0) + expect_equal( + standardized$d2_unbiased, + standardized$d2_fallback - standardized$bias_2p_over_H, + tolerance = 1e-15 + ) +}) + +test_that("category counts use model factor-level order after filtering", { + set.seed(102) + data <- data.frame( + category = factor( + c(rep("a", 20), rep("b", 80), "a"), + levels = c("b", "a", "unused") + ), + f1 = rnorm(101), + f2 = rnorm(101) + ) + data$f1[101] <- NA_real_ + + out <- pillai_overlap( + data, + c("f1", "f2"), + "category", + proportion_standardized = TRUE + ) + + expect_identical(out$n1, 80L) + expect_identical(out$n2, 20L) + expect_equal(out$H, 32) + expect_false(out$fragile_minority) +}) + +test_that("definedness and design failures are explicit", { + undefined <- data.frame( + category = c("a", "a", "b", "b", "b"), + f1 = c(0, 1, 0, 1, 2), + f2 = c(0, 0, 1, 2, 1) + ) + expect_warning( + out <- pillai_overlap( + undefined, + c("f1", "f2"), + "category", + proportion_standardized = TRUE + ), + "undefined" + ) + optional <- setdiff(names(out), c("pillai", "p_value")) + expect_true(all(vapply(out[optional], function(x) length(x) == 1L && is.na(x), + logical(1)))) + expect_type(out$n1, "integer") + expect_type(out$H, "double") + expect_type(out$pillai_eq_fallback, "logical") + + too_small <- data.frame( + category = c("a", rep("b", 5)), + f1 = c(0, 1:5) + ) + expect_error( + pillai_overlap( + too_small, + "f1", + "category", + proportion_standardized = TRUE + ), + "at least two" + ) + + singular <- data.frame( + category = rep(c("a", "b"), each = 5), + f1 = c(0:4, 1:5), + f2 = 2 * c(0:4, 1:5) + ) + expect_error( + pillai_overlap( + singular, + c("f1", "f2"), + "category", + proportion_standardized = TRUE + ), + "nonsingular" + ) + + three_classes <- data.frame( + category = rep(c("a", "b", "c"), each = 4), + f1 = seq_len(12) + ) + expect_error( + pillai_overlap( + three_classes, + "f1", + "category", + proportion_standardized = TRUE + ), + "exactly two" + ) + expect_error( + pillai_overlap( + undefined, + c("f1", "f2"), + "category", + proportion_standardized = NA + ), + "TRUE or FALSE" + ) +}) + +test_that("internal chain guards invalid scores and sample sizes", { + chain <- getFromNamespace(".pillai_standardized_fields", "phontrast") + + expect_error(chain(-0.01, 20, 20, 2), "\\[0, 1\\)") + expect_error(chain(1, 20, 20, 2), "\\[0, 1\\)") + expect_error(chain(Inf, 20, 20, 2), "\\[0, 1\\)") + expect_error(chain(NA_real_, 20, 20, 2), "\\[0, 1\\)") + expect_error(chain(0.2, 1, 20, 2), "at least 2") + expect_error(chain(0.2, 20.5, 20, 2), "positive integer") +}) + +test_that("end-to-end values agree with base-R MANOVA and the chain", { + set.seed(2026) + n1 <- 20L + n2 <- 80L + A <- cbind(rnorm(n1, 0, 1), rnorm(n1, 0, 1)) + B <- cbind(rnorm(n2, 0.9, 1), rnorm(n2, 0.4, 1)) + data <- data.frame( + class = rep(c("a", "b"), c(n1, n2)), + F1 = c(A[, 1], B[, 1]), + F2 = c(A[, 2], B[, 2]) + ) + + base_v <- summary( + stats::manova(cbind(F1, F2) ~ class, data = data), + test = "Pillai" + )$stats[1, "Pillai"] + out <- pillai_overlap( + data, + c("F1", "F2"), + "class", + proportion_standardized = TRUE + ) + expected <- getFromNamespace(".pillai_standardized_fields", "phontrast")( + base_v, n1, n2, 2 + ) + + expect_equal(out$pillai, base_v, tolerance = 1e-6) + expect_equal(base_v, 0.205456307, tolerance = 1e-6) + expect_relative_equal(out$d2_plugin, expected$d2_plugin, tolerance = 1e-12) + expect_relative_equal(out$d2_unbiased, expected$d2_unbiased, tolerance = 1e-12) + expect_identical(out$pillai_eq_fallback, FALSE) + expect_relative_equal(out$pillai_eq, expected$pillai_eq, tolerance = 1e-12) + expect_equal(out$d2_plugin, 1.58382716, tolerance = 1e-8) + expect_equal(out$d2_unbiased, 1.41034265, tolerance = 1e-8) + expect_equal(out$pillai_eq, 0.260675293, tolerance = 1e-8) +}) + +test_that("closed-form distribution and expectation anchors hold", { + expect_equal( + stats::qbeta(0.95, 1, 28.5), + 0.09977758022528305, + tolerance = 1e-14 + ) + expect_equal( + stats::qbeta(0.95, 1, 48.5), + 0.05989873024908154, + tolerance = 1e-14 + ) + expect_equal( + 1.44 / (4 + 1.44), + 0.2647058823529411, + tolerance = 1e-14 + ) + + fallback_probability <- function(N, p) { + nu_e <- N - 2 + stats::pf( + (nu_e - p + 1) / (nu_e - p - 1), + p, + nu_e - p + 1 + ) + } + expect_equal(fallback_probability(40, 2), 0.642293423, tolerance = 1e-9) + expect_equal(fallback_probability(100, 2), 0.635946052, tolerance = 1e-9) + expect_equal(fallback_probability(500, 2), 0.632862003, tolerance = 1e-9) + expect_equal(1 - 1 / exp(1), 0.632120559, tolerance = 1e-9) + expect_equal(fallback_probability(100, 3), 0.613214289, tolerance = 1e-9) + + expected_plugin_at_merger <- function(H) { + N <- 100 + p <- 2 + nu_e <- N - 2 + 2 * nu_e * p / (H * (nu_e - p - 1)) + } + expect_equal( + expected_plugin_at_merger(50), + 0.08252631578947368, + tolerance = 1e-14 + ) + expect_equal( + expected_plugin_at_merger(9.5), + 0.4343490304709141, + tolerance = 1e-14 + ) +}) + +test_that("merger fallback rate follows the split-independent law", { + chain <- getFromNamespace(".pillai_standardized_fields", "phontrast") + replicates <- 4000L + target <- stats::pf(97 / 95, 2, 97) + fallback_tolerance <- 3 * sqrt(target * (1 - target) / replicates) + size_tolerance <- 3 * sqrt(0.05 * 0.95 / replicates) + critical_v <- stats::qbeta(0.95, 1, 48.5) + + simulate_split <- function(n1, n2) { + N <- n1 + n2 + nu_e <- N - 2 + fallback <- logical(replicates) + d2_unbiased <- numeric(replicates) + reject <- logical(replicates) + + for (i in seq_len(replicates)) { + A <- matrix(rnorm(n1 * 2), ncol = 2) + B <- matrix(rnorm(n2 * 2), ncol = 2) + mean_A <- colMeans(A) + mean_B <- colMeans(B) + centered_A <- sweep(A, 2, mean_A) + centered_B <- sweep(B, 2, mean_B) + error_sscp <- crossprod(centered_A) + crossprod(centered_B) + delta <- mean_A - mean_B + t2 <- n1 * n2 / N * drop( + crossprod(delta, solve(error_sscp / nu_e, delta)) + ) + V <- t2 / (t2 + nu_e) + estimate <- chain(V, n1, n2, 2) + + fallback[i] <- estimate$pillai_eq_fallback + d2_unbiased[i] <- estimate$d2_unbiased + reject[i] <- V > critical_v + } + + list( + fallback_rate = mean(fallback), + d2_mean = mean(d2_unbiased), + d2_se = stats::sd(d2_unbiased) / sqrt(replicates), + size = mean(reject) + ) + } + + # Fixed R acceptance seed; do not change without recording new diagnostics. + set.seed(2026) + results <- list( + balanced = simulate_split(50L, 50L), + imbalanced = simulate_split(5L, 95L) + ) + + for (name in names(results)) { + result <- results[[name]] + expect_lt( + abs(result$fallback_rate - target), + fallback_tolerance + ) + expect_gt(result$fallback_rate, 0.5) + expect_lt(abs(result$d2_mean), 3 * result$d2_se) + expect_lt(abs(result$size - 0.05), size_tolerance) + } +})