Add generalizability theory G/D studies (crossed p x i, p x i x o)

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

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<usize>,
+) -> PyResult<Py<pyo3::types::PyDict>> {
+    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<Py<pyo3::types::PyDict>> {
+    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<Vec<Py<pyo3::types::PyDict>>> {
+    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)?)?;

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<f64>, Vec<f64>) {
     (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<f64>, Vec<f64>)> = std::sync::OnceLock::new();
-
-#[inline]
-fn gl64() -> &'static (Vec<f64>, Vec<f64>) {
-    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<f64>, Vec<f64>) {
 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::<f64>()
             * 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. <https://CRAN.R-project.org/package=betafunctions>
 ///
 /// 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,