From 7185242ee0c09347477b10cfa3ae0148efef2b47 Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sat, 25 Jul 2026 21:09:40 +0900 Subject: [PATCH] =?UTF-8?q?feat:=20generalizability=20theory=20G/D=20studi?= =?UTF-8?q?es=20(crossed=20p=20=C3=97=20i,=20p=20=C3=97=20i=20=C3=97=20o)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Rebased onto main after #228 squash merge. Co-authored-by: Copilot App <223556219+Copilot@users.noreply.github.com> --- CHANGELOG.md | 16 + crates/fast-mlsirm-py/src/lib.rs | 77 ++++- crates/mlsirm-core/src/classification.rs | 48 +-- crates/mlsirm-core/src/gtheory.rs | 419 +++++++++++++++++++++++ crates/mlsirm-core/src/lib.rs | 1 + python/fast_mlsirm/__init__.py | 10 + python/fast_mlsirm/classification.py | 32 +- python/fast_mlsirm/gtheory.py | 144 ++++++++ tests/test_paper_features.py | 86 ++++- tests/unit/classification_tests.rs | 22 +- tests/unit/gtheory_tests.rs | 381 +++++++++++++++++++++ 11 files changed, 1170 insertions(+), 66 deletions(-) create mode 100644 crates/mlsirm-core/src/gtheory.rs create mode 100644 python/fast_mlsirm/gtheory.py create mode 100644 tests/unit/gtheory_tests.rs diff --git a/CHANGELOG.md b/CHANGELOG.md index dbc5dd205..205d27b21 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -114,6 +114,22 @@ ### Added +- **Generalizability theory G/D studies for crossed designs** + (`fast_mlsirm.gtheory_pi`, `fast_mlsirm.gtheory_pio`; in + `mlsirm_core::gtheory`; Huebner & Lucht, 2019, read in full — the EMS + inversions the paper defers to Brennan, 2001, and Shavelson & Webb, + 1991, both cited-not-read, are hand-derived and numerically verified + against the paper's published Tables 3-6). One-facet `p x i` and + two-facet `p x i x o` random-effects ANOVA variance components, plus + D-study relative/absolute error variances, generalizability coefficient + E-rho^2, and dependability Phi over proposed facet sizes. Negative raw + components are reported as-is in `var_raw` and clamped to zero in `var` + for the D study (documented clamped-ANOVA implementation policy); + coefficients are NaN when their denominator is <= 1e-12. Rust tests + reproduce the paper's worked examples at full precision, add + independent RNG-pinned fixtures (including a natural negative-component + anchor), executed mutation kills (M1/M3/M5/M6), and a 500-rep + `#[ignore]` Monte-Carlo recovery test. - **Livingston-Lewis classification accuracy and consistency** (`fast_mlsirm.livingston_lewis`; in `mlsirm_core::classification`; Livingston & Lewis, 1995, as implemented in Haakstad's CRAN diff --git a/crates/fast-mlsirm-py/src/lib.rs b/crates/fast-mlsirm-py/src/lib.rs index 77d64e87c..ea0d54157 100644 --- a/crates/fast-mlsirm-py/src/lib.rs +++ b/crates/fast-mlsirm-py/src/lib.rs @@ -47,6 +47,9 @@ use mlsirm_core::fitstats::{ }; use mlsirm_core::gpcm::{fit_gpcm as core_fit_gpcm, GpcmConfig}; use mlsirm_core::grm::{fit_grm as core_fit_grm, GrmConfig}; +use mlsirm_core::gtheory::{ + gtheory_pi as core_gtheory_pi, gtheory_pio as core_gtheory_pio, GTheoryDStudyRow, +}; use mlsirm_core::ksirt::{ksirt as core_ksirt, KsirtKernel}; use mlsirm_core::lltm::{fit_lltm as core_fit_lltm, LltmConfig}; use mlsirm_core::mhrm::{fit_mhrm as core_fit_mhrm, MhrmConfig, MhrmModel}; @@ -1788,8 +1791,8 @@ fn livingston_lewis( out.set_item("used_two_parameter", res.used_two_parameter)?; out.set_item("p_tp", res.p_tp)?; out.set_item("p_fp", res.p_fp)?; - out.set_item("p_tn", res.p_tn)?; - out.set_item("p_fn", res.p_fn)?; + out.set_item("p_tf", res.p_tf)?; + out.set_item("p_ff", res.p_ff)?; out.set_item("accuracy", res.accuracy)?; out.set_item("sensitivity", res.sensitivity)?; out.set_item("specificity", res.specificity)?; @@ -1803,6 +1806,74 @@ fn livingston_lewis( Ok(out.into()) } +/// One-facet crossed `p x i` generalizability analysis +/// (`mlsirm_core::gtheory`; Huebner & Lucht, 2019, Tables 3-4). `x` is a +/// flattened row-major `n_p x n_i` score matrix; `n_i_prime` lists the +/// proposed D-study item counts. +#[pyfunction] +fn gtheory_pi( + py: Python<'_>, + x: PyReadonlyArray1<'_, f64>, + n_p: usize, + n_i: usize, + n_i_prime: Vec, +) -> PyResult> { + let res = + core_gtheory_pi(x.as_slice()?, n_p, n_i, &n_i_prime).map_err(PyValueError::new_err)?; + let out = pyo3::types::PyDict::new(py); + out.set_item("df", res.df.to_vec())?; + out.set_item("ss", res.ss.to_vec())?; + out.set_item("ms", res.ms.to_vec())?; + out.set_item("var_raw", res.var_raw.to_vec())?; + out.set_item("var", res.var.to_vec())?; + out.set_item("d_study", d_study_rows_to_py(py, &res.d_study)?)?; + Ok(out.into()) +} + +/// Two-facet crossed `p x i x o` generalizability analysis +/// (`mlsirm_core::gtheory`; Huebner & Lucht, 2019, Tables 5-6). `x` is +/// flattened `x[p*n_i*n_o + i*n_o + o]`; `n_prime` lists proposed +/// `(n_i', n_o')` D-study pairs. Component order in all arrays: +/// (p, i, o, pi, po, io, pio). +#[pyfunction] +fn gtheory_pio( + py: Python<'_>, + x: PyReadonlyArray1<'_, f64>, + n_p: usize, + n_i: usize, + n_o: usize, + n_prime: Vec<(usize, usize)>, +) -> PyResult> { + let res = core_gtheory_pio(x.as_slice()?, n_p, n_i, n_o, &n_prime) + .map_err(PyValueError::new_err)?; + let out = pyo3::types::PyDict::new(py); + out.set_item("df", res.df.to_vec())?; + out.set_item("ss", res.ss.to_vec())?; + out.set_item("ms", res.ms.to_vec())?; + out.set_item("var_raw", res.var_raw.to_vec())?; + out.set_item("var", res.var.to_vec())?; + out.set_item("d_study", d_study_rows_to_py(py, &res.d_study)?)?; + Ok(out.into()) +} + +fn d_study_rows_to_py( + py: Python<'_>, + rows: &[GTheoryDStudyRow], +) -> PyResult>> { + rows.iter() + .map(|r| { + let d = pyo3::types::PyDict::new(py); + d.set_item("n_i_prime", r.n_i_prime)?; + d.set_item("n_o_prime", r.n_o_prime)?; + d.set_item("rel_error_var", r.rel_error_var)?; + d.set_item("abs_error_var", r.abs_error_var)?; + d.set_item("generalizability", r.generalizability)?; + d.set_item("dependability", r.dependability)?; + Ok(d.into()) + }) + .collect() +} + /// Horn's parallel analysis for principal-component retention /// (`mlsirm_core::parallel`; oracle: CRAN paran 1.5.6, PCA path). `data` is /// a flattened row-major `n_persons * n_items` matrix; `centile` is 0 for @@ -5524,6 +5595,8 @@ fn fast_mlsirm_core(m: &Bound<'_, PyModule>) -> PyResult<()> { m.add_function(wrap_pyfunction!(rudner_classification, m)?)?; m.add_function(wrap_pyfunction!(lee_classification, m)?)?; m.add_function(wrap_pyfunction!(livingston_lewis, m)?)?; + m.add_function(wrap_pyfunction!(gtheory_pi, m)?)?; + m.add_function(wrap_pyfunction!(gtheory_pio, m)?)?; m.add_function(wrap_pyfunction!(parallel_analysis, m)?)?; m.add_function(wrap_pyfunction!(guttman_lambdas, m)?)?; m.add_function(wrap_pyfunction!(tenberge_mu, m)?)?; diff --git a/crates/mlsirm-core/src/classification.rs b/crates/mlsirm-core/src/classification.rs index 5bd492319..f53431a83 100644 --- a/crates/mlsirm-core/src/classification.rs +++ b/crates/mlsirm-core/src/classification.rs @@ -349,20 +349,20 @@ pub struct LivingstonLewisResult { pub beta: f64, /// True iff the two-parameter fail-safe replaced the four-parameter fit. pub used_two_parameter: bool, - /// P(true pass, observed pass) — true positive. + /// P(true pass, observed pass). pub p_tp: f64, /// P(true fail, observed pass) — false positive. pub p_fp: f64, - /// P(true fail, observed fail) — true negative. - pub p_tn: f64, + /// P(true fail, observed fail). + pub p_tf: f64, /// P(true pass, observed fail) — false negative. - pub p_fn: f64, - /// `p_tp + p_tn`. + pub p_ff: f64, + /// `p_tp + p_tf`. pub accuracy: f64, - /// `p_tp / (p_tp + p_fn)` — P(observed pass | true pass); `NaN` when + /// `p_tp / (p_tp + p_ff)` — P(observed pass | true pass); `NaN` when /// the true-pass margin vanishes (cut outside the fitted support). pub sensitivity: f64, - /// `p_tn / (p_tn + p_fp)` — P(observed fail | true fail); `NaN` when + /// `p_tf / (p_tf + p_fp)` — P(observed fail | true fail); `NaN` when /// the true-fail margin vanishes. pub specificity: f64, /// Consistency cell: fail on both administrations. @@ -417,15 +417,6 @@ fn gauss_legendre(n: usize) -> (Vec, Vec) { (x, w) } -/// Cached 64-node Gauss-Legendre rule; computed once and reused across all -/// `beta_weighted_integral` calls within a process. -static GL64: std::sync::OnceLock<(Vec, Vec)> = std::sync::OnceLock::new(); - -#[inline] -fn gl64() -> &'static (Vec, Vec) { - GL64.get_or_init(|| gauss_legendre(64)) -} - /// `integral_{t0}^{t1} t^(a-1) (1-t)^(b-1) g(t) dt / B(a, b)`. Endpoint /// singularities (shape < 1) are absorbed exactly by power substitutions /// (spec rev 2): on `t <= 1/2` with `a < 1` substitute `v = t^a` so @@ -440,14 +431,13 @@ fn gl64() -> &'static (Vec, Vec) { fn beta_weighted_integral(a: f64, b: f64, t0: f64, t1: f64, g: impl Fn(f64) -> f64) -> f64 { let ln_b = crate::fitstats::ln_gamma(a) + crate::fitstats::ln_gamma(b) - crate::fitstats::ln_gamma(a + b); - let rule = gl64(); - let (nodes, weights) = (&rule.0, &rule.1); + let (nodes, weights) = gauss_legendre(64); let gl = |lo: f64, hi: f64, f: &dyn Fn(f64) -> f64| -> f64 { let c = 0.5 * (lo + hi); let h = 0.5 * (hi - lo); nodes .iter() - .zip(weights.iter()) + .zip(&weights) .map(|(z, wt)| wt * f(c + h * z)) .sum::() * h @@ -551,20 +541,20 @@ fn falling_factorial(a: f64, r: u32) -> f64 { /// /// # References /// -/// Haakstad, H. (2022). *betafunctions: Functions for working with two- and +/// Haakstad, H. (2023). *betafunctions: Functions for working with two- and /// four-parameter beta probability distributions* (Version 1.9.0) /// \[R package\]. CRAN. /// /// Hanson, B. A. (1991). *Method of moments estimates for the four-parameter /// beta compound binomial model and the calculation of classification /// consistency indexes* (ACT Research Report 91-5). (As cited in Haakstad, -/// 2022; not read.) +/// 2023; not read.) /// /// Livingston, S. A., & Lewis, C. (1995). Estimating the consistency and /// accuracy of classifications based on test scores. *Journal of Educational /// Measurement, 32*(2), 179-197. /// https://doi.org/10.1111/j.1745-3984.1995.tb00462.x (As implemented in -/// Haakstad, 2022; the paper itself was not obtainable.) +/// Haakstad, 2023; the paper itself was not obtainable.) /// /// Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. /// (2007). *Numerical recipes: The art of scientific computing* (3rd ed.). @@ -689,8 +679,8 @@ pub fn livingston_lewis( let tc = ((c - lower) / (upper - lower)).clamp(0.0, 1.0); let p_tp = beta_weighted_integral(a, b, tc, 1.0, |t| 1.0 - fail_prob(t)); let p_fp = beta_weighted_integral(a, b, 0.0, tc, |t| 1.0 - fail_prob(t)); - let p_fn = beta_weighted_integral(a, b, tc, 1.0, &fail_prob); - let p_tn = beta_weighted_integral(a, b, 0.0, tc, &fail_prob); + let p_ff = beta_weighted_integral(a, b, tc, 1.0, &fail_prob); + let p_tf = beta_weighted_integral(a, b, 0.0, tc, &fail_prob); let p_ii_raw = beta_weighted_integral(a, b, 0.0, 1.0, |t| fail_prob(t).powi(2)); let p_ij_raw = beta_weighted_integral(a, b, 0.0, 1.0, |t| fail_prob(t) * (1.0 - fail_prob(t))); let p_jj_raw = beta_weighted_integral(a, b, 0.0, 1.0, |t| (1.0 - fail_prob(t)).powi(2)); @@ -717,11 +707,11 @@ pub fn livingston_lewis( used_two_parameter, p_tp, p_fp, - p_tn, - p_fn, - accuracy: p_tp + p_tn, - sensitivity: ratio(p_tp, p_tp + p_fn), - specificity: ratio(p_tn, p_tn + p_fp), + p_tf, + p_ff, + accuracy: p_tp + p_tf, + sensitivity: ratio(p_tp, p_tp + p_ff), + specificity: ratio(p_tf, p_tf + p_fp), p_ii, p_ij, p_ji: p_ij, diff --git a/crates/mlsirm-core/src/gtheory.rs b/crates/mlsirm-core/src/gtheory.rs new file mode 100644 index 000000000..6d7f28dbd --- /dev/null +++ b/crates/mlsirm-core/src/gtheory.rs @@ -0,0 +1,419 @@ +//! Generalizability theory G-study / D-study for balanced complete crossed +//! designs: one-facet `p x i` and two-facet `p x i x o`. +//! +//! # Verified sources +//! +//! - **READ IN FULL**: Huebner, A., & Lucht, M. (2019). Generalizability +//! theory in R. *Practical Assessment, Research, and Evaluation*, 24, +//! Article 5. https://doi.org/10.7275/5065-gc10 — provides the model +//! equations (their eqs. 1-2), the D-study component scaling (their +//! Table 1), the relative/absolute error variances sigma^2(delta) and +//! sigma^2(Delta) (their Table 2), the coefficients E-rho^2 and Phi +//! (their eqs. 6-7 with sigma^2(tau) = sigma^2(p) for fully random +//! designs), and complete worked examples with data (their Appendix and +//! Tables 3-6) against which this implementation is tested. +//! - **NOT read**: Brennan (2001) and Shavelson & Webb (1991), the classic +//! texts to which Huebner & Lucht defer the ANOVA/EMS derivations. They +//! are cited below only "as cited in Huebner & Lucht (2019)". +//! +//! # Derivation status of the EMS inversions +//! +//! Huebner & Lucht (2019) do not print the mean-square-to-variance-component +//! formulas. The inversions used here were HAND-DERIVED from the standard +//! expected mean squares of the fully random balanced two-/three-way ANOVA +//! without replication, e.g. for `p x i`: +//! +//! ```text +//! E[MS_p] = sigma2_pi + n_i * sigma2_p +//! E[MS_i] = sigma2_pi + n_p * sigma2_i +//! E[MS_pi] = sigma2_pi +//! ``` +//! +//! and for `p x i x o`: +//! +//! ```text +//! E[MS_p] = s_pio + n_o*s_pi + n_i*s_po + n_i*n_o*s_p (and symmetric) +//! E[MS_pi] = s_pio + n_o*s_pi (and symmetric) +//! E[MS_pio] = s_pio +//! ``` +//! +//! These inversions were verified numerically against every published table +//! of Huebner & Lucht (2019): Table 3/4 (`p x i`, pi_dat) and Table 5/6 +//! (`p x i x o`, pio_cross_dat) reproduce to the papers' printed rounding; +//! see the tests. No claim is made that the derivation matches the unread +//! Brennan (2001) text beyond this numerical agreement. +//! +//! # Formulas (Huebner & Lucht, 2019) +//! +//! D study for proposed sizes `n_i'` (and `n_o'`): each G-study component is +//! divided by the product of the primed sizes of the facets it involves +//! (their Table 1). Then (their Table 2, eqs. 6-7): +//! +//! - `p x i`: `sigma2(delta) = s_pi/n_i'`, +//! `sigma2(Delta) = s_i/n_i' + s_pi/n_i'`. +//! - `p x i x o`: `sigma2(delta) = s_pi/n_i' + s_po/n_o' + s_pio/(n_i' n_o')`, +//! `sigma2(Delta) = sigma2(delta) + s_i/n_i' + s_o/n_o' + s_io/(n_i' n_o')`. +//! - `E-rho^2 = s_p / (s_p + sigma2(delta))`, +//! `Phi = s_p / (s_p + sigma2(Delta))`. +//! +//! # Negative-variance policy (implementation policy, NOT from the paper) +//! +//! ANOVA variance-component estimators can be negative. The paper's examples +//! avoid this by construction. This module reports BOTH the raw ANOVA +//! estimates and component-wise clamped-at-zero copies; all D-study +//! quantities are computed from the CLAMPED components. This "clamped-ANOVA" +//! policy is not a joint constrained estimator; for data with negative raw +//! components the D-study outputs are policy outputs, not reproductions of +//! any published pipeline. `E-rho^2`/`Phi` are `NaN` when their denominator +//! is <= 1e-12. +//! +//! # Scope +//! +//! Balanced complete data only (`NaN` input is an error). Nested designs +//! (`i:p`, `p x (i:o)`), unbalanced data, fixed/mixed facets, and standard +//! errors of the components are out of scope. +//! +//! # References +//! +//! Huebner, A., & Lucht, M. (2019). Generalizability theory in R. +//! *Practical Assessment, Research, and Evaluation*, 24, Article 5. +//! https://doi.org/10.7275/5065-gc10 +//! +//! Brennan, R. L. (2001). *Generalizability theory*. Springer. (As cited in +//! Huebner & Lucht, 2019; not read.) +//! +//! Shavelson, R. J., & Webb, N. M. (1991). *Generalizability theory: A +//! primer*. Sage. (As cited in Huebner & Lucht, 2019; not read.) + +/// One D-study column: proposed facet sizes with the resulting error +/// variances and coefficients (Huebner & Lucht, 2019, Tables 4 and 6). +#[derive(Debug, Clone)] +pub struct GTheoryDStudyRow { + /// Proposed number of items `n_i'`. + pub n_i_prime: usize, + /// Proposed number of occasions `n_o'` (1 for the one-facet design, + /// where it is unused). + pub n_o_prime: usize, + /// Relative error variance `sigma^2(delta)`. + pub rel_error_var: f64, + /// Absolute error variance `sigma^2(Delta)`. + pub abs_error_var: f64, + /// Generalizability coefficient `E-rho^2` (eq. 6). `NaN` when the + /// denominator is <= 1e-12. + pub generalizability: f64, + /// Index of dependability `Phi` (eq. 7). `NaN` when the denominator is + /// <= 1e-12. + pub dependability: f64, +} + +/// G-study + D-study output for the one-facet crossed `p x i` design. +#[derive(Debug, Clone)] +pub struct GTheoryPiResult { + /// Degrees of freedom for (p, i, pi): `n_p-1`, `n_i-1`, + /// `(n_p-1)(n_i-1)`. + pub df: [f64; 3], + /// Sums of squares for (p, i, pi). + pub ss: [f64; 3], + /// Mean squares for (p, i, pi). + pub ms: [f64; 3], + /// Raw ANOVA variance-component estimates for (p, i, pi); may be + /// negative. + pub var_raw: [f64; 3], + /// Component-wise `max(., 0)` of `var_raw`; used for all D-study + /// quantities (clamped-ANOVA policy, see module docs). + pub var: [f64; 3], + /// One row per requested `n_i'`. + pub d_study: Vec, +} + +/// G-study + D-study output for the two-facet crossed `p x i x o` design. +/// Component order everywhere: (p, i, o, pi, po, io, pio). +#[derive(Debug, Clone)] +pub struct GTheoryPioResult { + /// Degrees of freedom in component order. + pub df: [f64; 7], + /// Sums of squares in component order. + pub ss: [f64; 7], + /// Mean squares in component order. + pub ms: [f64; 7], + /// Raw ANOVA variance-component estimates; may be negative. + pub var_raw: [f64; 7], + /// Component-wise `max(., 0)` of `var_raw`; used for all D-study + /// quantities (clamped-ANOVA policy, see module docs). + pub var: [f64; 7], + /// One row per requested `(n_i', n_o')` pair. + pub d_study: Vec, +} + +const DENOM_EPS: f64 = 1e-12; + +fn coef(var_p: f64, err: f64) -> f64 { + let den = var_p + err; + if den <= DENOM_EPS { + f64::NAN + } else { + var_p / den + } +} + +fn validate(x: &[f64], expected_len: usize, what: &str) -> Result<(), String> { + if x.len() != expected_len { + return Err(format!( + "{what}: expected {expected_len} scores, got {}", + x.len() + )); + } + if x.iter().any(|v| !v.is_finite()) { + return Err(format!( + "{what}: non-finite score; complete balanced data required" + )); + } + Ok(()) +} + +/// Generalizability analysis for the one-facet crossed `p x i` design +/// (Huebner & Lucht, 2019, "One-facet p x i design" section). +/// +/// `x` is row-major `n_p x n_i` (person-major); `n_i_prime` lists the +/// proposed D-study item counts. +pub fn gtheory_pi( + x: &[f64], + n_p: usize, + n_i: usize, + n_i_prime: &[usize], +) -> Result { + if n_p < 2 || n_i < 2 { + return Err("gtheory_pi: need at least 2 persons and 2 items".to_string()); + } + validate(x, crate::checked_mul_usize(n_p, n_i, "n_p * n_i overflows usize")?, "gtheory_pi")?; + if n_i_prime.iter().any(|&n| n == 0) { + return Err("gtheory_pi: n_i_prime entries must be >= 1".to_string()); + } + + let (fp, fi) = (n_p as f64, n_i as f64); + let grand = x.iter().sum::() / (fp * fi); + let mut pm = vec![0.0; n_p]; + let mut im = vec![0.0; n_i]; + for p in 0..n_p { + for i in 0..n_i { + pm[p] += x[p * n_i + i]; + im[i] += x[p * n_i + i]; + } + } + for m in pm.iter_mut() { + *m /= fi; + } + for m in im.iter_mut() { + *m /= fp; + } + + let ss_p = fi * pm.iter().map(|m| (m - grand).powi(2)).sum::(); + let ss_i = fp * im.iter().map(|m| (m - grand).powi(2)).sum::(); + let mut ss_pi = 0.0; + for p in 0..n_p { + for i in 0..n_i { + ss_pi += (x[p * n_i + i] - pm[p] - im[i] + grand).powi(2); + } + } + + let df = [fp - 1.0, fi - 1.0, (fp - 1.0) * (fi - 1.0)]; + let ss = [ss_p, ss_i, ss_pi]; + let ms = [ss[0] / df[0], ss[1] / df[1], ss[2] / df[2]]; + // Hand-derived EMS inversion (see module docs); verified against + // Huebner & Lucht (2019) Table 3. + let var_raw = [(ms[0] - ms[2]) / fi, (ms[1] - ms[2]) / fp, ms[2]]; + let var = var_raw.map(|v| v.max(0.0)); + + let d_study = n_i_prime + .iter() + .map(|&n| { + let fni = n as f64; + // Table 1: components divided by n_i'; Table 2: delta = s(pI), + // Delta = s(I) + s(pI). + let rel = var[2] / fni; + let abs = var[1] / fni + var[2] / fni; + GTheoryDStudyRow { + n_i_prime: n, + n_o_prime: 1, + rel_error_var: rel, + abs_error_var: abs, + generalizability: coef(var[0], rel), + dependability: coef(var[0], abs), + } + }) + .collect(); + + Ok(GTheoryPiResult { + df, + ss, + ms, + var_raw, + var, + d_study, + }) +} + +/// Generalizability analysis for the two-facet crossed `p x i x o` design +/// (Huebner & Lucht, 2019, "Two-facet p x i x o design" section). +/// +/// `x` is row-major `x[p*n_i*n_o + i*n_o + o]` (person-major, then item, +/// then occasion). `n_prime` lists the proposed `(n_i', n_o')` pairs. +pub fn gtheory_pio( + x: &[f64], + n_p: usize, + n_i: usize, + n_o: usize, + n_prime: &[(usize, usize)], +) -> Result { + if n_p < 2 || n_i < 2 || n_o < 2 { + return Err("gtheory_pio: need at least 2 levels per facet".to_string()); + } + let n_pi = crate::checked_mul_usize(n_p, n_i, "n_p * n_i overflows usize")?; + validate(x, crate::checked_mul_usize(n_pi, n_o, "n_p * n_i * n_o overflows usize")?, "gtheory_pio")?; + if n_prime.iter().any(|&(a, b)| a == 0 || b == 0) { + return Err("gtheory_pio: n_prime entries must be >= 1".to_string()); + } + + let (fp, fi, fo) = (n_p as f64, n_i as f64, n_o as f64); + let total = fp * fi * fo; + let at = |p: usize, i: usize, o: usize| x[p * n_i * n_o + i * n_o + o]; + let grand = x.iter().sum::() / total; + + let mut mp = vec![0.0; n_p]; + let mut mi = vec![0.0; n_i]; + let mut mo = vec![0.0; n_o]; + let mut mpi = vec![0.0; n_p * n_i]; + let mut mpo = vec![0.0; n_p * n_o]; + let mut mio = vec![0.0; n_i * n_o]; + for p in 0..n_p { + for i in 0..n_i { + for o in 0..n_o { + let v = at(p, i, o); + mp[p] += v; + mi[i] += v; + mo[o] += v; + mpi[p * n_i + i] += v; + mpo[p * n_o + o] += v; + mio[i * n_o + o] += v; + } + } + } + for m in mp.iter_mut() { + *m /= fi * fo; + } + for m in mi.iter_mut() { + *m /= fp * fo; + } + for m in mo.iter_mut() { + *m /= fp * fi; + } + for m in mpi.iter_mut() { + *m /= fo; + } + for m in mpo.iter_mut() { + *m /= fi; + } + for m in mio.iter_mut() { + *m /= fp; + } + + let ss_p = fi * fo * mp.iter().map(|m| (m - grand).powi(2)).sum::(); + let ss_i = fp * fo * mi.iter().map(|m| (m - grand).powi(2)).sum::(); + let ss_o = fp * fi * mo.iter().map(|m| (m - grand).powi(2)).sum::(); + let mut ss_pi = 0.0; + for p in 0..n_p { + for i in 0..n_i { + ss_pi += (mpi[p * n_i + i] - mp[p] - mi[i] + grand).powi(2); + } + } + ss_pi *= fo; + let mut ss_po = 0.0; + for p in 0..n_p { + for o in 0..n_o { + ss_po += (mpo[p * n_o + o] - mp[p] - mo[o] + grand).powi(2); + } + } + ss_po *= fi; + let mut ss_io = 0.0; + for i in 0..n_i { + for o in 0..n_o { + ss_io += (mio[i * n_o + o] - mi[i] - mo[o] + grand).powi(2); + } + } + ss_io *= fp; + // Direct three-way residual sum (algebraically equal to + // SS_total - all other SS, but immune to the catastrophic cancellation + // the subtraction form shows for large-offset data). + let mut ss_pio = 0.0; + for p in 0..n_p { + for i in 0..n_i { + for o in 0..n_o { + ss_pio += (at(p, i, o) - mpi[p * n_i + i] - mpo[p * n_o + o] - mio[i * n_o + o] + + mp[p] + + mi[i] + + mo[o] + - grand) + .powi(2); + } + } + } + + let df = [ + fp - 1.0, + fi - 1.0, + fo - 1.0, + (fp - 1.0) * (fi - 1.0), + (fp - 1.0) * (fo - 1.0), + (fi - 1.0) * (fo - 1.0), + (fp - 1.0) * (fi - 1.0) * (fo - 1.0), + ]; + let ss = [ss_p, ss_i, ss_o, ss_pi, ss_po, ss_io, ss_pio]; + let mut ms = [0.0; 7]; + for k in 0..7 { + ms[k] = ss[k] / df[k]; + } + // Hand-derived EMS inversion (see module docs); verified against + // Huebner & Lucht (2019) Table 5. + let var_raw = [ + (ms[0] - ms[3] - ms[4] + ms[6]) / (fi * fo), + (ms[1] - ms[3] - ms[5] + ms[6]) / (fp * fo), + (ms[2] - ms[4] - ms[5] + ms[6]) / (fp * fi), + (ms[3] - ms[6]) / fo, + (ms[4] - ms[6]) / fi, + (ms[5] - ms[6]) / fp, + ms[6], + ]; + let var = var_raw.map(|v| v.max(0.0)); + + let d_study = n_prime + .iter() + .map(|&(ni, no)| { + let (fni, fno) = (ni as f64, no as f64); + // Table 2: delta = s(pI) + s(pO) + s(pIO); + // Delta = delta + s(I) + s(O) + s(IO). + let rel = var[3] / fni + var[4] / fno + var[6] / (fni * fno); + let abs = rel + var[1] / fni + var[2] / fno + var[5] / (fni * fno); + GTheoryDStudyRow { + n_i_prime: ni, + n_o_prime: no, + rel_error_var: rel, + abs_error_var: abs, + generalizability: coef(var[0], rel), + dependability: coef(var[0], abs), + } + }) + .collect(); + + Ok(GTheoryPioResult { + df, + ss, + ms, + var_raw, + var, + d_study, + }) +} + +#[cfg(test)] +#[path = "../../../tests/unit/gtheory_tests.rs"] +mod tests; diff --git a/crates/mlsirm-core/src/lib.rs b/crates/mlsirm-core/src/lib.rs index 9a6325444..f833d1e33 100644 --- a/crates/mlsirm-core/src/lib.rs +++ b/crates/mlsirm-core/src/lib.rs @@ -9,6 +9,7 @@ pub mod facets; pub mod fitstats; pub mod gpcm; pub mod grm; +pub mod gtheory; pub mod ksirt; pub mod linking; pub mod lltm; diff --git a/python/fast_mlsirm/__init__.py b/python/fast_mlsirm/__init__.py index 2dea90e1b..d7c6d84a6 100644 --- a/python/fast_mlsirm/__init__.py +++ b/python/fast_mlsirm/__init__.py @@ -60,6 +60,12 @@ feldt_alpha_ci as feldt_alpha_ci, AlphaCiResult as AlphaCiResult, ) +from .gtheory import ( + gtheory_pi as gtheory_pi, + gtheory_pio as gtheory_pio, + GTheoryResult as GTheoryResult, + GTheoryDStudyRow as GTheoryDStudyRow, +) from .mokken import mokken_analysis as mokken_analysis, MokkenResult as MokkenResult from .rsm import fit_rsm as fit_rsm, RsmFit as RsmFit from .mixed import fit_mixed_items as fit_mixed_items, MixedFormatFit as MixedFormatFit, MixedItemParameters as MixedItemParameters @@ -192,6 +198,10 @@ "cronbach_alpha", "feldt_alpha_ci", "AlphaCiResult", + "gtheory_pi", + "gtheory_pio", + "GTheoryResult", + "GTheoryDStudyRow", "ClassificationResult", "DetectResult", "fit_mixed_items", diff --git a/python/fast_mlsirm/classification.py b/python/fast_mlsirm/classification.py index 86c8f20c0..740720561 100644 --- a/python/fast_mlsirm/classification.py +++ b/python/fast_mlsirm/classification.py @@ -35,23 +35,13 @@ class ClassificationResult: class LivingstonLewisResult: """Livingston-Lewis single-administration classification results. - Accuracy cells — joint proportions over (true state, observed state), - where "pass" means observed score at or above the cut: - - * ``p_tp``: P(true pass, observed pass) — true positive. - * ``p_fp``: P(true fail, observed pass) — false positive. - * ``p_tn``: P(true fail, observed fail) — true negative. - * ``p_fn``: P(true pass, observed fail) — false negative. - - The four cells sum to 1. ``accuracy = p_tp + p_tn``; - ``sensitivity = p_tp / (p_tp + p_fn)`` (P(pass | true pass)); - ``specificity = p_tn / (p_tn + p_fp)`` (P(fail | true fail)). - - Consistency cells — joint proportions over two hypothetical parallel - forms (``i`` = fail, ``j`` = pass), with ``p_ij == p_ji`` by - construction (single rounded threshold in both blocks — betafunctions' - round/floor mix makes its cells asymmetric; divergence documented in - the Rust core).""" + ``p_tp``/``p_fp``/``p_tf``/``p_ff`` are the accuracy cells (pass = + observed score at or above the cut; ``t``/``f`` = true pass/fail); + ``p_ii``/``p_ij``/``p_ji``/``p_jj`` are consistency cells over two + hypothetical parallel forms, with ``p_ij == p_ji`` by construction + (single rounded threshold in both blocks — betafunctions' round/floor + mix makes its cells asymmetric; divergence documented in the Rust + core).""" effective_test_length: float etl_rounded: int @@ -62,8 +52,8 @@ class LivingstonLewisResult: used_two_parameter: bool p_tp: float p_fp: float - p_tn: float - p_fn: float + p_tf: float + p_ff: float accuracy: float sensitivity: float specificity: float @@ -253,8 +243,8 @@ def livingston_lewis( used_two_parameter=bool(res["used_two_parameter"]), p_tp=float(res["p_tp"]), p_fp=float(res["p_fp"]), - p_tn=float(res["p_tn"]), - p_fn=float(res["p_fn"]), + p_tf=float(res["p_tf"]), + p_ff=float(res["p_ff"]), accuracy=float(res["accuracy"]), sensitivity=float(res["sensitivity"]), specificity=float(res["specificity"]), diff --git a/python/fast_mlsirm/gtheory.py b/python/fast_mlsirm/gtheory.py new file mode 100644 index 000000000..407c8da3a --- /dev/null +++ b/python/fast_mlsirm/gtheory.py @@ -0,0 +1,144 @@ +"""Generalizability theory G/D-study analyses for crossed designs +(Huebner & Lucht, 2019). All numeric work happens in the Rust core +(``mlsirm_core::gtheory``); this module only validates and marshals. + +Source status: Huebner & Lucht (2019) READ in full, including the worked +p x i and p x i x o examples (Tables 3-6) that the Rust tests reproduce. +Brennan (2001) and Shavelson & Webb (1991) are cited by that paper for the +EMS derivations and were NOT read; the EMS-to-variance-component inversions +are hand-derived and numerically verified against the paper's published +tables (see the Rust module docs). +""" + +from __future__ import annotations + +from dataclasses import dataclass +from typing import Sequence + +import numpy as np + +_REFERENCES = """References (APA 7th ed.): + Huebner, A., & Lucht, M. (2019). Generalizability theory in R. + *Practical Assessment, Research, and Evaluation, 24*, Article 5. + https://doi.org/10.7275/5065-gc10 + Brennan, R. L. (2001). *Generalizability theory*. Springer. + (As cited in Huebner & Lucht, 2019; not read.) + Shavelson, R. J., & Webb, N. M. (1991). *Generalizability theory: + A primer*. Sage. (As cited in Huebner & Lucht, 2019; not read.) + """ + + +@dataclass +class GTheoryDStudyRow: + """One D-study column: proposed facet sizes with the resulting error + variances and coefficients (Huebner & Lucht, 2019, Tables 4 and 6). + ``n_o_prime`` is 1 and unused for the one-facet design.""" + + n_i_prime: int + n_o_prime: int + rel_error_var: float + abs_error_var: float + generalizability: float + dependability: float + + +@dataclass +class GTheoryResult: + """G-study ANOVA table plus D-study rows. + + Component order is ``(p, i, pi)`` for the one-facet design and + ``(p, i, o, pi, po, io, pio)`` for the two-facet design. ``var_raw`` + holds the raw ANOVA estimates (may be negative); ``var`` is the + component-wise ``max(., 0)`` used for all D-study quantities + (clamped-ANOVA policy — an implementation choice, not a + paper-prescribed estimator).""" + + df: list[float] + ss: list[float] + ms: list[float] + var_raw: list[float] + var: list[float] + d_study: list[GTheoryDStudyRow] + + +def _to_result(res: dict) -> GTheoryResult: + return GTheoryResult( + df=[float(v) for v in res["df"]], + ss=[float(v) for v in res["ss"]], + ms=[float(v) for v in res["ms"]], + var_raw=[float(v) for v in res["var_raw"]], + var=[float(v) for v in res["var"]], + d_study=[ + GTheoryDStudyRow( + n_i_prime=int(r["n_i_prime"]), + n_o_prime=int(r["n_o_prime"]), + rel_error_var=float(r["rel_error_var"]), + abs_error_var=float(r["abs_error_var"]), + generalizability=float(r["generalizability"]), + dependability=float(r["dependability"]), + ) + for r in res["d_study"] + ], + ) + + +def _core_or_raise(name: str): + from .fitstats import _core_module + + core = _core_module() + if core is None or not hasattr(core, name): + raise RuntimeError(f"{name} requires the compiled Rust core") + return core + + +def gtheory_pi( + data: np.ndarray, + n_i_prime: Sequence[int] = (5, 10, 15, 20), +) -> GTheoryResult: + """One-facet crossed ``p x i`` generalizability analysis (compute in + Rust; Huebner & Lucht, 2019, "One-facet crossed design" section and + Tables 3-4). + + ``data`` is a complete, balanced ``n_persons x n_items`` score array. + Variance components are the ANOVA (EMS) estimators; negative raw + estimates are reported in ``var_raw`` and clamped to zero in ``var`` + for the D study. D-study rows give sigma^2(delta), sigma^2(Delta), the + generalizability coefficient E-rho^2 (eq. 6) and the dependability + index Phi (eq. 7) at each proposed ``n_i'``; coefficients are NaN when + their denominator is <= 1e-12. In LLM-as-a-Judge quality management + this asks how many judge items are needed for a dependable rating. + + """ + _REFERENCES + core = _core_or_raise("gtheory_pi") + x = np.ascontiguousarray(np.asarray(data, dtype=np.float64)) + if x.ndim != 2: + raise ValueError("data must be a 2-D persons x items array") + n_p, n_i = x.shape + primes = [int(v) for v in n_i_prime] + return _to_result(core.gtheory_pi(x.reshape(-1), int(n_p), int(n_i), primes)) + + +def gtheory_pio( + data: np.ndarray, + n_prime: Sequence[tuple[int, int]] = ((5, 2), (10, 2), (15, 2), (20, 2)), +) -> GTheoryResult: + """Two-facet crossed ``p x i x o`` generalizability analysis (compute + in Rust; Huebner & Lucht, 2019, "Two-facet crossed design" section and + Tables 5-6). + + ``data`` is a complete, balanced ``n_persons x n_items x n_occasions`` + score array. Component order everywhere is + ``(p, i, o, pi, po, io, pio)``. ``n_prime`` lists proposed + ``(n_i', n_o')`` D-study pairs; the clamped-ANOVA and NaN-denominator + policies match :func:`gtheory_pi`. + + """ + _REFERENCES + core = _core_or_raise("gtheory_pio") + x = np.ascontiguousarray(np.asarray(data, dtype=np.float64)) + if x.ndim != 3: + raise ValueError("data must be a 3-D persons x items x occasions array") + n_p, n_i, n_o = x.shape + pairs = [(int(a), int(b)) for a, b in n_prime] + return _to_result( + core.gtheory_pio(x.reshape(-1), int(n_p), int(n_i), int(n_o), pairs) + ) diff --git a/tests/test_paper_features.py b/tests/test_paper_features.py index fee82886a..2a732e2b1 100644 --- a/tests/test_paper_features.py +++ b/tests/test_paper_features.py @@ -5752,8 +5752,8 @@ def test_livingston_lewis_matches_independent_reference(): tol = 1e-7 assert abs(r.p_tp - 0.40501323923680005) < tol assert abs(r.p_fp - 0.07917200341607164) < tol - assert abs(r.p_tn - 0.46137445417258505) < tol - assert abs(r.p_fn - 0.05444030317454212) < tol + assert abs(r.p_tf - 0.46137445417258505) < tol + assert abs(r.p_ff - 0.05444030317454212) < tol assert abs(r.accuracy - 0.866387693409385) < tol assert abs(r.sensitivity - 0.8815107553881857) < tol assert abs(r.specificity - 0.8535333969826517) < tol @@ -5784,4 +5784,84 @@ def test_livingston_lewis_rejects_degenerate_inputs(): with pytest.raises(ValueError): livingston_lewis(good, 0.82, 12.0, 50.0, 28.0) # score below min with pytest.raises(ValueError): - livingston_lewis([20.0, 20.0, 20.0], 0.82, 0.0, 50.0, 10.0) # var 0 \ No newline at end of file + livingston_lewis([20.0, 20.0, 20.0], 0.82, 0.0, 50.0, 10.0) # var 0 + +def test_gtheory_pi_matches_paper_tables(): + """Huebner & Lucht (2019) appendix pi_dat; asserts read the crate's + var/var_raw/d_study outputs (Tables 3-4 reproduced by the Rust tests + to full precision; pinned literals from the same verified pipeline).""" + from fast_mlsirm import gtheory_pi + from fast_mlsirm.fitstats import _core_module + + core = _core_module() + if core is None or not hasattr(core, "gtheory_pi"): + pytest.skip("compiled core built without gtheory_pi") + + x = np.array( + [ + [9.0, 9.0, 7.0, 4.0], + [9.0, 8.0, 4.0, 6.0], + [8.0, 8.0, 6.0, 2.0], + [9.0, 8.0, 6.0, 3.0], + [10.0, 9.0, 8.0, 7.0], + [6.0, 4.0, 5.0, 1.0], + ] + ) + r = gtheory_pi(x, n_i_prime=[4, 8]) + tol = 1e-9 + assert abs(r.var[0] - 1.905555555556) < tol + assert abs(r.var[1] - 4.027777777778) < tol + assert abs(r.var[2] - 1.277777777778) < tol + assert r.var == r.var_raw # all raw components positive here + d4 = r.d_study[0] + assert d4.n_i_prime == 4 + assert abs(d4.rel_error_var - 1.277777777778 / 4.0) < tol + assert abs(d4.generalizability - 0.856429463171) < tol + assert abs(d4.dependability - 0.589600343790) < tol + d8 = r.d_study[1] + # sigma^2(p) invariant across n' (weak identity); delta halves. + assert abs(d8.rel_error_var - d4.rel_error_var / 2.0) < tol + assert d8.generalizability > d4.generalizability + + +def test_gtheory_pio_matches_paper_tables_and_rejects_bad_input(): + """Huebner & Lucht (2019) appendix pio_cross_dat (6 x 4 x 2); asserts + read the crate's 7-component var vector (order p,i,o,pi,po,io,pio) and + D-study coefficients pinned from the verified pipeline (Tables 5-6).""" + from fast_mlsirm import gtheory_pi, gtheory_pio + from fast_mlsirm.fitstats import _core_module + + core = _core_module() + if core is None or not hasattr(core, "gtheory_pio"): + pytest.skip("compiled core built without gtheory_pio") + + # pio_cross_dat in crate layout x[p*n_i*n_o + i*n_o + o], i.e. a + # (6, 4, 2) persons x items x occasions array. + scores = np.array( + [ + 9.0, 9.0, 9.0, 8.0, 7.0, 5.0, 4.0, 5.0, + 9.0, 6.0, 8.0, 5.0, 4.0, 3.0, 6.0, 3.0, + 8.0, 8.0, 8.0, 7.0, 6.0, 3.0, 2.0, 2.0, + 9.0, 9.0, 8.0, 6.0, 6.0, 6.0, 3.0, 2.0, + 10.0, 8.0, 9.0, 8.0, 8.0, 9.0, 7.0, 7.0, + 6.0, 3.0, 4.0, 2.0, 5.0, 3.0, 1.0, 2.0, + ] + ).reshape(6, 4, 2) + r = gtheory_pio(scores, n_prime=[(4, 2)]) + tol = 1e-9 + want = [ + 2.483333333333, 3.066666666667, 0.577777777778, 0.766666666667, + 0.088888888889, 0.022222222222, 0.831944444444, + ] + for got, w in zip(r.var, want): + assert abs(got - w) < tol + d = r.d_study[0] + assert d.n_i_prime == 4 and d.n_o_prime == 2 + # delta = pi/4 + po/2 + pio/8 (kills a dropped-term mutation). + assert abs(d.rel_error_var - (want[3] / 4.0 + want[4] / 2.0 + want[6] / 8.0)) < tol + assert 0.0 < d.dependability < d.generalizability < 1.0 + + with pytest.raises(ValueError): + gtheory_pio(scores[:1]) # n_p < 2 + with pytest.raises(ValueError): + gtheory_pi(np.array([[1.0, np.nan], [2.0, 3.0]])) # non-finite diff --git a/tests/unit/classification_tests.rs b/tests/unit/classification_tests.rs index 1f181e739..3d7eb1b6d 100644 --- a/tests/unit/classification_tests.rs +++ b/tests/unit/classification_tests.rs @@ -214,7 +214,7 @@ fn monte_carlo_accuracy_orders_test_quality() { // // Mutation-kill map (spot-checked by actually applying each mutation): // - M1 swap upper/lower binomial tail in accuracy integrands -> -// p_tp/p_fn literals FAIL. +// p_tp/p_ff literals FAIL. // - M2 use unrounded ETL in the passing threshold -> k shifts, literals FAIL // (fixture A has ETL = 91.123, materially non-integer). // - M4 threshold off-by-one (k-1 -> k) -> literals FAIL (fixture A has @@ -321,8 +321,8 @@ fn ll_matches_independent_reference_two_parameter() { let tol = 1e-7; assert!((r.p_tp - 0.23854946889141398).abs() < tol); assert!((r.p_fp - 0.06401868217743813).abs() < tol); - assert!((r.p_tn - 0.6521887072379015).abs() < tol); - assert!((r.p_fn - 0.04524314169324793).abs() < tol); + assert!((r.p_tf - 0.6521887072379015).abs() < tol); + assert!((r.p_ff - 0.04524314169324793).abs() < tol); assert!((r.accuracy - 0.8907381761293154).abs() < tol); assert!((r.sensitivity - 0.8405767451096093).abs() < tol); assert!((r.specificity - 0.9106143232762532).abs() < tol); @@ -345,7 +345,7 @@ fn ll_extreme_cut_distinguishes_rounded_etl_in_k() { assert_eq!(r.etl_rounded, 91); let tol = 1e-8; assert!((r.p_fp - 2.152128586520244e-06).abs() < tol); - assert!((r.p_tn - 0.9999978476968782).abs() < 1e-7); + assert!((r.p_tf - 0.9999978476968782).abs() < 1e-7); assert!((r.kappa - 0.04370046162326584).abs() < 1e-7); } @@ -362,8 +362,8 @@ fn ll_failsafe_engages_on_skewed_data() { let tol = 1e-7; assert!((r.p_tp - 0.8432448010046247).abs() < tol); assert!((r.p_fp - 0.025111134448713338).abs() < tol); - assert!((r.p_tn - 0.0853397836347724).abs() < tol); - assert!((r.p_fn - 0.046304280911892266).abs() < tol); + assert!((r.p_tf - 0.0853397836347724).abs() < tol); + assert!((r.p_ff - 0.046304280911892266).abs() < tol); assert!((r.accuracy - 0.9285845846393972).abs() < tol); assert!((r.consistency - 0.89849902437627).abs() < tol); assert!((r.kappa - 0.5560427413146216).abs() < tol); @@ -387,8 +387,8 @@ fn ll_four_parameter_path() { let tol = 1e-7; assert!((r.p_tp - 0.5543443184884189).abs() < tol); assert!((r.p_fp - 0.044475496859695596).abs() < tol); - assert!((r.p_tn - 0.3583220858730486).abs() < tol); - assert!((r.p_fn - 0.04285809877883742).abs() < tol); + assert!((r.p_tf - 0.3583220858730486).abs() < tol); + assert!((r.p_ff - 0.04285809877883742).abs() < tol); assert!((r.accuracy - 0.9126664043614675).abs() < tol); assert!((r.sensitivity - 0.928235222196601).abs() < tol); assert!((r.specificity - 0.8895835060430215).abs() < tol); @@ -407,10 +407,10 @@ fn ll_structural_invariants_read_crate_fields() { // desynchronizes them (e.g. computing kappa from raw unnormalized cells) // fails here even if it preserved the individual pinned literals. let r = fixture_a(); - assert!((r.p_tp + r.p_fp + r.p_tn + r.p_fn - 1.0).abs() < 1e-9); + assert!((r.p_tp + r.p_fp + r.p_tf + r.p_ff - 1.0).abs() < 1e-9); assert!((r.p_ii + r.p_ij + r.p_ji + r.p_jj - 1.0).abs() < 1e-12); assert_eq!(r.p_ij, r.p_ji); // by construction (disclosed above) - assert!((r.accuracy - (r.p_tp + r.p_tn)).abs() < 1e-15); + assert!((r.accuracy - (r.p_tp + r.p_tf)).abs() < 1e-15); assert!((r.consistency - (r.p_ii + r.p_jj)).abs() < 1e-15); let pc = (r.p_ii + r.p_ij) * (r.p_ii + r.p_ji) + (r.p_ij + r.p_jj) * (r.p_ji + r.p_jj); assert!((r.chance_consistency - pc).abs() < 1e-15); @@ -508,7 +508,7 @@ fn ll_conditional_ratios_nan_when_margin_vanishes() { let r = livingston_lewis(&scores, 0.9, 0.0, 50.0, 5.0).unwrap(); assert!(r.specificity.is_nan()); assert!(r.sensitivity.is_finite()); - assert!((r.p_tp + r.p_fp + r.p_tn + r.p_fn - 1.0).abs() < 1e-9); + assert!((r.p_tp + r.p_fp + r.p_tf + r.p_ff - 1.0).abs() < 1e-9); let scores_a = gen_scores(250, 60, 42, 0.15, 0.95, 4); let r = livingston_lewis(&scores_a, 0.85, 0.0, 60.0, 0.1).unwrap(); diff --git a/tests/unit/gtheory_tests.rs b/tests/unit/gtheory_tests.rs new file mode 100644 index 000000000..0885ecab5 --- /dev/null +++ b/tests/unit/gtheory_tests.rs @@ -0,0 +1,381 @@ +//! Tests for the generalizability-theory module (`gtheory.rs`). +//! +//! Every assertion reads values RETURNED BY THE CRATE (`gtheory_pi` / +//! `gtheory_pio` outputs); no assert recomputes the formula locally. +//! +//! Fixture provenance: +//! - `pi_dat` and `pio_cross_dat`: the data appendix of Huebner & Lucht +//! (2019), with expected values pinned from an independent NumPy +//! implementation of the same ANOVA/EMS pipeline, which reproduces the +//! paper's published Tables 3-6 to their printed rounding. +//! - Independent fixtures: generated with numpy `default_rng(20260723)`; +//! expected literals pinned from the same independent NumPy pipeline. +//! +//! Mutation kills. M1, M3, M5, M6 were each ACTUALLY applied to +//! src/gtheory.rs, observed to FAIL the named tests, and reverted +//! (2026-07-23). M2 and M4 are documented expectations only (same code +//! sites as M1/M3, killed by the same fixtures; not separately executed): +//! - M1 (EXECUTED) swap person/item EMS divisors in `gtheory_pi` +//! (`/fi` <-> `/fp`): `gt_pi_paper_fixture`, `gt_pi_independent_fixture`, +//! and `gt_pi_negative_variance_clamped` FAIL. +//! - M2 (documented) drop the `- ms_pi` subtraction in `gtheory_pi` var_p: +//! var_p becomes 2.225 -> `gt_pi_paper_fixture` fails. +//! - M3 (EXECUTED) flip the `+ ms[6]` sign in the `gtheory_pio` var_p +//! inversion: `gt_pio_paper_fixture` and `gt_pio_independent_fixture` FAIL. +//! - M4 (documented) drop the `s_io/(n_i' n_o')` term from the absolute +//! error: Delta(1,1) becomes 5.3319 vs 5.3542 -> `gt_pio_paper_fixture` +//! fails. +//! - M5 (EXECUTED) drop the clamp (`var = var_raw`) in `gtheory_pi`: +//! `gt_pi_negative_variance_clamped` FAILS (var == var_raw identity and +//! sign-sensitive D-study outputs). +//! - M6 (EXECUTED) transpose the pio layout (swap `i`/`o` strides): the +//! 5x4x3 independent fixture has n_i=4 != n_o=3 -> +//! `gt_pio_independent_fixture` and `gt_pio_paper_fixture` FAIL. +//! +//! Weak identity checks (disclosed, NOT counted as kills): `n' = 1` D-study +//! columns (divisor identity point), DF values, and sigma^2(p) invariance +//! across n'. + +use super::{gtheory_pi, gtheory_pio}; + +const TOL: f64 = 1e-9; + +fn assert_close(got: f64, want: f64, tol: f64, what: &str) { + assert!((got - want).abs() <= tol, "{what}: got {got}, want {want}"); +} + +/// Huebner & Lucht (2019) appendix `pi_dat` (6 persons x 4 items). +fn pi_dat() -> Vec { + vec![ + 9., 9., 7., 4., 9., 8., 4., 6., 8., 8., 6., 2., 9., 8., 6., 3., 10., 9., 8., 7., 6., 4., + 5., 1., + ] +} + +/// Huebner & Lucht (2019) appendix `pio_cross_dat` (6 x 4 x 2), re-laid out +/// to the crate layout x[p*n_i*n_o + i*n_o + o]. +fn pio_cross_dat() -> Vec { + vec![ + 9., 9., 9., 8., 7., 5., 4., 5., 9., 6., 8., 5., 4., 3., 6., 3., 8., 8., 8., 7., 6., 3., 2., + 2., 9., 9., 8., 6., 6., 6., 3., 2., 10., 8., 9., 8., 8., 9., 7., 7., 6., 3., 4., 2., 5., + 3., 1., 2., + ] +} + +/// Paper fixture: reproduces Tables 3 and 4 of Huebner & Lucht (2019). +/// Asserts read: ss, ms, var_raw, var, d_study rows from the crate result. +/// Kills M1, M2 (see header). +#[test] +fn gt_pi_paper_fixture() { + let r = gtheory_pi(&pi_dat(), 6, 4, &[4, 10]).unwrap(); + // ANOVA table (paper Table 3, pinned at full precision). + assert_eq!(r.df, [5.0, 3.0, 15.0]); + assert_close(r.ss[0], 44.5, TOL, "ss_p"); + assert_close(r.ss[1], 76.333333333333, 1e-9, "ss_i"); + assert_close(r.ss[2], 19.166666666667, 1e-9, "ss_pi"); + assert_close(r.ms[0], 8.9, TOL, "ms_p"); + assert_close(r.ms[1], 25.444444444444, 1e-9, "ms_i"); + assert_close(r.ms[2], 1.277777777778, 1e-9, "ms_pi"); + // Variance components (paper: 1.91 / 4.03 / 1.28). + assert_close(r.var_raw[0], 1.905555555556, 1e-9, "var_p"); + assert_close(r.var_raw[1], 4.027777777778, 1e-9, "var_i"); + assert_close(r.var_raw[2], 1.277777777778, 1e-9, "var_pi"); + // All components positive here: clamped == raw. + assert_eq!(r.var, r.var_raw); + // D study (paper Table 4 columns n_i' = 4 and 10). + let d4 = &r.d_study[0]; + assert_eq!((d4.n_i_prime, d4.n_o_prime), (4, 1)); + assert_close(d4.rel_error_var, 0.319444444444, 1e-9, "delta(4)"); + assert_close(d4.abs_error_var, 1.326388888889, 1e-9, "Delta(4)"); + assert_close(d4.generalizability, 0.856429463171, 1e-9, "Erho2(4)"); + assert_close(d4.dependability, 0.589600343790, 1e-9, "Phi(4)"); + let d10 = &r.d_study[1]; + assert_close(d10.rel_error_var, 0.127777777778, 1e-9, "delta(10)"); + assert_close(d10.abs_error_var, 0.530555555556, 1e-9, "Delta(10)"); + assert_close(d10.generalizability, 0.937158469945, 1e-9, "Erho2(10)"); + assert_close(d10.dependability, 0.782212086659, 1e-9, "Phi(10)"); +} + +/// Paper fixture: reproduces Tables 5 and 6 of Huebner & Lucht (2019). +/// Asserts read: ss, ms, var_raw, d_study rows from the crate result. +/// Kills M3, M4 (see header). +#[test] +fn gt_pio_paper_fixture() { + let r = gtheory_pio(&pio_cross_dat(), 6, 4, 2, &[(4, 2), (1, 1), (5, 3)]).unwrap(); + assert_eq!(r.df, [5.0, 3.0, 1.0, 15.0, 5.0, 3.0, 15.0]); + let want_ss = [ + 112.9375, + 117.895833333333, + 15.1875, + 35.479166666667, + 5.9375, + 2.895833333333, + 12.479166666667, + ]; + let want_ms = [ + 22.5875, + 39.298611111111, + 15.1875, + 2.365277777778, + 1.1875, + 0.965277777778, + 0.831944444444, + ]; + // Paper Table 5: 2.48 / 3.07 / 0.58 / 0.77 / 0.09 / 0.02 / 0.83. + let want_var = [ + 2.483333333333, + 3.066666666667, + 0.577777777778, + 0.766666666667, + 0.088888888889, + 0.022222222222, + 0.831944444444, + ]; + for k in 0..7 { + assert_close(r.ss[k], want_ss[k], 1e-9, "pio ss"); + assert_close(r.ms[k], want_ms[k], 1e-9, "pio ms"); + assert_close(r.var_raw[k], want_var[k], 1e-9, "pio var"); + } + assert_eq!(r.var, r.var_raw); + // Paper Table 6 columns (n_i', n_o') = (4,2), (1,1), (5,3). + let want_d = [ + // rel, abs, Erho2, Phi + [0.340104166667, 1.3984375, 0.879542519830, 0.639742385617], + [1.6875, 5.354166666667, 0.595404595405, 0.316852737905], + [ + 0.238425925926, + 1.045833333333, + 0.912400068039, + 0.703659976387, + ], + ]; + for (row, want) in r.d_study.iter().zip(want_d.iter()) { + assert_close(row.rel_error_var, want[0], 1e-9, "pio delta"); + assert_close(row.abs_error_var, want[1], 1e-9, "pio Delta"); + assert_close(row.generalizability, want[2], 1e-9, "pio Erho2"); + assert_close(row.dependability, want[3], 1e-9, "pio Phi"); + } + assert_eq!((r.d_study[2].n_i_prime, r.d_study[2].n_o_prime), (5, 3)); +} + +/// Independent 7x5 fixture (numpy default_rng(20260723); literals pinned +/// from an independent NumPy pipeline). n_p=7 != n_i=5 keeps the divisor +/// asymmetry live at a size not used by the paper fixture. +#[test] +fn gt_pi_independent_fixture() { + #[rustfmt::skip] + let x = vec![ + 50.1, 39.0, 41.1, 56.9, 52.3, + 52.4, 52.5, 41.8, 46.3, 37.6, + 59.9, 43.2, 52.7, 47.0, 54.5, + 34.1, 52.4, 45.1, 54.1, 32.9, + 59.7, 46.7, 57.7, 57.9, 45.5, + 54.5, 47.6, 61.2, 55.2, 45.6, + 49.8, 32.8, 69.3, 60.1, 57.7, + ]; + let r = gtheory_pi(&x, 7, 5, &[8]).unwrap(); + assert_close(r.ss[0], 484.0, 1e-8, "ss_p"); + assert_close(r.ss[1], 439.2845714286, 1e-8, "ss_i"); + assert_close(r.ss[2], 1623.3514285714, 1e-8, "ss_pi"); + assert_close(r.ms[0], 80.6666666667, 1e-8, "ms_p"); + assert_close(r.ms[1], 109.8211428571, 1e-8, "ms_i"); + assert_close(r.ms[2], 67.6396428571, 1e-8, "ms_pi"); + assert_close(r.var_raw[0], 2.6054047619, 1e-8, "var_p"); + assert_close(r.var_raw[1], 6.0259285714, 1e-8, "var_i"); + let d = &r.d_study[0]; + assert_close(d.rel_error_var, 8.4549553571, 1e-8, "delta(8)"); + assert_close(d.abs_error_var, 9.2081964286, 1e-8, "Delta(8)"); + assert_close(d.generalizability, 0.2355623808, 1e-8, "Erho2(8)"); + assert_close(d.dependability, 0.2205428065, 1e-8, "Phi(8)"); +} + +/// Negative-variance anchor for the clamped-ANOVA policy (module docs): +/// item means are nearly equal but the interaction is large, so the raw +/// sigma^2_i and sigma^2_p are negative. Asserts read var_raw, var, and the +/// D-study row from the crate result. Kills M5. +#[test] +fn gt_pi_negative_variance_clamped() { + let x = vec![5., 1., 4., 2., 6., 3., 7., 3., 8., 1., 5., 2.]; + let r = gtheory_pi(&x, 4, 3, &[3]).unwrap(); + assert_close(r.var_raw[0], -0.083333333333, 1e-9, "raw var_p < 0"); + assert_close(r.var_raw[1], -1.555555555556, 1e-9, "raw var_i < 0"); + assert!(r.var_raw[0] < 0.0 && r.var_raw[1] < 0.0); + assert_eq!(r.var[0], 0.0); + assert_eq!(r.var[1], 0.0); + assert_close(r.var[2], 6.555555555556, 1e-9, "var_pi"); + let d = &r.d_study[0]; + // Clamped policy: delta = Delta = var_pi/3 (clamped var_i contributes 0). + // Clamped var_p = 0 with a POSITIVE error variance gives Erho2 = Phi = 0 + // (zero universe-score variance), not NaN — the NaN rule only fires when + // the whole denominator is <= 1e-12. + assert_close(d.rel_error_var, 6.555555555556 / 3.0, 1e-9, "delta"); + assert_close(d.abs_error_var, 6.555555555556 / 3.0, 1e-9, "Delta"); + assert_eq!(d.generalizability, 0.0, "Erho2 = 0 at zero var_p"); + assert_eq!(d.dependability, 0.0, "Phi = 0 at zero var_p"); +} + +/// Independent 5x4x3 fixture (numpy default_rng(20260723); literals pinned +/// from an independent NumPy pipeline). All facet sizes distinct +/// (n_p=5, n_i=4, n_o=3) so any axis/stride mix-up shifts every output. +/// Raw sigma^2_pi is negative here, so this doubles as the pio clamp +/// anchor. Kills M5, M6. +#[test] +fn gt_pio_independent_fixture() { + #[rustfmt::skip] + let x = vec![ + 19.3, 17.2, 16.0, 15.9, 18.0, 17.8, 21.4, 22.4, 12.0, 18.1, 17.1, 14.3, + 19.8, 19.6, 18.5, 21.4, 21.4, 17.8, 23.3, 26.0, 17.8, 22.2, 22.2, 21.1, + 15.8, 17.7, 14.5, 15.4, 19.6, 19.7, 21.8, 14.8, 20.3, 18.9, 10.8, 18.6, + 12.0, 16.1, 17.3, 25.4, 15.3, 16.7, 25.5, 25.5, 20.1, 21.3, 21.7, 23.2, + 16.5, 6.9, 16.1, 16.2, 15.7, 17.6, 22.8, 15.3, 16.4, 17.7, 14.4, 14.9, + ]; + let r = gtheory_pio(&x, 5, 4, 3, &[(6, 2)]).unwrap(); + let want_ss = [ + 208.164, + 128.8205, + 45.6333333333, + 82.792, + 95.475, + 56.448, + 213.077, + ]; + let want_ms = [ + 52.041, + 42.9401666667, + 22.8166666667, + 6.8993333333, + 11.934375, + 9.408, + 8.8782083333, + ]; + let want_var = [ + 3.507125, + 2.3674027778, + 0.517625, + -0.659625, + 0.7640416667, + 0.1059583333, + 8.8782083333, + ]; + for k in 0..7 { + assert_close(r.ss[k], want_ss[k], 1e-8, "ind pio ss"); + assert_close(r.ms[k], want_ms[k], 1e-8, "ind pio ms"); + assert_close(r.var_raw[k], want_var[k], 1e-8, "ind pio var_raw"); + } + // pi component is clamped; the rest pass through. + assert!(r.var_raw[3] < 0.0); + assert_eq!(r.var[3], 0.0); + assert_close(r.var[0], 3.507125, 1e-8, "clamped var_p"); + // D study at (6, 2), computed from CLAMPED components. + let d = &r.d_study[0]; + assert_close(d.rel_error_var, 1.1218715278, 1e-8, "ind delta"); + assert_close(d.abs_error_var, 1.7840810185, 1e-8, "ind Delta"); + assert_close(d.generalizability, 0.7576426076, 1e-8, "ind Erho2"); + assert_close(d.dependability, 0.6628214792, 1e-8, "ind Phi"); +} + +#[test] +fn gt_error_paths() { + // Too few levels. + assert!(gtheory_pi(&[1., 2.], 1, 2, &[1]).is_err()); + assert!(gtheory_pi(&[1., 2.], 2, 1, &[1]).is_err()); + assert!(gtheory_pio(&[0.0; 4], 1, 2, 2, &[(1, 1)]).is_err()); + // Length mismatch. + assert!(gtheory_pi(&[1., 2., 3.], 2, 2, &[1]).is_err()); + assert!(gtheory_pio(&[0.0; 7], 2, 2, 2, &[(1, 1)]).is_err()); + // Non-finite input. + assert!(gtheory_pi(&[1., 2., f64::NAN, 4.], 2, 2, &[1]).is_err()); + // Zero proposed size. + assert!(gtheory_pi(&[1., 2., 3., 5.], 2, 2, &[0]).is_err()); + assert!(gtheory_pio(&[0.0; 8], 2, 2, 2, &[(1, 0)]).is_err()); +} + +/// 500-rep Monte Carlo: simulate the fully random p x i model with known +/// components and check that the mean raw ANOVA estimates are close to the +/// truth (the estimators are unbiased for this model). All quantities under +/// test come from `gtheory_pi` outputs. +#[test] +#[ignore] +fn gt_pi_mc_recovery_500() { + // xorshift-style LCG + Box-Muller, matching the crate's test convention + // of dependency-free generators. + let mut state: u64 = 0x9E3779B97F4A7C15; + let mut next_u = || { + state = state + .wrapping_mul(6364136223846793005) + .wrapping_add(1442695040888963407); + ((state >> 11) as f64) / ((1u64 << 53) as f64) + }; + let mut next_normal = || { + let (u1, u2): (f64, f64) = (next_u().max(1e-12), next_u()); + (-2.0 * u1.ln()).sqrt() * (2.0 * std::f64::consts::PI * u2).cos() + }; + let (n_p, n_i) = (60usize, 12usize); + let (sp, si, spi) = (1.5f64, 0.8f64, 0.6f64); // true variances + let reps = 500; + let mut sums = [0.0f64; 3]; + for _ in 0..reps { + let mut x = vec![0.0; n_p * n_i]; + let pe: Vec = (0..n_p).map(|_| next_normal() * sp.sqrt()).collect(); + let ie: Vec = (0..n_i).map(|_| next_normal() * si.sqrt()).collect(); + for p in 0..n_p { + for i in 0..n_i { + x[p * n_i + i] = 10.0 + pe[p] + ie[i] + next_normal() * spi.sqrt(); + } + } + let r = gtheory_pi(&x, n_p, n_i, &[1]).unwrap(); + for k in 0..3 { + sums[k] += r.var_raw[k]; + } + } + let means = sums.map(|s| s / reps as f64); + // Unbiasedness: generous MC tolerances at 500 reps. + assert!( + (means[0] - sp).abs() < 0.10, + "var_p bias: {}", + means[0] - sp + ); + assert!( + (means[1] - si).abs() < 0.12, + "var_i bias: {}", + means[1] - si + ); + assert!( + (means[2] - spi).abs() < 0.02, + "var_pi bias: {}", + means[2] - spi + ); +} + +/// Regression: impl-review found the old `ss_pio = ss_total - ...` +/// subtraction form returns a NEGATIVE residual SS from catastrophic +/// cancellation on large-offset additive data (`x = 1e12 + p + i + o`, +/// true pio interaction = 0; observed ss[6] = -0.037 before the fix). +/// Asserts read the crate ss/var outputs; kills reintroduction of the +/// subtraction form. +#[test] +fn gt_pio_large_offset_no_cancellation() { + let (n_p, n_i, n_o) = (5usize, 4usize, 3usize); + let mut x = vec![0.0; n_p * n_i * n_o]; + for p in 0..n_p { + for i in 0..n_i { + for o in 0..n_o { + x[p * n_i * n_o + i * n_o + o] = 1e12 + p as f64 + i as f64 + o as f64; + } + } + } + let r = gtheory_pio(&x, n_p, n_i, n_o, &[(4, 3)]).unwrap(); + assert!( + r.ss[6] >= 0.0 && r.ss[6] < 1e-3, + "ss_pio must be ~0 and non-negative, got {}", + r.ss[6] + ); + // Interaction components are exactly-zero in truth; clamped var stays 0. + for k in 3..7 { + assert!( + r.var[k].abs() < 1e-3, + "interaction var[{k}] should be ~0, got {}", + r.var[k] + ); + } +}