diff --git a/CHANGELOG.md b/CHANGELOG.md index bf870b754..27b107a91 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -114,6 +114,26 @@ ### Added +- **Mokken scale analysis** (`fast_mlsirm.mokken_analysis`; new + `mlsirm_core::mokken`; Mokken, 1971, as cited in van der Ark, 2007). + Computes the Loevinger scalability coefficients `Hij`, `Hi`, `H` and their + Mokken Z statistics, and partitions items into Mokken scales with the + automated item selection procedure (AISP, "search normal"), with sample + statistics and selection mechanics verified line-by-line against the mokken + R package source (van der Ark, 2007; Straat et al., 2013): `Hij = + S_ij/Smax_ij` with `Smax` from the comonotone (sorted-column) coupling, + `Hi`/`H` as ratios of pairwise sums, and per-scale Bonferroni-adjusted Z + gates. For LLM-as-a-Judge item-quality management this flags evaluation + items that fail to scale (label 0) and detects multidimensional item pools + before parametric calibration. Complete integer data required (dichotomous + or polytomous). Rust-only numerics; the Python wrapper validates and + marshals. Tests include a brute-force covariance oracle, an exact Guttman + `H = 1` anchor, a hand-computed Z fixture, a Z-gate anchor whose deletion + seeds a spurious scale (this test caught a real sign error in the normal + quantile during development), a hand-constructed Criterion-1 design whose + negative-`Hij` exclusion is the only active gate (mutation-verified), a + two-cluster AISP recovery, and an `#[ignore]` 500-replicate Monte Carlo + (normal + skewed traits). - **Many-Facet Rasch Model (MFRM) rater-severity calibration** (`fast_mlsirm.fit_facets`; new `mlsirm_core::facets`; Linacre, 1989; Eckes, 2015). Fits `ln[P(k)/P(k-1)] = theta_p - d_i - c_j - f_k` — the rating scale model diff --git a/crates/fast-mlsirm-py/src/lib.rs b/crates/fast-mlsirm-py/src/lib.rs index 437e33ea8..5c8d98852 100644 --- a/crates/fast-mlsirm-py/src/lib.rs +++ b/crates/fast-mlsirm-py/src/lib.rs @@ -60,6 +60,7 @@ use mlsirm_core::rasch_cml::{ andersen_lr_test as core_andersen_lr, fit_rasch_cml as core_fit_rasch_cml, }; use mlsirm_core::facets::fit_facets as core_fit_facets; +use mlsirm_core::mokken::{aisp as core_mokken_aisp, coef_h as core_mokken_coef_h}; use mlsirm_core::rsm::fit_rsm as core_fit_rsm; use mlsirm_core::rt::{ fit_rt_lognormal as core_fit_rt, rt_person_fit as core_rt_person_fit, RtConfig, @@ -1470,6 +1471,49 @@ fn fit_facets( Ok(out.into()) } +/// Mokken scalability coefficients (`mlsirm_core::mokken::coef_h`). +/// `x` is a row-major complete `n_persons * n_items` integer score matrix. +/// Returns a dict with `hij`/`zij` (flattened `J*J`, NaN diagonal), `hi`, +/// `zi` (`J`), and scalars `h`, `z`. Sample statistics follow the mokken R +/// package (van der Ark, 2007, https://doi.org/10.18637/jss.v020.i11). +#[pyfunction] +#[pyo3(signature = (x, n_persons, n_items))] +fn mokken_coef_h( + py: Python<'_>, + x: PyReadonlyArray1<'_, i64>, + n_persons: usize, + n_items: usize, +) -> PyResult> { + let res = core_mokken_coef_h(x.as_slice()?, n_persons, n_items) + .map_err(PyValueError::new_err)?; + let out = pyo3::types::PyDict::new(py); + out.set_item("hij", res.hij)?; + out.set_item("hi", res.hi)?; + out.set_item("h", res.h)?; + out.set_item("zij", res.zij)?; + out.set_item("zi", res.zi)?; + out.set_item("z", res.z)?; + Ok(out.into()) +} + +/// Mokken automated item selection procedure (`mlsirm_core::mokken::aisp`, +/// the "search normal" algorithm of the mokken R package). Returns per-item +/// scale labels: 0 = unscalable, 1, 2, ... in formation order. `c` is the +/// scalability lower bound (rule of thumb 0.3), `alpha` the nominal +/// significance level. +#[pyfunction] +#[pyo3(signature = (x, n_persons, n_items, c = 0.3, alpha = 0.05))] +fn mokken_aisp( + x: PyReadonlyArray1<'_, i64>, + n_persons: usize, + n_items: usize, + c: f64, + alpha: f64, +) -> PyResult> { + core_mokken_aisp(x.as_slice()?, n_persons, n_items, c, alpha).map_err(PyValueError::new_err) +} + + /// Marginal-EM fit of a mixed Rasch / mixture-IRT model (`mlsirm_core::mixture`, Rost, /// 1990). `y`/`observed` are row-major `n_persons * n_items`; `model` is "rasch" or /// "2pl". `n_classes` latent classes each get their own item parameters. Returns a dict @@ -5066,6 +5110,8 @@ fn fast_mlsirm_core(m: &Bound<'_, PyModule>) -> PyResult<()> { m.add_function(wrap_pyfunction!(fit_crm, m)?)?; m.add_function(wrap_pyfunction!(fit_rsm, m)?)?; m.add_function(wrap_pyfunction!(fit_facets, m)?)?; + m.add_function(wrap_pyfunction!(mokken_coef_h, m)?)?; + m.add_function(wrap_pyfunction!(mokken_aisp, m)?)?; m.add_function(wrap_pyfunction!(fit_mixture, m)?)?; m.add_function(wrap_pyfunction!(fit_lltm, m)?)?; m.add_function(wrap_pyfunction!(fit_testlet, m)?)?; diff --git a/crates/mlsirm-core/src/lib.rs b/crates/mlsirm-core/src/lib.rs index e55b47ff6..748091f05 100644 --- a/crates/mlsirm-core/src/lib.rs +++ b/crates/mlsirm-core/src/lib.rs @@ -14,6 +14,7 @@ pub mod mhrm; pub mod mixed; pub mod mixture; pub mod mmle; +pub mod mokken; pub mod nodes; pub mod nominal; pub mod oakes; diff --git a/crates/mlsirm-core/src/mokken.rs b/crates/mlsirm-core/src/mokken.rs new file mode 100644 index 000000000..934d96746 --- /dev/null +++ b/crates/mlsirm-core/src/mokken.rs @@ -0,0 +1,377 @@ +//! Mokken scale analysis: Loevinger scalability coefficients and the +//! automated item selection procedure (AISP). +//! +//! Implements, for a complete integer response matrix (dichotomous or +//! polytomous), the sample scalability coefficients +//! +//! ```text +//! Hij = S_ij / Smax_ij +//! Hi = sum_{j != i} S_ij / sum_{j != i} Smax_ij +//! H = sum_{i < j} S_ij / sum_{i < j} Smax_ij +//! ``` +//! +//! where `S` is the sample covariance matrix (denominator N-1) and +//! `Smax_ij = cov(sort(X_i), sort(X_j))` is the maximum covariance +//! attainable given the two items' marginal score distributions — the +//! comonotone (sorted-sorted) coupling maximizes `sum x_p y_p` by the +//! rearrangement inequality, and the means are marginal-fixed, so it +//! maximizes the covariance; the N-1 denominators cancel in every ratio. +//! +//! `Hi` uses the ratio of PAIRWISE sums, exactly as the mokken R package +//! computes it (`coefHTiny`). Verified caveat: this is NOT generally equal +//! to a "max Cov(X_j, R_-j) holding the realized rest-score marginal fixed" +//! reading of van der Ark (2007, Eq. 2) — counterexample: X1=X2=[0,0,1,1], +//! X3=[0,1,0,1] gives fixed-marginal max 1/3 for item 1 but pairwise-sum +//! denominator 2/3. The pairwise-sum form is the de-facto MSA standard and +//! is what this module implements. +//! +//! Mokken's Z statistics (null hypothesis of inter-item independence) follow +//! the mokken package's `coefZ` (`type.z = "Z"`): +//! +//! ```text +//! Zij = S_ij * sqrt(N-1) / sqrt(s_ii * s_jj) +//! Zi = (sum_{j != i} S_ij) * sqrt(N-1) / sqrt(sum_{j != i} s_ii * s_jj) +//! Z = (sum_{i < j} S_ij) * sqrt(N-1) / sqrt(sum_{i < j} s_ii * s_jj) +//! ``` +//! +//! The AISP ("search normal") partitions items into Mokken scales: a start +//! pair maximizing `Hij` among pairs significantly positive (`|Zij| >= Z_c`) +//! with pair `H >= c`, then repeatedly adds the free item that (1) has no +//! negative `Hij` with any selected item (nonnegative allowed), (2) has +//! within-augmented-set `Hi >= c`, (3) has `Zi >= Z_c`, and (4) maximizes the +//! augmented set's total `H`; the scale closes when the best augmented-set +//! `H < c`, and further scales are formed from leftover items. The +//! significance level is Bonferroni-adjusted per scale as +//! `alpha / (K1*(K1-1)/2 + sum of later step candidate counts)`, with the +//! candidate-count vector resetting at each new scale, matching +//! `search.normal.R` (`adjusted.alpha`). +//! +//! Verification status: the coefficient definitions, rules of thumb, and the +//! Mokken-scale definition (all inter-item covariances nonnegative in the +//! selection sense and `Hi >= c > 0`) were read in van der Ark (2007) and +//! Straat et al. (2013); the exact sample statistics, Z forms, tie-breaking, +//! and AISP mechanics were verified line-by-line against the mokken R package +//! source (CRAN, `R/internalFunctions.R::coefHTiny`, `R/coefZ.R`, +//! `R/search.normal.R`). Mokken (1971) and Sijtsma & Molenaar (2002) were NOT +//! read directly; claims from them are relayed via the above sources. No +//! primary-source derivation of the Z normal approximation was verified; +//! it is implementation-verified only. +//! +//! References (APA 7th ed.): +//! - van der Ark, L. A. (2007). Mokken scale analysis in R. *Journal of +//! Statistical Software, 20*(11), 1-19. https://doi.org/10.18637/jss.v020.i11 +//! - Straat, J. H., van der Ark, L. A., & Sijtsma, K. (2013). Comparing +//! optimization algorithms for item selection in Mokken scale analysis. +//! *Journal of Classification, 30*(1), 75-99. +//! https://doi.org/10.1007/s00357-013-9122-y +//! - Mokken, R. J. (1971). *A theory and procedure of scale analysis*. +//! De Gruyter. (as cited in van der Ark, 2007, and Straat et al., 2013) +//! - Sijtsma, K., & Molenaar, I. W. (2002). *Introduction to nonparametric +//! item response theory*. Sage. (as cited in Straat et al., 2013) + +/// Scalability coefficients and Mokken Z statistics for one item set. +#[derive(Debug, Clone)] +pub struct MokkenH { + /// Row-major `n_items x n_items`; `hij[i*J + j] = Hij`, diagonal = NaN. + pub hij: Vec, + /// Per-item scalability `Hi`. + pub hi: Vec, + /// Total scale coefficient `H`. + pub h: f64, + /// Row-major `n_items x n_items` Mokken Z; diagonal = NaN. + pub zij: Vec, + /// Per-item Z. + pub zi: Vec, + /// Total Z. + pub z: f64, +} + +fn validate(x: &[i64], n_persons: usize, n_items: usize) -> Result<(), String> { + if n_persons < 3 { + return Err("mokken requires at least 3 persons".to_string()); + } + if n_items < 2 { + return Err("mokken requires at least 2 items".to_string()); + } + let expected = crate::checked_mul_usize(n_persons, n_items, "n_persons * n_items overflows usize")?; + if x.len() != expected { + return Err(format!( + "responses length {} != n_persons*n_items {}", + x.len(), + expected + )); + } + if x.iter().any(|&v| v < 0) { + return Err("scores must be nonnegative integers".to_string()); + } + Ok(()) +} + +/// Pairwise machinery shared by `coef_h` and `aisp`: covariance matrix `s`, +/// sorted-column max-covariance matrix `smax`, and per-column variances +/// (diagonal of `s`). All with denominator N-1. +fn pairwise(x: &[i64], n_persons: usize, n_items: usize) -> Result<(Vec, Vec), String> { + let n = n_persons as f64; + let j = n_items; + // column means and centered columns; sorted centered columns for smax + let mut cols: Vec> = Vec::with_capacity(j); + let mut sorted: Vec> = Vec::with_capacity(j); + for it in 0..j { + let mut c: Vec = (0..n_persons).map(|p| x[p * j + it] as f64).collect(); + let mean = c.iter().sum::() / n; + for v in c.iter_mut() { + *v -= mean; + } + let mut s = c.clone(); + s.sort_by(|a, b| a.partial_cmp(b).expect("finite")); + cols.push(c); + sorted.push(s); + } + let denom = n - 1.0; + let mut s = vec![0.0; j * j]; + let mut smax = vec![0.0; j * j]; + for a in 0..j { + for b in a..j { + let cov = cols[a] + .iter() + .zip(cols[b].iter()) + .map(|(u, v)| u * v) + .sum::() + / denom; + let cmx = sorted[a] + .iter() + .zip(sorted[b].iter()) + .map(|(u, v)| u * v) + .sum::() + / denom; + s[a * j + b] = cov; + s[b * j + a] = cov; + smax[a * j + b] = cmx; + smax[b * j + a] = cmx; + } + if s[a * j + a] <= 0.0 { + return Err(format!("item {a} has zero variance")); + } + } + Ok((s, smax)) +} + +/// H and Z coefficients for the item subset `idx` (crate-internal; `idx` +/// indexes into the full `s`/`smax` matrices of width `j_full`). +fn h_subset(s: &[f64], smax: &[f64], j_full: usize, idx: &[usize], n_persons: usize) -> (Vec, f64, Vec, f64) { + let k = idx.len(); + let sqn = ((n_persons - 1) as f64).sqrt(); + let mut hi = vec![0.0; k]; + let mut zi = vec![0.0; k]; + let (mut num, mut den, mut vsum) = (0.0, 0.0, 0.0); + for (a, &ia) in idx.iter().enumerate() { + let (mut na, mut da, mut va) = (0.0, 0.0, 0.0); + for &ib in idx.iter() { + if ia == ib { + continue; + } + na += s[ia * j_full + ib]; + da += smax[ia * j_full + ib]; + va += s[ia * j_full + ia] * s[ib * j_full + ib]; + } + hi[a] = na / da; + zi[a] = na * sqn / va.sqrt(); + num += na; + den += da; + vsum += va; + } + // each unordered pair counted twice in the row sums + (hi, num / den, zi, (num / 2.0) * sqn / (vsum / 2.0).sqrt()) +} + +/// Compute `Hij`, `Hi`, `H` and the Mokken Z statistics for a complete +/// `n_persons x n_items` row-major integer score matrix. +pub fn coef_h(x: &[i64], n_persons: usize, n_items: usize) -> Result { + validate(x, n_persons, n_items)?; + let (s, smax) = pairwise(x, n_persons, n_items)?; + let j = n_items; + let sqn = ((n_persons - 1) as f64).sqrt(); + let mut hij = vec![f64::NAN; j * j]; + let mut zij = vec![f64::NAN; j * j]; + for a in 0..j { + for b in 0..j { + if a != b { + hij[a * j + b] = s[a * j + b] / smax[a * j + b]; + zij[a * j + b] = s[a * j + b] * sqn / (s[a * j + a] * s[b * j + b]).sqrt(); + } + } + } + let all: Vec = (0..j).collect(); + let (hi, h, zi, z) = h_subset(&s, &smax, j, &all, n_persons); + Ok(MokkenH { hij, hi, h, zij, zi, z }) +} + +/// Standard-normal upper quantile via inverse complementary error function +/// (Acklam-style rational approximation; |error| < 1.15e-9, sufficient for +/// an alpha cut-off). Returns z such that P(N(0,1) > z) = p. +fn normal_upper_quantile(p: f64) -> f64 { + // invert the CDF at 1 - p using Peter Acklam's approximation + let q = 1.0 - p; + debug_assert!(q > 0.0 && q < 1.0); + const A: [f64; 6] = [ + -3.969683028665376e+01, + 2.209460984245205e+02, + -2.759285104469687e+02, + 1.383577518672690e+02, + -3.066479806614716e+01, + 2.506628277459239e+00, + ]; + const B: [f64; 5] = [ + -5.447609879822406e+01, + 1.615858368580409e+02, + -1.556989798598866e+02, + 6.680131188771972e+01, + -1.328068155288572e+01, + ]; + const C: [f64; 6] = [ + -7.784894002430293e-03, + -3.223964580411365e-01, + -2.400758277161838e+00, + -2.549732539343734e+00, + 4.374664141464968e+00, + 2.938163982698783e+00, + ]; + const D: [f64; 4] = [ + 7.784695709041462e-03, + 3.224671290700398e-01, + 2.445134137142996e+00, + 3.754408661907416e+00, + ]; + let plow = 0.02425; + // standard Acklam sign convention: lower branch yields negative values, + // central passes through, upper is the negated lower expression. + if q < plow { + let r = (-2.0 * q.ln()).sqrt(); + (((((C[0] * r + C[1]) * r + C[2]) * r + C[3]) * r + C[4]) * r + C[5]) + / ((((D[0] * r + D[1]) * r + D[2]) * r + D[3]) * r + 1.0) + } else if q <= 1.0 - plow { + let r = q - 0.5; + let t = r * r; + (((((A[0] * t + A[1]) * t + A[2]) * t + A[3]) * t + A[4]) * t + A[5]) * r + / (((((B[0] * t + B[1]) * t + B[2]) * t + B[3]) * t + B[4]) * t + 1.0) + } else { + let r = (-2.0 * (1.0 - q).ln()).sqrt(); + -((((((C[0] * r + C[1]) * r + C[2]) * r + C[3]) * r + C[4]) * r + C[5]) + / ((((D[0] * r + D[1]) * r + D[2]) * r + D[3]) * r + 1.0)) + } +} + +/// Automated item selection procedure (Mokken's "search normal" AISP). +/// +/// Returns a per-item scale label: 0 = unscalable, 1, 2, ... in formation +/// order. `c` is the scalability lower bound (rule of thumb 0.3); `alpha` the +/// nominal significance level (default 0.05 in the literature). +pub fn aisp( + x: &[i64], + n_persons: usize, + n_items: usize, + c: f64, + alpha: f64, +) -> Result, String> { + validate(x, n_persons, n_items)?; + if !(0.0..1.0).contains(&c) { + return Err("lower bound c must be in [0, 1)".to_string()); + } + if !(alpha > 0.0 && alpha < 1.0) { + return Err("alpha must be in (0, 1)".to_string()); + } + let (s, smax) = pairwise(x, n_persons, n_items)?; + let j = n_items; + let sqn = ((n_persons - 1) as f64).sqrt(); + let hij = |a: usize, b: usize| s[a * j + b] / smax[a * j + b]; + let zij = |a: usize, b: usize| s[a * j + b] * sqn / (s[a * j + a] * s[b * j + b]).sqrt(); + + let mut in_set = vec![0u32; j]; + let mut scale = 0u32; + loop { + scale += 1; + let free: Vec = (0..j).filter(|&i| in_set[i] == 0).collect(); + if free.len() < 2 { + break; + } + // Bonferroni accumulation: k_counts[0] = K1 = #free at scale start; + // later entries are candidate counts of each add step (resets per scale). + let k1 = free.len() as f64; + let mut k_rest = 0.0f64; + let z_c = |k_rest: f64| { + let adj = alpha / (k1 * (k1 - 1.0) * 0.5 + k_rest); + normal_upper_quantile(adj) + }; + // start pair: max Hij among free pairs with |Zij| >= Z_c. Ties mirror + // mokken's eps rule (search.normal.R subtracts row*1e-10 where row is + // the LARGER member index): smaller larger-member index wins, then + // smaller smaller-member index. + let zc0 = z_c(0.0); + let mut best: Option<(usize, usize, f64)> = None; + for (ai, &a) in free.iter().enumerate() { + for &b in free.iter().skip(ai + 1) { + if zij(a, b).abs() < zc0 { + continue; + } + let h = hij(a, b); + let better = match best { + None => true, + Some((ba, bb, bh)) => h > bh || (h == bh && (b, a) < (bb, ba)), + }; + if better { + best = Some((a, b, h)); + } + } + } + let Some((a0, b0, h0)) = best else { break }; + // pair Hi == Hij for both members; require >= c + if h0 < c { + break; + } + let mut selected = vec![a0, b0]; + in_set[a0] = scale; + in_set[b0] = scale; + // add loop + loop { + let candidates: Vec = (0..j) + .filter(|&i| in_set[i] == 0) + .filter(|&i| selected.iter().all(|&sj| hij(i, sj) >= 0.0)) + .collect(); + if candidates.is_empty() { + break; + } + k_rest += candidates.len() as f64; + let zc = z_c(k_rest); + let mut best_h = f64::NEG_INFINITY; + let mut best_item = None; + for &cand in &candidates { + let mut aug = selected.clone(); + aug.push(cand); + let (hi, h_total, zi, _) = h_subset(&s, &smax, j, &aug, n_persons); + // candidate is last in aug + if hi[aug.len() - 1] < c { + continue; + } + if zi[aug.len() - 1] < zc { + continue; + } + if h_total > best_h { + best_h = h_total; + best_item = Some(cand); + } + } + match best_item { + Some(it) if best_h >= c => { + in_set[it] = scale; + selected.push(it); + } + _ => break, + } + } + } + Ok(in_set) +} + +#[cfg(test)] +#[path = "../../../tests/unit/mokken_tests.rs"] +mod tests; diff --git a/python/fast_mlsirm/__init__.py b/python/fast_mlsirm/__init__.py index d3121b1e9..b41f94d50 100644 --- a/python/fast_mlsirm/__init__.py +++ b/python/fast_mlsirm/__init__.py @@ -33,6 +33,7 @@ from .grm import fit_grm as fit_grm, GrmFit as GrmFit from .gpcm import fit_gpcm as fit_gpcm, GpcmFit as GpcmFit from .facets import fit_facets as fit_facets, FacetsFit as FacetsFit +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 from .lltm import fit_lltm as fit_lltm, LltmFit as LltmFit @@ -146,6 +147,8 @@ "RsmFit", "fit_facets", "FacetsFit", + "mokken_analysis", + "MokkenResult", "fit_mixed_items", "MixedFormatFit", "MixedItemParameters", diff --git a/python/fast_mlsirm/mokken.py b/python/fast_mlsirm/mokken.py new file mode 100644 index 000000000..f7b095aa2 --- /dev/null +++ b/python/fast_mlsirm/mokken.py @@ -0,0 +1,109 @@ +"""Mokken scale analysis: Loevinger scalability coefficients and the automated +item selection procedure (AISP). All numeric work happens in the Rust core; +this module only validates and marshals arrays.""" + +from __future__ import annotations + +from dataclasses import dataclass + +import numpy as np + +from .config import MAX_POLYTOMOUS_CATEGORIES + + +@dataclass +class MokkenResult: + """Mokken scalability coefficients and (optionally) AISP scale labels. + + ``hij`` is the ``items x items`` matrix of pairwise scalability + coefficients (NaN diagonal), ``hi`` the per-item coefficients, ``h`` the + total scale coefficient; ``zij``/``zi``/``z`` are the matching Mokken Z + statistics for the null hypothesis of inter-item independence. + ``scale`` holds per-item AISP labels: 0 = unscalable, 1, 2, ... in + formation order. Sample statistics follow the mokken R package + (van der Ark, 2007).""" + + hij: np.ndarray + hi: np.ndarray + h: float + zij: np.ndarray + zi: np.ndarray + z: float + scale: np.ndarray + + +def mokken_analysis( + responses: np.ndarray, + lower_bound: float = 0.3, + alpha: float = 0.05, +) -> MokkenResult: + """Mokken scale analysis (compute in Rust; Mokken, 1971, as cited in + van der Ark, 2007). + + Computes the Loevinger scalability coefficients ``Hij``, ``Hi``, ``H`` + with their Mokken Z statistics, and partitions the items into Mokken + scales with the automated item selection procedure (AISP), following the + sample statistics and "search normal" algorithm of the mokken R package + (van der Ark, 2007): ``Hij = S_ij / Smax_ij`` where ``S`` is the sample + covariance matrix and ``Smax_ij`` the maximum covariance given the two + items' marginals (sorted-column coupling); ``Hi`` and ``H`` are ratios of + the corresponding pairwise sums. A Mokken scale at lower bound ``c`` + requires nonnegative inter-item covariances and ``Hi >= c`` (rule of + thumb ``c = 0.3``; Straat et al., 2013). + + In LLM-as-a-Judge item-quality management, AISP flags evaluation items + that do not scale with the rest (label 0) and detects multidimensional + item pools before parametric IRT calibration. + + ``responses`` is a complete ``persons x items`` array of integer scores + (dichotomous 0/1 or polytomous); missing values are not supported — + Mokken sample statistics assume complete data (van der Ark, 2007). + + References (APA 7th ed.): + van der Ark, L. A. (2007). Mokken scale analysis in R. *Journal of + Statistical Software, 20*(11), 1-19. + https://doi.org/10.18637/jss.v020.i11 + Straat, J. H., van der Ark, L. A., & Sijtsma, K. (2013). Comparing + optimization algorithms for item selection in Mokken scale + analysis. *Journal of Classification, 30*(1), 75-99. + https://doi.org/10.1007/s00357-013-9122-y + Mokken, R. J. (1971). *A theory and procedure of scale analysis*. + De Gruyter. (as cited in van der Ark, 2007) + """ + from .fitstats import _core_module + + core = _core_module() + if core is None or not hasattr(core, "mokken_coef_h"): + raise RuntimeError("mokken_analysis requires the compiled Rust core") + + if not np.isfinite(lower_bound) or not (0.0 <= lower_bound < 1.0): + raise ValueError("lower_bound must be in [0, 1)") + if not np.isfinite(alpha) or not (0.0 < alpha < 1.0): + raise ValueError("alpha must be in (0, 1)") + + y = np.asarray(responses, dtype=np.float64) + if y.ndim != 2: + raise ValueError("responses must be a 2-D persons x items array") + n_persons, n_items = y.shape + if not np.all(np.isfinite(y)): + raise ValueError("responses must be complete (no missing values)") + if np.any(y != np.floor(y)) or np.any(y < 0): + raise ValueError("responses must be non-negative integer scores") + if y.size and int(y.max()) + 1 > MAX_POLYTOMOUS_CATEGORIES: + raise ValueError( + f"responses imply more than {MAX_POLYTOMOUS_CATEGORIES} categories" + ) + x = y.astype(np.int64).reshape(-1) + res = core.mokken_coef_h(x, int(n_persons), int(n_items)) + scale = core.mokken_aisp( + x, int(n_persons), int(n_items), float(lower_bound), float(alpha) + ) + return MokkenResult( + hij=np.asarray(res["hij"], dtype=np.float64).reshape(n_items, n_items), + hi=np.asarray(res["hi"], dtype=np.float64), + h=float(res["h"]), + zij=np.asarray(res["zij"], dtype=np.float64).reshape(n_items, n_items), + zi=np.asarray(res["zi"], dtype=np.float64), + z=float(res["z"]), + scale=np.asarray(scale, dtype=np.int64), + ) diff --git a/tests/test_paper_features.py b/tests/test_paper_features.py index c0a6b2964..9835d9c0f 100644 --- a/tests/test_paper_features.py +++ b/tests/test_paper_features.py @@ -5008,3 +5008,79 @@ def test_fit_facets_rejects_malformed_and_flags_disconnected(): yb[6, 1, 1] = 0.0 resb = fit_facets(yb, n_cat=2, max_iter=50) assert resb.connected is True +def test_mokken_analysis_coefficients_and_cluster_recovery(): + """Mokken scale analysis (van der Ark, 2007): every assert reads crate + outputs (hij/hi/h/zij/scale from mokken_coef_h / mokken_aisp via the + wrapper). A perfect Guttman scalogram must give H = 1 exactly (kills + covmax mutants), and a two-cluster simulation must be partitioned into + exactly two AISP scales (kills selection-logic mutants).""" + import numpy as np + import pytest + from fast_mlsirm import MokkenResult, mokken_analysis + from fast_mlsirm.fitstats import _core_module + + core = _core_module() + if core is None or not hasattr(core, "mokken_coef_h"): + pytest.skip("compiled core built without mokken_coef_h") + + # perfect Guttman scalogram -> H exactly 1 + guttman = np.array( + [[0, 0, 0], [1, 0, 0], [1, 1, 0], [1, 1, 1], [1, 1, 1]], dtype=float + ) + g = mokken_analysis(guttman) + assert isinstance(g, MokkenResult) + assert abs(g.h - 1.0) < 1e-12 + off = ~np.eye(3, dtype=bool) + assert np.allclose(g.hij[off], 1.0) + assert np.all(np.isnan(np.diag(g.hij))) + + # two independent Rasch clusters -> two scales + rng = np.random.default_rng(2013) + n, per = 1500, 4 + bs = np.array([-0.8, -0.3, 0.3, 0.8]) + t1 = rng.normal(size=(n, 1)) * 1.6 + t2 = rng.normal(size=(n, 1)) * 1.6 + xa = (rng.random((n, per)) < 1.0 / (1.0 + np.exp(-(t1 - bs)))).astype(int) + xb = (rng.random((n, per)) < 1.0 / (1.0 + np.exp(-(t2 - bs)))).astype(int) + res = mokken_analysis(np.hstack([xa, xb]), lower_bound=0.3) + a, b = res.scale[0], res.scale[4] + assert a > 0 and b > 0 and a != b, res.scale + assert np.all(res.scale[:4] == a) and np.all(res.scale[4:] == b), res.scale + # crate zij symmetry read-back on real data + assert np.allclose(res.zij[off_idx := ~np.eye(8, dtype=bool)], + res.zij.T[off_idx]) + + +def test_mokken_analysis_rejects_malformed_input(): + """Wrapper validation: incomplete, non-integer, negative, non-2D data and + out-of-range c/alpha raise ValueError (each assert exercises the wrapper + guard in front of the crate; kills guard-deletion mutants).""" + import numpy as np + import pytest + from fast_mlsirm import mokken_analysis + from fast_mlsirm.fitstats import _core_module + + core = _core_module() + if core is None or not hasattr(core, "mokken_coef_h"): + pytest.skip("compiled core built without mokken_coef_h") + + ok = np.array([[0, 1], [1, 0], [1, 1], [0, 0], [1, 0], [0, 1]], dtype=float) + with pytest.raises(ValueError, match="complete"): + bad = ok.copy() + bad[0, 0] = np.nan + mokken_analysis(bad) + with pytest.raises(ValueError, match="integer"): + mokken_analysis(ok + 0.5) + with pytest.raises(ValueError, match="non-negative"): + mokken_analysis(ok - 1.0) + with pytest.raises(ValueError, match="2-D"): + mokken_analysis(ok.reshape(-1)) + with pytest.raises(ValueError, match="lower_bound"): + mokken_analysis(ok, lower_bound=1.0) + with pytest.raises(ValueError, match="alpha"): + mokken_analysis(ok, alpha=0.0) + # crate-side guard surfaces as ValueError too (zero-variance item) + cst = ok.copy() + cst[:, 0] = 1.0 + with pytest.raises(ValueError, match="zero variance"): + mokken_analysis(cst) \ No newline at end of file diff --git a/tests/unit/mokken_tests.rs b/tests/unit/mokken_tests.rs new file mode 100644 index 000000000..f93aeef35 --- /dev/null +++ b/tests/unit/mokken_tests.rs @@ -0,0 +1,391 @@ +//! Tests for Mokken scale analysis (`mlsirm_core::mokken`). +//! +//! Every assert reads values returned by the crate (`MokkenH` fields or the +//! `aisp` label vector). Each test names the crate value it reads and a +//! mutant it kills. + +use super::{aisp, coef_h, normal_upper_quantile}; + +/// Reads: `normal_upper_quantile` directly against published anchors +/// Phi^-1(0.95) = 1.6448536..., Phi^-1(0.999) = 3.0902323... . +/// Kills: sign/branch flips in the Acklam approximation (one such flip was +/// caught by `aisp_z_gate_blocks_insignificant_pair` during development). +#[test] +fn normal_quantile_matches_published_anchors() { + assert!((normal_upper_quantile(0.05) - 1.6448536269514722).abs() < 1e-8); + assert!((normal_upper_quantile(0.001) - 3.090232306167813).abs() < 1e-8); + assert!((normal_upper_quantile(0.5)).abs() < 1e-8); +} + +/// Deterministic xorshift for simulation without external deps. +struct Rng(u64); +impl Rng { + fn new(seed: u64) -> Self { + Rng(seed.max(1)) + } + fn next_f64(&mut self) -> f64 { + let mut x = self.0; + x ^= x << 13; + x ^= x >> 7; + x ^= x << 17; + self.0 = x; + (x >> 11) as f64 / (1u64 << 53) as f64 + } + /// Standard normal via Box-Muller. + fn next_normal(&mut self) -> f64 { + let u1 = self.next_f64().max(1e-12); + let u2 = self.next_f64(); + (-2.0 * u1.ln()).sqrt() * (std::f64::consts::TAU * u2).cos() + } +} + +/// Simulate Rasch data: P(X=1) = logistic(theta - b). +fn simulate_rasch(rng: &mut Rng, n: usize, bs: &[f64], theta_scale: f64) -> Vec { + let j = bs.len(); + let mut x = vec![0i64; n * j]; + for p in 0..n { + let th = rng.next_normal() * theta_scale; + for (i, &b) in bs.iter().enumerate() { + let pr = 1.0 / (1.0 + (-(th - b)).exp()); + x[p * j + i] = if rng.next_f64() < pr { 1 } else { 0 }; + } + } + x +} + +/// Brute-force oracle: two-pass covariance and sorted-column covariance, +/// computed with a DIFFERENT code path (per-pair, f64 accumulation in a +/// different order) than the crate's matrix construction. +fn oracle_pair(x: &[i64], n: usize, j: usize, a: usize, b: usize) -> (f64, f64) { + let col = |it: usize| -> Vec { (0..n).map(|p| x[p * j + it] as f64).collect() }; + let (ca, cb) = (col(a), col(b)); + let (ma, mb) = ( + ca.iter().sum::() / n as f64, + cb.iter().sum::() / n as f64, + ); + let cov = (0..n).map(|p| (ca[p] - ma) * (cb[p] - mb)).sum::() / (n as f64 - 1.0); + let mut sa = ca.clone(); + let mut sb = cb.clone(); + sa.sort_by(|u, v| u.partial_cmp(v).unwrap()); + sb.sort_by(|u, v| u.partial_cmp(v).unwrap()); + let cmx = (0..n).map(|p| (sa[p] - ma) * (sb[p] - mb)).sum::() / (n as f64 - 1.0); + (cov, cmx) +} + +/// Reads: `MokkenH::hij`, `hi`, `h` on random polytomous data, compared to a +/// brute-force oracle computed by an independent path. +/// Kills: any algebra mutant in `pairwise`/`h_subset` (wrong denominator, +/// mean subtraction, sorted-column pairing, aggregation order). +#[test] +fn coefficients_match_brute_force_oracle() { + let mut rng = Rng::new(42); + let (n, j) = (60, 4); + // polytomous 0..=3 with item-varying marginals + let mut x = vec![0i64; n * j]; + for p in 0..n { + let th = rng.next_normal(); + for i in 0..j { + let mut score = 0i64; + for k in 0..3 { + let cut = -1.0 + i as f64 * 0.4 + k as f64 * 0.8; + if th + rng.next_normal() * 0.7 > cut { + score += 1; + } + } + x[p * j + i] = score; + } + } + let res = coef_h(&x, n, j).expect("fit"); + let mut num_tot = 0.0; + let mut den_tot = 0.0; + for a in 0..j { + let mut num_i = 0.0; + let mut den_i = 0.0; + for b in 0..j { + if a == b { + assert!(res.hij[a * j + b].is_nan()); + continue; + } + let (cov, cmx) = oracle_pair(&x, n, j, a, b); + assert!( + (res.hij[a * j + b] - cov / cmx).abs() < 1e-12, + "Hij[{a},{b}] crate {} oracle {}", + res.hij[a * j + b], + cov / cmx + ); + num_i += cov; + den_i += cmx; + if b > a { + num_tot += cov; + den_tot += cmx; + } + } + assert!((res.hi[a] - num_i / den_i).abs() < 1e-12, "Hi[{a}]"); + } + assert!((res.h - num_tot / den_tot).abs() < 1e-12, "H"); +} + +/// Reads: `MokkenH::h` and `hij` on a perfect Guttman scalogram. +/// Kills: covmax mutants — any error in the sorted-column max covariance +/// breaks the exact H = 1 identity (for a nested dichotomous scalogram every +/// observed pair is already comonotone, so S_ij = Smax_ij). +#[test] +fn perfect_guttman_scalogram_has_h_one() { + // 5 persons x 3 items, nested pattern + let x = vec![ + 0, 0, 0, // + 1, 0, 0, // + 1, 1, 0, // + 1, 1, 1, // + 1, 1, 1, // + ]; + let res = coef_h(&x, 5, 3).expect("fit"); + assert!((res.h - 1.0).abs() < 1e-12, "H = {}", res.h); + for a in 0..3 { + for b in 0..3 { + if a != b { + assert!((res.hij[a * 3 + b] - 1.0).abs() < 1e-12); + } + } + } +} + +/// Reads: `MokkenH::zij`, `z` on a hand-computed 2-item fixture +/// (X = [0,0,0,1,1,1], Y = [0,0,1,0,1,1], N = 6); the exact hand derivation +/// is in the test body. +/// Kills: sqrt(N-1) and variance-product mutants in the Z formula. +#[test] +fn z_statistic_matches_hand_computation() { + let x = vec![ + 0, 0, // + 0, 0, // + 0, 1, // + 1, 0, // + 1, 1, // + 1, 1, // + ]; + let res = coef_h(&x, 6, 2).expect("fit"); + // hand: means .5/.5; centered cross products: + // (-.5)(-.5)*2 + (-.5)(.5) + (.5)(-.5) + (.5)(.5)*2 = .5 - .5 + .5 = 0.5 + // S_xy = 0.5/5 = 0.1 ; s_xx = s_yy = (6*.25)/5 = 0.3 + // Smax: sorted-sorted = comonotone = 1.5/5 = 0.3 -> Hij = 1/3 + // Zij = 0.1*sqrt(5)/sqrt(0.09) = 0.1*2.23606.../0.3 = 0.745355... + let expect_z = 0.1 * 5f64.sqrt() / 0.3; + assert!((res.hij[1] - 1.0 / 3.0).abs() < 1e-12, "Hij = {}", res.hij[1]); + assert!((res.zij[1] - expect_z).abs() < 1e-12, "Zij = {}", res.zij[1]); + // total Z for 2 items equals Zij + assert!((res.z - expect_z).abs() < 1e-12); +} + +/// Reads: `aisp` labels on the same fixture: Hij = 1/3 > c = 0.3 but +/// Zij ~ 0.745 < z_crit(0.05) = 1.645, so NO scale may form. +/// Kills: deleting the start-pair Z significance gate (mutant seeds a scale +/// because Hij exceeds c). +#[test] +fn aisp_z_gate_blocks_insignificant_pair() { + let x = vec![ + 0, 0, // + 0, 0, // + 0, 1, // + 1, 0, // + 1, 1, // + 1, 1, // + ]; + let labels = aisp(&x, 6, 2, 0.3, 0.05).expect("aisp"); + assert_eq!(labels, vec![0, 0], "Z-gate must block the scale"); +} + +/// Reads: `aisp` labels plus `MokkenH::hij`/`hi` on a hand-constructed 80x3 +/// contingency design (profile counts: 19x(1,1,1), 11x(1,1,0), 10x(1,0,0), +/// 10x(0,1,1), 11x(0,0,1), 19x(0,0,0); all marginals 0.5). By construction +/// H01 = 0.5 (start pair), H12 = 0.45, H02 = -0.05, and candidate item 2 +/// passes every other gate at c = 0.15: Hi(2) = 0.2 >= c, Zi(2) ~ 2.51 > +/// z_crit, augmented H = 0.3 >= c. Only the negative-Hij (Criterion 1) +/// exclusion keeps it out. +/// Kills: removing the `hij >= 0` candidate filter in the add loop (the +/// mutant then admits item 2, flipping labels to [1,1,1]). +#[test] +fn aisp_excludes_candidate_with_negative_hij() { + let profiles: [([i64; 3], usize); 6] = [ + ([1, 1, 1], 19), + ([1, 1, 0], 11), + ([1, 0, 0], 10), + ([0, 1, 1], 10), + ([0, 0, 1], 11), + ([0, 0, 0], 19), + ]; + let mut x = Vec::with_capacity(80 * 3); + for (row, count) in profiles { + for _ in 0..count { + x.extend_from_slice(&row); + } + } + let res = coef_h(&x, 80, 3).expect("fit"); + // verify the construction via crate values: item 2 negative with item 0 + // yet passes the Hi gate at c = 0.15 + assert!(res.hij[2 * 3] < 0.0, "Hij(2,0) = {}", res.hij[2 * 3]); + assert!((res.hij[2 * 3] - (-0.05)).abs() < 1e-12); + assert!((res.hi[2] - 0.2).abs() < 1e-12, "Hi(2) = {}", res.hi[2]); + assert!((res.hij[1] - 0.5).abs() < 1e-12, "Hij(0,1) = {}", res.hij[1]); + let labels = aisp(&x, 80, 3, 0.15, 0.05).expect("aisp"); + assert_eq!(labels, vec![1, 1, 0], "Criterion 1 must exclude item 2"); +} + +/// Reads: `aisp` labels on a two-cluster simulation (two independent Rasch +/// dimensions). AISP at c = 0.3 must recover the two clusters exactly. +/// Kills: selection-logic mutants (wrong argmax, wrong exclusion of previous +/// scales, missing multi-scale restart). +#[test] +fn aisp_recovers_two_clusters() { + let mut rng = Rng::new(2013); + let n = 1500; + let bs = [-0.8, -0.3, 0.3, 0.8]; + // cluster A: items 0..=3 driven by theta1; cluster B: items 4..=7 by theta2 + let j = 8; + let mut x = vec![0i64; n * j]; + for p in 0..n { + let t1 = rng.next_normal() * 1.6; + let t2 = rng.next_normal() * 1.6; + for (i, &b) in bs.iter().enumerate() { + let pr1 = 1.0 / (1.0 + (-(t1 - b)).exp()); + let pr2 = 1.0 / (1.0 + (-(t2 - b)).exp()); + x[p * j + i] = if rng.next_f64() < pr1 { 1 } else { 0 }; + x[p * j + 4 + i] = if rng.next_f64() < pr2 { 1 } else { 0 }; + } + } + let labels = aisp(&x, n, j, 0.3, 0.05).expect("aisp"); + let first = labels[0]; + let second = labels[4]; + assert!(first > 0 && second > 0 && first != second, "labels = {labels:?}"); + assert!(labels[..4].iter().all(|&l| l == first), "{labels:?}"); + assert!(labels[4..].iter().all(|&l| l == second), "{labels:?}"); +} + +/// Reads: `MokkenH` fields for score-translation invariance: adding a +/// constant to every score of an item must not change any coefficient +/// (covariances are translation-invariant). +/// Kills: accidental use of raw (uncentered) moments. +#[test] +fn coefficients_invariant_to_score_translation() { + let mut rng = Rng::new(99); + let n = 200; + let x = simulate_rasch(&mut rng, n, &[-0.5, 0.0, 0.5], 1.3); + let mut shifted = x.clone(); + for p in 0..n { + shifted[p * 3 + 1] += 3; // item 1 scored 3..4 instead of 0..1 + } + let a = coef_h(&x, n, 3).expect("fit"); + let b = coef_h(&shifted, n, 3).expect("fit"); + assert!((a.h - b.h).abs() < 1e-12); + for i in 0..3 { + assert!((a.hi[i] - b.hi[i]).abs() < 1e-12); + assert!((a.zi[i] - b.zi[i]).abs() < 1e-12); + } +} + +/// Reads: error `Result`s from both entry points. +/// Kills: deletion of the validation guards. +#[test] +fn rejects_bad_inputs() { + let ok = vec![0, 1, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1]; + assert!(coef_h(&ok, 2, 2).is_err(), "n_persons < 3"); + assert!(coef_h(&ok[..4], 4, 2).is_err(), "length mismatch"); + assert!(coef_h(&[0, -1, 1, 0, 1, 1], 3, 2).is_err(), "negative score"); + assert!(coef_h(&[1, 0, 1, 1, 1, 0], 3, 2).is_err(), "zero variance item 0"); + assert!(coef_h(&ok, 6, 1).is_err(), "single item"); + assert!(aisp(&ok, 6, 2, 1.2, 0.05).is_err(), "c out of range"); + assert!(aisp(&ok, 6, 2, 0.3, 0.0).is_err(), "alpha out of range"); +} + +/// Reads: `aisp` labels on an exact-tie design: X0 == X3 and X1 == X2 +/// (identical columns, Hij = 1 for both pairs) with the two blocks exactly +/// uncorrelated (balanced half-split vs alternating pattern gives sample +/// cov = 0). mokken's eps tie-break (search.normal.R: penalty row*1e-10 on +/// the LARGER member index) must pick pair {1,2} first, so labels are +/// [2, 1, 1, 2]. +/// Kills: reverting to first-encountered lexicographic tie-breaking, which +/// would start with pair {0,3} and yield [1, 2, 2, 1]. +#[test] +fn aisp_tie_break_matches_mokken_eps_rule() { + let n = 40; + let j = 4; + let mut x = vec![0i64; n * j]; + for p in 0..n { + let a = if p < 20 { 1 } else { 0 }; // half-split + let b = (p % 2) as i64; // alternating; sample cov(a, b) = 0 exactly + x[p * j] = a; + x[p * j + 1] = b; + x[p * j + 2] = b; + x[p * j + 3] = a; + } + let labels = aisp(&x, n, j, 0.3, 0.05).expect("aisp"); + assert_eq!(labels, vec![2, 1, 1, 2], "eps tie-break must favor pair {{1,2}}"); +} + +/// Reads: `aisp` labels; independent items (no common trait) must all remain +/// unscalable at c = 0.3. Smoke check of overall gating (not attributed to a +/// single mutant; the Z-gate kill lives in `aisp_z_gate_blocks_insignificant_pair`). +#[test] +fn aisp_leaves_independent_items_unscaled() { + let mut rng = Rng::new(5); + let n = 500; + let j = 5; + let mut x = vec![0i64; n * j]; + for v in x.iter_mut() { + *v = if rng.next_f64() < 0.5 { 1 } else { 0 }; + } + let labels = aisp(&x, n, j, 0.3, 0.05).expect("aisp"); + assert_eq!(labels, vec![0; j], "labels = {labels:?}"); +} + +/// Monte Carlo: >= 500 replications of a unidimensional Rasch scale +/// (normal and skew-positive traits). Reads crate `h` and `aisp` labels. +/// Asserts distributional behavior: mean H within a plausible band and +/// one-scale full recovery in >= 95% of replications. +/// Limitations stated: this cannot pin exact constants; the algebra anchors +/// live in `coefficients_match_brute_force_oracle` and +/// `z_statistic_matches_hand_computation`. +#[test] +#[ignore] +fn monte_carlo_unidimensional_recovery() { + let bs = [-1.0, -0.5, 0.0, 0.5, 1.0]; + let n = 500; + for (label, skew) in [("normal", false), ("skew", true)] { + let mut full = 0usize; + let mut h_sum = 0.0; + let reps = 500; + for rep in 0..reps { + let mut rng = Rng::new(1000 + rep as u64); + let j = bs.len(); + let mut x = vec![0i64; n * j]; + for p in 0..n { + let mut th = rng.next_normal(); + if skew { + // half-normal shifted: skewed positive trait + th = th.abs() * 1.2 - 0.9; + } + th *= 1.5; + for (i, &b) in bs.iter().enumerate() { + let pr = 1.0 / (1.0 + (-(th - b)).exp()); + x[p * j + i] = if rng.next_f64() < pr { 1 } else { 0 }; + } + } + let res = coef_h(&x, n, bs.len()).expect("fit"); + h_sum += res.h; + let labels = aisp(&x, n, bs.len(), 0.3, 0.05).expect("aisp"); + if labels.iter().all(|&l| l == 1) { + full += 1; + } + } + let mean_h = h_sum / reps as f64; + assert!( + mean_h > 0.35 && mean_h < 0.75, + "{label}: mean H = {mean_h}" + ); + assert!( + full as f64 / reps as f64 >= 0.95, + "{label}: full-recovery rate = {}", + full as f64 / reps as f64 + ); + } +}