diff --git a/CHANGELOG.md b/CHANGELOG.md index 27b107a91..478dde4e2 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -114,6 +114,52 @@ ### Added +- **Haberman subscore added-value analysis** (`fast_mlsirm.subscore_analysis`; + new `mlsirm_core::subscores`; Haberman, 2008, as cited in Sinharay, 2010). + For each subscale of a disjoint, exhaustive item partition computes the + PRMSEs of the three classical-test-theory true-subscore estimators — from + the observed subscore (`= Cronbach alpha`), from the observed total + (`rho^2(s_t, x_t) * alpha_x` with the true-score covariance row sum over + subscore columns only), and from both jointly (Wainer-style augmentation via + `tau`/`beta`/`gamma`) — plus per-person estimator matrices, the + `(K+1)^2` score correlation matrix, disattenuated subscore correlations, and + added-value decisions (Haberman's `PRMSE_s > PRMSE_x`; Sinharay's 2010 + `+ 0.01` margin for augmentation, labeled — CRAN `CTTsub`'s relative rule is + documented but not implemented). Formulas verified against the Appendix of + Sinharay (2010, ETS RR-10-16) and the CRAN `subscore` R source read + line-by-line; degenerate samples (alpha outside `(0, 1]`, zero variance, + subscore collinear with the total) are rejected instead of propagating NaN. + For LLM-as-a-Judge item-quality management this decides whether per-domain + judge subscores add diagnostic value over the overall score. Rust-only + numerics; the Python wrapper validates and marshals. Tests pin every + reported statistic against literals from an independent NumPy transcription + of the R semantics on an asymmetric fixture with mixed added-value + outcomes, include rejection tests for the structural and degeneracy guards + (the defensive computed-PRMSE-range guard is not separately exercised), a + conditional dominance + sweep on guard-passing random data, and a 500-rep `#[ignore]` Monte Carlo + MSE comparison; three mutation spot-checks (dropped `m/(m-1)`, rowsum + including the total column, `tau` numerator sign flip) were run and killed. +- **Kernel-smoothing nonparametric IRT** (`fast_mlsirm.ksirt_analysis`; new + `mlsirm_core::ksirt`; Ramsay, 1991, as cited in Mazza et al., 2014). + Estimates option characteristic curves by Nadaraya-Watson kernel regression + (gaussian/quadratic/uniform kernels) of option indicators on rank-based + ordinal ability estimates `qnorm(rank/(n+1))`, on an equally spaced + evaluation grid, with Silverman-rule default bandwidths, plus expected item + score and expected total score curves. Formulas verified against the + KernSmoothIRT JSS paper (Mazza et al., 2014, Sections 2-2.3) and the + KernSmoothIRT R/C++ package source read line-by-line; standard errors and + cross-validation bandwidth selection are deliberately out of scope (the R + implementation's SE accumulator is order-dependent and unverifiable from + read sources). For LLM-as-a-Judge item-quality management this reveals + non-monotone or poorly discriminating evaluation items without a parametric + model. Rust-only numerics; the Python wrapper validates and marshals. Tests + pin a hand-computed 4-person fixture (rank->theta qnorm literals, grid + endpoints, Silverman constant), enforce structural invariants + (row-sums-to-one with positive denominators, compact-support zeros, + zero-denominator fallback), and include a 500-replication Monte Carlo + recovery study (`#[ignore]`) under normal and skewed ability generation + using the rank-invariance composition oracle. - **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 diff --git a/crates/fast-mlsirm-py/src/lib.rs b/crates/fast-mlsirm-py/src/lib.rs index 5c8d98852..966b49bfe 100644 --- a/crates/fast-mlsirm-py/src/lib.rs +++ b/crates/fast-mlsirm-py/src/lib.rs @@ -60,6 +60,8 @@ 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::ksirt::{ksirt as core_ksirt, KsirtKernel}; +use mlsirm_core::subscores::subscores as core_subscores; 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::{ @@ -1513,6 +1515,124 @@ fn mokken_aisp( core_mokken_aisp(x.as_slice()?, n_persons, n_items, c, alpha).map_err(PyValueError::new_err) } +/// Kernel-smoothing nonparametric option characteristic curves +/// (`mlsirm_core::ksirt`; Ramsay, 1991, as cited in Mazza, Punzo, & +/// McGuire, 2014, https://doi.org/10.18637/jss.v058.i06). `x` is a +/// row-major complete `n_persons * n_items` pre-scored response matrix. +/// Returns a dict with `theta` (`N`), `grid` (`Q`), `bandwidth` (`J`), and +/// per-item lists `options`, `occ` (flattened `m_j * Q`, row-major by +/// option), `expected` (`Q`), plus `expected_total` (`Q`). +#[pyfunction] +#[pyo3(signature = (x, n_persons, n_items, kernel = "gaussian", nevalpoints = 51, bandwidth = None))] +fn ksirt_occ( + py: Python<'_>, + x: PyReadonlyArray1<'_, f64>, + n_persons: usize, + n_items: usize, + kernel: &str, + nevalpoints: usize, + bandwidth: Option>, +) -> PyResult> { + let flat = x.as_slice()?; + if flat.len() != n_persons * n_items { + return Err(PyValueError::new_err(format!( + "x has {} entries, expected n_persons * n_items = {}", + flat.len(), + n_persons * n_items + ))); + } + let kern = match kernel { + "gaussian" => KsirtKernel::Gaussian, + "quadratic" => KsirtKernel::Quadratic, + "uniform" => KsirtKernel::Uniform, + other => { + return Err(PyValueError::new_err(format!( + "unknown kernel '{other}' (expected gaussian, quadratic, or uniform)" + ))) + } + }; + let rows: Vec<&[f64]> = flat.chunks_exact(n_items).collect(); + let res = core_ksirt(&rows, kern, nevalpoints, bandwidth.as_deref()) + .map_err(PyValueError::new_err)?; + let out = pyo3::types::PyDict::new(py); + out.set_item("theta", res.theta)?; + out.set_item("grid", res.grid)?; + out.set_item("bandwidth", res.bandwidth)?; + out.set_item("expected_total", res.expected_total)?; + let options: Vec> = res.items.iter().map(|it| it.options.clone()).collect(); + let occ: Vec> = res + .items + .iter() + .map(|it| it.occ.iter().flatten().copied().collect()) + .collect(); + let expected: Vec> = res.items.iter().map(|it| it.expected.clone()).collect(); + out.set_item("options", options)?; + out.set_item("occ", occ)?; + out.set_item("expected", expected)?; + Ok(out.into()) +} + +/// Haberman subscore added-value analysis (`mlsirm_core::subscores`; +/// Haberman, 2008, as cited in Sinharay, 2010, +/// ETS RR-10-16). `x` is a row-major complete `n_persons * n_items` +/// scored response matrix; `groups[j]` in `0..K` assigns item `j` to a +/// subscale. Returns a dict with per-subscale `alpha`, `prmse_s`, +/// `prmse_x`, `prmse_sx`, `tau`, `beta`, `gamma`, `added_value_s`, +/// `added_value_sx`, `alpha_total`, the `(K+1)^2` flattened `corr`, the +/// `K*K` flattened `disattenuated_corr` (NaN diagonal), and the `n*K` +/// flattened estimator matrices `observed`, `subscore_s`, `subscore_x`, +/// `subscore_sx` plus `total` (`n`). +#[pyfunction] +fn subscore_analysis( + py: Python<'_>, + x: PyReadonlyArray1<'_, f64>, + n_persons: usize, + n_items: usize, + groups: Vec, +) -> PyResult> { + let flat = x.as_slice()?; + // Validate BEFORE allocating rows: unchecked n_persons * n_items can + // wrap on 64-bit (e.g. 2^63 * 2 == 0, matching an empty array) and then + // panic with capacity overflow inside the row allocation. + if n_persons < 3 || n_items < 4 || groups.len() != n_items { + return Err(PyValueError::new_err( + "need n_persons >= 3, n_items >= 4, and one group index per item", + )); + } + let expected = n_persons + .checked_mul(n_items) + .ok_or_else(|| PyValueError::new_err("n_persons * n_items overflows"))?; + if flat.len() != expected { + return Err(PyValueError::new_err(format!( + "x has {} entries, expected n_persons * n_items = {expected}", + flat.len(), + ))); + } + let rows: Vec> = (0..n_persons) + .map(|i| flat[i * n_items..(i + 1) * n_items].to_vec()) + .collect(); + let res = core_subscores(&rows, &groups).map_err(PyValueError::new_err)?; + let out = pyo3::types::PyDict::new(py); + out.set_item("alpha", res.alpha)?; + out.set_item("alpha_total", res.alpha_total)?; + out.set_item("prmse_s", res.prmse_s)?; + out.set_item("prmse_x", res.prmse_x)?; + out.set_item("prmse_sx", res.prmse_sx)?; + out.set_item("tau", res.tau)?; + out.set_item("beta", res.beta)?; + out.set_item("gamma", res.gamma)?; + out.set_item("added_value_s", res.added_value_s)?; + out.set_item("added_value_sx", res.added_value_sx)?; + out.set_item("total", res.total)?; + let flatten = |m: Vec>| -> Vec { m.into_iter().flatten().collect() }; + out.set_item("corr", flatten(res.corr))?; + out.set_item("disattenuated_corr", flatten(res.disattenuated_corr))?; + out.set_item("observed", flatten(res.observed))?; + out.set_item("subscore_s", flatten(res.subscore_s))?; + out.set_item("subscore_x", flatten(res.subscore_x))?; + out.set_item("subscore_sx", flatten(res.subscore_sx))?; + Ok(out.into()) +} /// 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 @@ -5112,6 +5232,8 @@ fn fast_mlsirm_core(m: &Bound<'_, PyModule>) -> PyResult<()> { 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!(ksirt_occ, m)?)?; + m.add_function(wrap_pyfunction!(subscore_analysis, 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/ksirt.rs b/crates/mlsirm-core/src/ksirt.rs new file mode 100644 index 000000000..240436acc --- /dev/null +++ b/crates/mlsirm-core/src/ksirt.rs @@ -0,0 +1,282 @@ +//! Kernel-smoothing nonparametric item response theory (Ramsay-style OCCs). +//! +//! Estimates option characteristic curves (OCCs) by Nadaraya-Watson kernel +//! regression of option indicators on rank-based ordinal ability surrogates, +//! the approach popularized by Ramsay (1991, as cited in Mazza et al., 2014) +//! and implemented in TestGraf and the R package KernSmoothIRT. +//! +//! # Verified sources (citation discipline) +//! +//! Every formula below was verified against sources actually read: +//! +//! - Mazza, Punzo, and McGuire (2014), Sections 2, 2.1, 2.2, 2.3 (full PDF +//! read): rank transform `r_i = rank(t_i)/(n+1)`, ordinal ability +//! `theta_i = F^{-1}(r_i)`, Nadaraya-Watson weights, Gaussian / +//! quadratic / uniform kernels, Silverman bandwidth (their Eq. 5), and +//! the expected item score `e_j(theta) = sum_l x_{jl} p_{jl}(theta)`. +//! - KernSmoothIRT 1.0.3 R/C++ source, read line by line +//! (github.com/cran/KernSmoothIRT): `R/ksIRT.R` (ties.method="first", +//! denominator `n+1`, grid endpoints `F^{-1}(1/(n+1))` to +//! `F^{-1}(n/(n+1))` with 51 default points, `h = 1.06 * sigma * n^{-1/5}` +//! with `sigma = 1` for the normal ability metric) and +//! `src/smoother3.cpp` (per-subject NW smoothing; zero-denominator +//! fallback returns all-zero weights). +//! +//! Ramsay (1991) itself was NOT obtainable and is cited only through Mazza +//! et al. (2014); no formula here is attributed to it directly. +//! +//! # Deliberate scope reductions and divergences +//! +//! - Pointwise standard errors are omitted: the R package's `stderr` +//! accumulates `p(1-p)` from a *partially summed* running estimate +//! (smoother3.cpp lines 148-150, order-dependent), and the JSS paper's +//! Eq. 6 uses a different per-subject form; neither yields a closed form +//! verifiable from the read sources, so v1 ships without SEs. +//! - Cross-validation bandwidth selection is omitted (the R implementation +//! subsamples 10% of subjects at random, making it nondeterministic). +//! - Option lists are reported in ascending score order, whereas R keeps +//! first-seen order; the estimated curves are unaffected. +//! - Responses must be complete and pre-scored (numeric option scores); +//! missing-data handling, answer keys, DIF groups, and non-normal ability +//! metrics are out of scope for v1. +//! +//! # References +//! +//! Mazza, A., Punzo, A., & McGuire, B. (2014). KernSmoothIRT: An R package +//! for kernel smoothing in item response theory. *Journal of Statistical +//! Software, 58*(6), 1-34. https://doi.org/10.18637/jss.v058.i06 +//! +//! Nadaraya, E. A. (1964). On estimating regression. *Theory of Probability +//! & Its Applications, 9*(1), 141-142. (As cited in Mazza et al., 2014.) +//! +//! Ramsay, J. O. (1991). Kernel smoothing approaches to nonparametric item +//! characteristic curve estimation. *Psychometrika, 56*(4), 611-630. +//! https://doi.org/10.1007/BF02294494 (As cited in Mazza et al., 2014.) +//! +//! Silverman, B. W. (1986). *Density estimation for statistics and data +//! analysis*. Chapman & Hall. (As cited in Mazza et al., 2014.) +//! +//! Watson, G. S. (1964). Smooth regression analysis. *Sankhya A, 26*(4), +//! 359-372. (As cited in Mazza et al., 2014.) + +use crate::mokken::normal_upper_quantile; + +/// Kernel function for the Nadaraya-Watson smoother. +/// +/// Formulas verified against Mazza et al. (2014, Section 2) and +/// smoother3.cpp: Gaussian `exp(-u^2/2)`, quadratic `(1-u^2)` on `[-1,1]`, +/// uniform indicator on `[-1,1]`. Multiplicative kernel constants cancel in +/// the NW normalization, so the unnormalized forms match the R package. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum KsirtKernel { + Gaussian, + Quadratic, + Uniform, +} + +impl KsirtKernel { + fn eval(self, u: f64) -> f64 { + match self { + KsirtKernel::Gaussian => (-0.5 * u * u).exp(), + KsirtKernel::Quadratic => { + if u.abs() <= 1.0 { + 1.0 - u * u + } else { + 0.0 + } + } + KsirtKernel::Uniform => { + if u.abs() <= 1.0 { + 1.0 + } else { + 0.0 + } + } + } + } +} + +/// Per-item kernel-smoothing output. +#[derive(Debug, Clone)] +pub struct KsirtItem { + /// Distinct observed option scores, ascending. + pub options: Vec, + /// `options.len() x grid.len()` option characteristic curves; row `l` + /// gives `p_hat_{jl}(theta_s)` over the evaluation grid. Rows sum to 1 + /// at any grid point with a positive weight denominator and to 0 where + /// all kernel weights vanish (compact-support kernels far from data). + pub occ: Vec>, + /// Expected item score curve `sum_l x_{jl} * p_hat_{jl}(theta_s)`. + pub expected: Vec, +} + +/// Result of [`ksirt`]. +#[derive(Debug, Clone)] +pub struct KsirtResult { + /// Ordinal ability surrogates `Phi^{-1}(rank(t_i)/(n+1))`, subject order. + pub theta: Vec, + /// Evaluation grid (equally spaced, `Phi^{-1}(1/(n+1))` to + /// `Phi^{-1}(n/(n+1))`). + pub grid: Vec, + /// Per-item bandwidths actually used. + pub bandwidth: Vec, + /// Per-item OCC estimates. + pub items: Vec, + /// Expected total score curve (sum of per-item expected curves). + pub expected_total: Vec, +} + +/// Kernel smoothing of option characteristic curves. +/// +/// `x[i][j]` is the observed (pre-scored, finite) response of subject `i` +/// to item `j`; the distinct values of column `j` form the option set. +/// `bandwidth` overrides the Silverman default `1.06 * n^{-1/5}` (per-item +/// values, all > 0). See the module docs for the algorithm and sources. +pub fn ksirt( + x: &[R], + kernel: KsirtKernel, + nevalpoints: usize, + bandwidth: Option<&[f64]>, +) -> Result +where + R: AsRef<[f64]>, +{ + let n = x.len(); + if n < 2 { + return Err("ksirt requires at least 2 subjects".to_string()); + } + let k = x[0].as_ref().len(); + if k == 0 { + return Err("ksirt requires at least 1 item".to_string()); + } + for (i, row) in x.iter().enumerate() { + let row = row.as_ref(); + if row.len() != k { + return Err(format!( + "ragged response matrix: row {i} has {} items, expected {k}", + row.len() + )); + } + for (j, &v) in row.iter().enumerate() { + if !v.is_finite() { + return Err(format!("non-finite response at subject {i}, item {j}")); + } + } + } + if nevalpoints < 2 { + return Err("nevalpoints must be at least 2".to_string()); + } + let h: Vec = match bandwidth { + Some(b) => { + if b.len() != k { + return Err(format!( + "bandwidth length {} does not match {} items", + b.len(), + k + )); + } + if b.iter().any(|&v| !(v > 0.0) || !v.is_finite()) { + return Err("bandwidths must be finite and positive".to_string()); + } + b.to_vec() + } + None => { + // Silverman rule, sigma = 1 on the normal ability metric + // (Mazza et al., 2014, Eq. 5; ksIRT.R lines 171-179). + let hs = 1.06 * (n as f64).powf(-0.2); + vec![hs; k] + } + }; + + // Step 1: total scores -> ranks (ties by first occurrence, matching + // R's ties.method="first"; ksIRT.R line 121) -> normal quantiles. + let totals: Vec = x.iter().map(|row| row.as_ref().iter().sum()).collect(); + let mut order: Vec = (0..n).collect(); + // stable sort keeps original subject order within ties => "first" + order.sort_by(|&a, &b| totals[a].partial_cmp(&totals[b]).unwrap()); + let mut rank = vec![0usize; n]; + for (pos, &subj) in order.iter().enumerate() { + rank[subj] = pos + 1; + } + let np1 = (n + 1) as f64; + // Phi^{-1}(r) = normal_upper_quantile(1 - r): the helper returns z with + // P(N(0,1) > z) = p, so upper tail 1-r gives the lower quantile at r. + let theta: Vec = rank + .iter() + .map(|&r| normal_upper_quantile(1.0 - r as f64 / np1)) + .collect(); + + // Step 2: evaluation grid (ksIRT.R lines 134-141). + let lim1 = normal_upper_quantile(1.0 - 1.0 / np1); + let lim2 = normal_upper_quantile(1.0 - n as f64 / np1); + let q = nevalpoints; + let step = (lim2 - lim1) / (q - 1) as f64; + let grid: Vec = (0..q).map(|s| lim1 + step * s as f64).collect(); + + // Steps 3-4: per grid point, NW weights shared across the item's + // options (smoother3.cpp lines 76-154, incl. zero-denominator fallback). + let mut items = Vec::with_capacity(k); + let mut expected_total = vec![0.0; q]; + for j in 0..k { + let mut options: Vec = x.iter().map(|row| row.as_ref()[j]).collect(); + options.sort_by(|a, b| a.partial_cmp(b).unwrap()); + options.dedup(); + let option_index: Vec = x + .iter() + .enumerate() + .map(|(i, row)| { + let value = row.as_ref()[j]; + options + .binary_search_by(|option| { + option.partial_cmp(&value).expect( + "internal invariant violation: ksirt compares only finite \ + responses validated at function entry", + ) + }) + .unwrap_or_else(|_| { + panic!( + "option for subject {i} item {j} value {value} not found \ + in sorted/deduped options (internal invariant violation)" + ) + }) + }) + .collect(); + let m = options.len(); + let mut occ = vec![vec![0.0; q]; m]; + let mut expected = vec![0.0; q]; + for s in 0..q { + let kw: Vec = theta + .iter() + .map(|&t| kernel.eval((grid[s] - t) / h[j])) + .collect(); + let denom: f64 = kw.iter().sum(); + if denom <= 0.0 { + continue; // all weights zero: occ stays 0 (R fallback) + } + for (&l, &w) in option_index.iter().zip(&kw) { + occ[l][s] += w / denom; + } + for l in 0..m { + expected[s] += options[l] * occ[l][s]; + } + expected_total[s] += expected[s]; + } + items.push(KsirtItem { + options, + occ, + expected, + }); + } + + Ok(KsirtResult { + theta, + grid, + bandwidth: h, + items, + expected_total, + }) +} + +#[cfg(test)] +#[path = "../../../tests/unit/ksirt_tests.rs"] +mod tests; diff --git a/crates/mlsirm-core/src/lib.rs b/crates/mlsirm-core/src/lib.rs index 748091f05..213a558aa 100644 --- a/crates/mlsirm-core/src/lib.rs +++ b/crates/mlsirm-core/src/lib.rs @@ -7,6 +7,7 @@ pub mod facets; pub mod fitstats; pub mod gpcm; pub mod grm; +pub mod ksirt; pub mod linking; pub mod lltm; pub mod marginal; @@ -26,6 +27,7 @@ pub mod rsm; pub mod rt; pub mod rt_joint; pub mod scoring; +pub mod subscores; pub mod testlet; pub mod twopl; diff --git a/crates/mlsirm-core/src/mokken.rs b/crates/mlsirm-core/src/mokken.rs index 934d96746..78fa54432 100644 --- a/crates/mlsirm-core/src/mokken.rs +++ b/crates/mlsirm-core/src/mokken.rs @@ -209,7 +209,7 @@ pub fn coef_h(x: &[i64], n_persons: usize, n_items: usize) -> Result z) = p. -fn normal_upper_quantile(p: f64) -> f64 { +pub(crate) 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); diff --git a/crates/mlsirm-core/src/subscores.rs b/crates/mlsirm-core/src/subscores.rs new file mode 100644 index 000000000..7ae13bd5a --- /dev/null +++ b/crates/mlsirm-core/src/subscores.rs @@ -0,0 +1,379 @@ +//! Haberman (2008) subscore added-value analysis via the proportional +//! reduction in mean squared error (PRMSE). +//! +//! Given a persons x items matrix of scored responses partitioned into `K` +//! disjoint, exhaustive subscales, this module computes for each subscale the +//! three Haberman estimators of the true subscore `s_t` and their PRMSEs: +//! +//! ```text +//! s_hat_s = E(s) + rho_s (s - E(s)) PRMSE_s = rho_s +//! s_hat_x = E(s) + sqrt(PRMSE_x) (sigma_t / sigma_x)(x - E(x)) +//! PRMSE_x = rho^2(s_t, x_t) rho_x +//! s_hat_sx = E(s) + beta (s - E(s)) + gamma (x - E(x)) PRMSE_sx = rho_s + tau^2 (1 - r^2) +//! ``` +//! +//! with `rho_s`/`rho_x` the Cronbach-alpha reliabilities of the subscale and +//! the total test, `r = corr(s, x)`, +//! `rho^2(s_t, x_t) = cov_k^2 / (V(s_t) V(x_t))` where `cov_k` is the row sum +//! of the true-subscore covariance matrix over the `K` subscore columns +//! (observed covariances off the diagonal, `alpha * observed variance` on the +//! diagonal; the total-score column is EXCLUDED), and +//! +//! ```text +//! tau = (sqrt(rho_x) sqrt(rho^2(s_t, x_t)) - r sqrt(rho_s)) / (1 - r^2) +//! beta = sqrt(rho_s) (sqrt(rho_s) - r tau) +//! gamma = sqrt(rho_s) tau (sigma_s / sigma_x) +//! ``` +//! +//! Added-value decisions: a subscore has added value iff +//! `PRMSE_s > PRMSE_x` (Haberman's rule); an augmented subscore has added +//! value iff `PRMSE_sx > max(PRMSE_s, PRMSE_x) + 0.01`, which is the +//! operational convention of Sinharay (2010) — the CRAN `subscore` package's +//! `CTTsub` uses a different relative rule (`0.1 * (1 - max)`) that is NOT +//! implemented here. +//! +//! # Verified sources +//! +//! Formulas were verified against (a) the Appendix of Sinharay (2010, ETS +//! RR-10-16), which reproduces the Haberman (2008) methodology, and (b) the +//! CRAN `subscore` package R source (`subscore.s.r`, `subscore.x.R`, +//! `subscore.sx.R`, `data.prep.R`) read line by line; the R code was used to +//! disambiguate Greek symbols lost in PDF extraction. Haberman (2008) itself +//! and Wainer et al. (2001) were NOT read (paywalled) and are cited only as +//! cited in Sinharay (2010). +//! +//! # Divergences from the R package (deliberate) +//! +//! - The partition is validated (every item in exactly one subscale, each +//! subscale with >= 2 items); CRAN `data.prep()` silently allows totals +//! that are not the union of the subscales. +//! - Degenerate inputs are rejected instead of propagating NaN or emitting a +//! warning: any Cronbach alpha outside `(0, 1]`, any non-positive observed +//! variance, `|corr(s_k, x)| >= 1 - 1e-12`, non-finite moments, or a +//! computed PRMSE outside `[0, 1 + 1e-9]`. +//! - `s_hat_x` uses the nonnegative root `sqrt(PRMSE_x)` exactly as the R +//! code does, even when the signed correlation form of the Sinharay +//! appendix would be negative (i.e. when `cov(s_t, x_t) < 0`, which the +//! guards do NOT rule out); this follows CRAN's convention. +//! - Missing data are not supported (the R code uses `na.rm`); v1 requires +//! complete data. +//! +//! All moments are unbiased (`n - 1`), matching R's `var`/`cov`/`cor`. +//! +//! In LLM-as-a-Judge item-quality management this decides whether per-domain +//! judge subscores carry diagnostic information beyond the overall score, or +//! whether reporting them would be statistically misleading. +//! +//! # References (APA 7th ed.) +//! +//! Haberman, S. J. (2008). When can subscores have value? *Journal of +//! Educational and Behavioral Statistics, 33*(2), 204-229. +//! https://doi.org/10.3102/1076998607302636 (as cited in Sinharay, 2010) +//! +//! Sinharay, S. (2010). *When can subscores be expected to have added value? +//! Results from operational and simulated data* (ETS Research Rep. No. +//! RR-10-16). Educational Testing Service. +//! +//! Wainer, H., Vevea, J. L., Camacho, F., Reeve, B. B., Rosa, K., & Nelson, +//! L. (2001). Augmented scores — "borrowing strength" to compute scores based +//! on small numbers of items. In D. Thissen & H. Wainer (Eds.), *Test +//! scoring* (pp. 343-387). Lawrence Erlbaum. (as cited in Sinharay, 2010) + +/// Result of the Haberman subscore added-value analysis. All vectors are +/// indexed by subscale `k = 0..K`; person-level estimator matrices are +/// `n_persons x K` in row-major nested `Vec`s. +#[derive(Debug, Clone)] +pub struct SubscoreResult { + /// Cronbach alpha of each subscale ( = PRMSE_s). + pub alpha: Vec, + /// Cronbach alpha of the total test. + pub alpha_total: f64, + /// `(K+1) x (K+1)` correlation matrix of `(s_1..s_K, x)`, total last. + pub corr: Vec>, + /// `K x K` disattenuated subscore correlations + /// `corr_kl / sqrt(alpha_k alpha_l)`; diagonal is NaN. + pub disattenuated_corr: Vec>, + /// PRMSE of the observed-subscore estimator ( = subscale reliability). + pub prmse_s: Vec, + /// PRMSE of the total-score estimator `rho^2(s_t, x_t) rho_x`. + pub prmse_x: Vec, + /// PRMSE of the augmented estimator `rho_s + tau^2 (1 - r^2)`. + pub prmse_sx: Vec, + /// Regression helpers for the augmented estimator. + pub tau: Vec, + pub beta: Vec, + pub gamma: Vec, + /// `PRMSE_s > PRMSE_x` (Haberman's added-value rule). + pub added_value_s: Vec, + /// `PRMSE_sx > max(PRMSE_s, PRMSE_x) + 0.01` (Sinharay 2010 convention). + pub added_value_sx: Vec, + /// Observed subscores `s_k` per person (`n x K`). + pub observed: Vec>, + /// Observed total score per person. + pub total: Vec, + /// Estimated true subscores from the observed subscore (`n x K`). + pub subscore_s: Vec>, + /// Estimated true subscores from the observed total (`n x K`). + pub subscore_x: Vec>, + /// Estimated true subscores from both (`n x K`). + pub subscore_sx: Vec>, +} + +/// Unbiased sample variance (`n - 1`); caller guarantees `v.len() >= 2`. +fn var(v: &[f64]) -> f64 { + let n = v.len() as f64; + let m = v.iter().sum::() / n; + v.iter().map(|&a| (a - m) * (a - m)).sum::() / (n - 1.0) +} + +/// Unbiased sample covariance (`n - 1`). +fn cov(a: &[f64], b: &[f64]) -> f64 { + let n = a.len() as f64; + let (ma, mb) = ( + a.iter().sum::() / n, + b.iter().sum::() / n, + ); + a.iter() + .zip(b) + .map(|(&x, &y)| (x - ma) * (y - mb)) + .sum::() + / (n - 1.0) +} + +/// Cronbach alpha `m/(m-1) (1 - sum_j var(y_j) / var(sum_j y_j))` over the +/// given item columns. Same unbiased convention in numerator and denominator. +fn cronbach_alpha(x: &[Vec], items: &[usize]) -> f64 { + let m = items.len() as f64; + let totals: Vec = x + .iter() + .map(|row| items.iter().map(|&j| row[j]).sum()) + .collect(); + let item_var_sum: f64 = items + .iter() + .map(|&j| { + let col: Vec = x.iter().map(|row| row[j]).collect(); + var(&col) + }) + .sum(); + let tv = var(&totals); + m / (m - 1.0) * (1.0 - item_var_sum / tv) +} + +/// Haberman (2008) subscore added-value analysis (see module docs). +/// +/// `x` is a complete `n_persons x n_items` matrix of finite scored responses; +/// `groups[j]` in `0..K` assigns item `j` to a subscale. The partition must be +/// exhaustive with at least 2 items per subscale, `n_persons >= 3`, and +/// `K >= 2`. Returns an error on any validation or degeneracy failure (see +/// the module-level guard list). +pub fn subscores(x: &[Vec], groups: &[usize]) -> Result { + let n = x.len(); + if n < 3 { + return Err("subscores requires at least 3 persons".into()); + } + let n_items = x[0].len(); + if groups.len() != n_items { + return Err("groups must assign every item to a subscale".into()); + } + for row in x { + if row.len() != n_items { + return Err("ragged response matrix".into()); + } + if row.iter().any(|v| !v.is_finite()) { + return Err("responses must be complete and finite".into()); + } + } + let k_count = match groups.iter().max() { + Some(&g) => g + 1, + None => return Err("groups must assign every item to a subscale".into()), + }; + if k_count < 2 { + return Err("at least 2 subscales are required".into()); + } + // Bound BEFORE allocating items_of: a hostile sparse index (e.g. 10^9) + // would otherwise drive a huge allocation. K <= n_items/2 given >= 2 + // items per subscale. + if k_count > n_items / 2 { + return Err("subscale indices must be dense: 0..K with K <= n_items / 2".into()); + } + let mut items_of: Vec> = vec![Vec::new(); k_count]; + for (j, &g) in groups.iter().enumerate() { + items_of[g].push(j); + } + if items_of.iter().any(|v| v.len() < 2) { + return Err("every subscale needs at least 2 items".into()); + } + + // Observed subscores and total (the partition makes x = sum_k s_k). + let observed: Vec> = x + .iter() + .map(|row| { + items_of + .iter() + .map(|items| items.iter().map(|&j| row[j]).sum()) + .collect() + }) + .collect(); + let total: Vec = observed.iter().map(|s| s.iter().sum()).collect(); + + // Columns of the (K+1)-variate score matrix, total last. + let mut cols: Vec> = (0..k_count) + .map(|k| observed.iter().map(|s| s[k]).collect()) + .collect(); + cols.push(total.clone()); + let kk = k_count + 1; + + let all_items: Vec = (0..n_items).collect(); + let mut alpha: Vec = items_of + .iter() + .map(|items| cronbach_alpha(x, items)) + .collect(); + let alpha_total = cronbach_alpha(x, &all_items); + alpha.push(alpha_total); + for (i, &a) in alpha.iter().enumerate() { + if !a.is_finite() || a <= 0.0 || a > 1.0 { + return Err(format!( + "Cronbach alpha of {} is {a:.6}, outside (0, 1]; the Haberman \ + analysis is undefined (zero variance or negatively \ + correlated items)", + if i < k_count { "a subscale" } else { "the total test" } + )); + } + } + + let var_obs: Vec = cols.iter().map(|c| var(c)).collect(); + if var_obs.iter().any(|&v| !(v.is_finite() && v > 0.0)) { + return Err("zero-variance subscore or total score".into()); + } + let mut c_obs = vec![vec![0.0f64; kk]; kk]; + for a in 0..kk { + for b in a..kk { + let v = cov(&cols[a], &cols[b]); + c_obs[a][b] = v; + c_obs[b][a] = v; + } + } + let corr: Vec> = (0..kk) + .map(|a| { + (0..kk) + .map(|b| c_obs[a][b] / (var_obs[a] * var_obs[b]).sqrt()) + .collect() + }) + .collect(); + for k in 0..k_count { + if corr[k][k_count].abs() >= 1.0 - 1e-12 { + return Err( + "a subscore is (numerically) collinear with the total score; \ + the augmented regression is undefined" + .into(), + ); + } + } + + // True-score covariance matrix C_T: observed off-diagonals, diagonal + // alpha_k * V(s_k). cov_k sums row k over the K subscore columns ONLY + // (the total column is excluded). + let var_true: Vec = (0..kk).map(|i| var_obs[i] * alpha[i]).collect(); + let cov_rowsum: Vec = (0..k_count) + .map(|k| { + (0..k_count) + .map(|l| if l == k { var_true[k] } else { c_obs[k][l] }) + .sum() + }) + .collect(); + + let mut prmse_s = Vec::with_capacity(k_count); + let mut prmse_x = Vec::with_capacity(k_count); + let mut prmse_sx = Vec::with_capacity(k_count); + let (mut tau, mut beta, mut gamma) = ( + Vec::with_capacity(k_count), + Vec::with_capacity(k_count), + Vec::with_capacity(k_count), + ); + for k in 0..k_count { + let r_stxt = cov_rowsum[k] * cov_rowsum[k] / (var_true[k] * var_true[k_count]); + let r = corr[k][k_count]; + let t = (alpha_total.sqrt() * r_stxt.sqrt() - r * alpha[k].sqrt()) / (1.0 - r * r); + let ps = alpha[k]; + let px = r_stxt * alpha_total; + let psx = alpha[k] + t * t * (1.0 - r * r); + for (name, v) in [("PRMSE_s", ps), ("PRMSE_x", px), ("PRMSE_sx", psx)] { + if !v.is_finite() || v < 0.0 || v > 1.0 + 1e-9 { + return Err(format!( + "computed {name} = {v:.6} outside [0, 1]; the sample \ + moments are inconsistent with the CTT assumptions" + )); + } + } + prmse_s.push(ps); + prmse_x.push(px); + prmse_sx.push(psx); + tau.push(t); + beta.push(alpha[k].sqrt() * (alpha[k].sqrt() - r * t)); + gamma.push(alpha[k].sqrt() * t * (var_obs[k].sqrt() / var_obs[k_count].sqrt())); + } + + let disattenuated_corr: Vec> = (0..k_count) + .map(|a| { + (0..k_count) + .map(|b| { + if a == b { + f64::NAN + } else { + corr[a][b] / (alpha[a] * alpha[b]).sqrt() + } + }) + .collect() + }) + .collect(); + + let mean_s: Vec = (0..k_count) + .map(|k| cols[k].iter().sum::() / n as f64) + .collect(); + let mean_x = total.iter().sum::() / n as f64; + let sd_x = var_obs[k_count].sqrt(); + + let mut subscore_s = vec![vec![0.0f64; k_count]; n]; + let mut subscore_x = vec![vec![0.0f64; k_count]; n]; + let mut subscore_sx = vec![vec![0.0f64; k_count]; n]; + for p in 0..n { + for k in 0..k_count { + let ds = observed[p][k] - mean_s[k]; + let dx = total[p] - mean_x; + subscore_s[p][k] = mean_s[k] + alpha[k] * ds; + subscore_x[p][k] = + mean_s[k] + prmse_x[k].sqrt() * (var_true[k].sqrt() / sd_x) * dx; + subscore_sx[p][k] = mean_s[k] + beta[k] * ds + gamma[k] * dx; + } + } + + let added_value_s: Vec = (0..k_count).map(|k| prmse_s[k] > prmse_x[k]).collect(); + let added_value_sx: Vec = (0..k_count) + .map(|k| prmse_sx[k] > prmse_s[k].max(prmse_x[k]) + 0.01) + .collect(); + + alpha.pop(); + Ok(SubscoreResult { + alpha, + alpha_total, + corr, + disattenuated_corr, + prmse_s, + prmse_x, + prmse_sx, + tau, + beta, + gamma, + added_value_s, + added_value_sx, + observed, + total, + subscore_s, + subscore_x, + subscore_sx, + }) +} + +#[cfg(test)] +#[path = "../../../tests/unit/subscores_tests.rs"] +mod tests; diff --git a/python/fast_mlsirm/__init__.py b/python/fast_mlsirm/__init__.py index b41f94d50..8af312542 100644 --- a/python/fast_mlsirm/__init__.py +++ b/python/fast_mlsirm/__init__.py @@ -33,6 +33,11 @@ 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 .ksirt import ksirt_analysis as ksirt_analysis, KsirtResult as KsirtResult +from .subscores import ( + subscore_analysis as subscore_analysis, + SubscoreResult as SubscoreResult, +) 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 @@ -149,6 +154,10 @@ "FacetsFit", "mokken_analysis", "MokkenResult", + "ksirt_analysis", + "KsirtResult", + "subscore_analysis", + "SubscoreResult", "fit_mixed_items", "MixedFormatFit", "MixedItemParameters", diff --git a/python/fast_mlsirm/ksirt.py b/python/fast_mlsirm/ksirt.py new file mode 100644 index 000000000..537138b75 --- /dev/null +++ b/python/fast_mlsirm/ksirt.py @@ -0,0 +1,125 @@ +"""Kernel-smoothing nonparametric IRT: option characteristic curves by +Nadaraya-Watson regression on rank-based ordinal ability estimates (Ramsay, +1991, as cited in Mazza et al., 2014). 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 + + +@dataclass +class KsirtResult: + """Kernel-smoothed option characteristic curves. + + ``theta`` are the rank-based ordinal ability estimates + ``Phi^-1(rank(total_i)/(n+1))`` in subject order; ``grid`` the + equally-spaced evaluation points; ``bandwidth`` the per-item bandwidths + used. ``options[j]`` lists item ``j``'s distinct observed scores + (ascending), ``occ[j]`` is the matching ``m_j x len(grid)`` option + characteristic curve matrix, ``expected[j]`` the expected item score + curve, and ``expected_total`` their sum over items.""" + + theta: np.ndarray + grid: np.ndarray + bandwidth: np.ndarray + options: list[np.ndarray] + occ: list[np.ndarray] + expected: list[np.ndarray] + expected_total: np.ndarray + + +def ksirt_analysis( + responses: np.ndarray, + kernel: str = "gaussian", + nevalpoints: int = 51, + bandwidth: np.ndarray | None = None, +) -> KsirtResult: + """Kernel smoothing of option characteristic curves (compute in Rust; + Ramsay, 1991, as cited in Mazza et al., 2014). + + Estimates each item's option characteristic curves by Nadaraya-Watson + kernel regression of the option indicators on ordinal ability estimates + ``Phi^-1(rank(total score)/(n+1))`` (ties broken by subject order), + evaluated on an equally-spaced grid from ``Phi^-1(1/(n+1))`` to + ``Phi^-1(n/(n+1))``. The default bandwidth is Silverman's rule + ``1.06 * n^(-1/5)`` on the standard-normal ability metric. Formulas + follow Mazza et al. (2014, Sections 2-2.3) and the KernSmoothIRT R + package source (both read); Ramsay (1991) itself is cited only through + Mazza et al. (2014). Standard errors and cross-validation bandwidths + are deliberately not implemented (see the Rust module docs). + + In LLM-as-a-Judge item-quality management, nonparametric OCCs reveal + non-monotone or poorly discriminating evaluation items without assuming + a parametric response model. + + ``responses`` is a complete ``persons x items`` array of pre-scored + numeric responses; each column's distinct values form that item's + options. ``kernel`` is ``"gaussian"``, ``"quadratic"``, or + ``"uniform"``. ``bandwidth`` optionally gives one positive value per + item. + + References (APA 7th ed.): + Mazza, A., Punzo, A., & McGuire, B. (2014). KernSmoothIRT: An R + package for kernel smoothing in item response theory. *Journal + of Statistical Software, 58*(6), 1-34. + https://doi.org/10.18637/jss.v058.i06 + Ramsay, J. O. (1991). Kernel smoothing approaches to nonparametric + item characteristic curve estimation. *Psychometrika, 56*(4), + 611-630. https://doi.org/10.1007/BF02294494 (as cited in Mazza + et al., 2014) + """ + from .fitstats import _core_module + + core = _core_module() + if core is None or not hasattr(core, "ksirt_occ"): + raise RuntimeError("ksirt_analysis requires the compiled Rust core") + + if kernel not in ("gaussian", "quadratic", "uniform"): + raise ValueError("kernel must be gaussian, quadratic, or uniform") + nevalpoints = int(nevalpoints) + if nevalpoints < 2: + raise ValueError("nevalpoints must be at least 2") + if nevalpoints > 100_000: + # trust boundary: nevalpoints drives Rust-side allocations + raise ValueError("nevalpoints must be at most 100000") + + 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 n_persons < 2 or n_items < 1: + raise ValueError("responses needs at least 2 persons and 1 item") + if not np.all(np.isfinite(y)): + raise ValueError("responses must be complete (no missing values)") + + bw = None + if bandwidth is not None: + bw_arr = np.asarray(bandwidth, dtype=np.float64).reshape(-1) + if bw_arr.shape[0] != n_items: + raise ValueError("bandwidth must supply one value per item") + if not np.all(np.isfinite(bw_arr)) or np.any(bw_arr <= 0.0): + raise ValueError("bandwidths must be finite and positive") + bw = [float(v) for v in bw_arr] + + res = core.ksirt_occ( + y.reshape(-1), int(n_persons), int(n_items), kernel, nevalpoints, bw + ) + grid = np.asarray(res["grid"], dtype=np.float64) + q = grid.shape[0] + options = [np.asarray(o, dtype=np.float64) for o in res["options"]] + occ = [ + np.asarray(flat, dtype=np.float64).reshape(len(opts), q) + for flat, opts in zip(res["occ"], options) + ] + return KsirtResult( + theta=np.asarray(res["theta"], dtype=np.float64), + grid=grid, + bandwidth=np.asarray(res["bandwidth"], dtype=np.float64), + options=options, + occ=occ, + expected=[np.asarray(e, dtype=np.float64) for e in res["expected"]], + expected_total=np.asarray(res["expected_total"], dtype=np.float64), + ) diff --git a/python/fast_mlsirm/subscores.py b/python/fast_mlsirm/subscores.py new file mode 100644 index 000000000..23595a137 --- /dev/null +++ b/python/fast_mlsirm/subscores.py @@ -0,0 +1,154 @@ +"""Haberman subscore added-value analysis via proportional reduction in mean +squared error (PRMSE; Haberman, 2008, as cited in Sinharay, 2010). 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 + + +@dataclass +class SubscoreResult: + """Haberman subscore added-value analysis for ``K`` subscales. + + ``alpha`` are the per-subscale Cronbach alphas ( = ``prmse_s``); + ``alpha_total`` the total-test alpha. ``corr`` is the ``(K+1) x (K+1)`` + correlation matrix of the observed subscores with the total score last; + ``disattenuated_corr`` the ``K x K`` disattenuated subscore correlations + (NaN diagonal). ``prmse_s``/``prmse_x``/``prmse_sx`` are the PRMSEs of + predicting the true subscore from the observed subscore, the observed + total, and both; ``tau``/``beta``/``gamma`` the augmented-regression + weights. ``added_value_s`` is Haberman's rule ``PRMSE_s > PRMSE_x``; + ``added_value_sx`` uses Sinharay's (2010) ``+ 0.01`` margin. + ``observed`` (``n x K``), ``total`` (``n``), and the three estimator + matrices ``subscore_s``/``subscore_x``/``subscore_sx`` (each ``n x K``) + give per-person scores.""" + + alpha: np.ndarray + alpha_total: float + corr: np.ndarray + disattenuated_corr: np.ndarray + prmse_s: np.ndarray + prmse_x: np.ndarray + prmse_sx: np.ndarray + tau: np.ndarray + beta: np.ndarray + gamma: np.ndarray + added_value_s: np.ndarray + added_value_sx: np.ndarray + observed: np.ndarray + total: np.ndarray + subscore_s: np.ndarray + subscore_x: np.ndarray + subscore_sx: np.ndarray + + +def subscore_analysis( + responses: np.ndarray, + groups: np.ndarray, +) -> SubscoreResult: + """Haberman subscore added-value analysis (compute in Rust; Haberman, + 2008, as cited in Sinharay, 2010). + + Decides, for each subscale of a test, whether reporting its subscore adds + value over reporting the total score alone, by comparing the PRMSEs of + three classical-test-theory estimators of the true subscore (from the + observed subscore, from the observed total, and from both jointly). + Formulas follow the Appendix of Sinharay (2010) and the CRAN ``subscore`` + package R source (both read); Haberman (2008) and Wainer et al. (2001) + are cited only through Sinharay (2010). Degenerate samples (any Cronbach + alpha outside ``(0, 1]``, zero-variance scores, a subscore collinear with + the total) are rejected with ``ValueError`` rather than propagating NaN. + + In LLM-as-a-Judge item-quality management this decides whether + per-domain judge subscores carry diagnostic information beyond the + overall score, or whether reporting them would be statistically + misleading. + + ``responses`` is a complete ``persons x items`` array of scored + responses (``n >= 3``). ``groups`` assigns each item an integer subscale + index in ``0..K`` (``K >= 2``, every subscale with at least 2 items, + partition exhaustive by construction). + + References (APA 7th ed.): + Haberman, S. J. (2008). When can subscores have value? *Journal of + Educational and Behavioral Statistics, 33*(2), 204-229. + https://doi.org/10.3102/1076998607302636 (as cited in Sinharay, + 2010) + Sinharay, S. (2010). *When can subscores be expected to have added + value? Results from operational and simulated data* (ETS + Research Rep. No. RR-10-16). Educational Testing Service. + Wainer, H., Vevea, J. L., Camacho, F., Reeve, B. B., Rosa, K., & + Nelson, L. (2001). Augmented scores — "borrowing strength" to + compute scores based on small numbers of items. In D. Thissen & + H. Wainer (Eds.), *Test scoring* (pp. 343-387). Lawrence + Erlbaum. (as cited in Sinharay, 2010) + """ + from .fitstats import _core_module + + core = _core_module() + if core is None or not hasattr(core, "subscore_analysis"): + raise RuntimeError("subscore_analysis requires the compiled Rust core") + + 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 n_persons < 3 or n_items < 4: + # K >= 2 subscales with >= 2 items each needs at least 4 items; + # rejecting here keeps degenerate shapes (e.g. huge zero-column + # arrays) from crossing the Rust boundary at all. + raise ValueError("responses needs at least 3 persons and 4 items") + if not np.all(np.isfinite(y)): + raise ValueError("responses must be complete (no missing values)") + + g = np.asarray(groups).reshape(-1) + if g.shape[0] != n_items: + raise ValueError("groups must assign one subscale index per item") + if not np.issubdtype(g.dtype, np.integer): + gf = np.asarray(groups, dtype=np.float64).reshape(-1) + if not np.all(np.isfinite(gf)) or np.any(gf != np.round(gf)): + raise ValueError("groups must be integer subscale indices") + g = gf.astype(np.int64) + if np.any(g < 0): + raise ValueError("groups must be nonnegative subscale indices") + if np.any(g >= n_items): + # trust boundary: the subscale count drives Rust-side allocations + raise ValueError("groups indices must be < n_items") + + res = core.subscore_analysis( + y.reshape(-1), int(n_persons), int(n_items), [int(v) for v in g] + ) + k = len(res["alpha"]) + return SubscoreResult( + alpha=np.asarray(res["alpha"], dtype=np.float64), + alpha_total=float(res["alpha_total"]), + corr=np.asarray(res["corr"], dtype=np.float64).reshape(k + 1, k + 1), + disattenuated_corr=np.asarray( + res["disattenuated_corr"], dtype=np.float64 + ).reshape(k, k), + prmse_s=np.asarray(res["prmse_s"], dtype=np.float64), + prmse_x=np.asarray(res["prmse_x"], dtype=np.float64), + prmse_sx=np.asarray(res["prmse_sx"], dtype=np.float64), + tau=np.asarray(res["tau"], dtype=np.float64), + beta=np.asarray(res["beta"], dtype=np.float64), + gamma=np.asarray(res["gamma"], dtype=np.float64), + added_value_s=np.asarray(res["added_value_s"], dtype=bool), + added_value_sx=np.asarray(res["added_value_sx"], dtype=bool), + observed=np.asarray(res["observed"], dtype=np.float64).reshape( + n_persons, k + ), + total=np.asarray(res["total"], dtype=np.float64), + subscore_s=np.asarray(res["subscore_s"], dtype=np.float64).reshape( + n_persons, k + ), + subscore_x=np.asarray(res["subscore_x"], dtype=np.float64).reshape( + n_persons, k + ), + subscore_sx=np.asarray(res["subscore_sx"], dtype=np.float64).reshape( + n_persons, k + ), + ) diff --git a/tests/test_paper_features.py b/tests/test_paper_features.py index 9835d9c0f..2a86176f2 100644 --- a/tests/test_paper_features.py +++ b/tests/test_paper_features.py @@ -5083,4 +5083,194 @@ def test_mokken_analysis_rejects_malformed_input(): cst = ok.copy() cst[:, 0] = 1.0 with pytest.raises(ValueError, match="zero variance"): - mokken_analysis(cst) \ No newline at end of file + mokken_analysis(cst) + +def test_ksirt_analysis_fixture_exact_values(): + """Kernel-smoothing IRT (Ramsay, 1991, as cited in Mazza et al., 2014): + every assert reads crate outputs (theta/grid/bandwidth/occ/expected/ + expected_total from ksirt_occ via the wrapper) against independently + derived constants. theta must equal qnorm literals (kills rank-map + mutants), grid endpoints qnorm(1/5)/qnorm(4/5) (kills grid mutants), + bandwidth the Silverman literal 1.06*4^(-1/5) (kills bandwidth mutants), + and a huge bandwidth must flatten every OCC to the marginal proportion + 0.5 (kills NW-normalization mutants).""" + import numpy as np + import pytest + from fast_mlsirm import KsirtResult, ksirt_analysis + from fast_mlsirm.fitstats import _core_module + + core = _core_module() + if core is None or not hasattr(core, "ksirt_occ"): + pytest.skip("compiled core built without ksirt_occ") + + x = np.array([[0.0], [1.0], [0.0], [1.0]]) + q02, q04 = -0.8416212335729143, -0.2533471031357997 # qnorm(.2), qnorm(.4) + r = ksirt_analysis(x, nevalpoints=5) + assert isinstance(r, KsirtResult) + # ties.method="first": totals [0,1,0,1] -> ranks [1,3,2,4] + assert np.allclose(r.theta, [q02, -q04, q04, -q02], atol=1e-8) + assert abs(r.grid[0] - q02) < 1e-8 and abs(r.grid[-1] + q02) < 1e-8 + assert np.allclose(np.diff(r.grid), np.diff(r.grid)[0], atol=1e-12) + assert abs(r.bandwidth[0] - 1.06 * 4.0 ** (-0.2)) < 1e-12 + assert len(r.options) == 1 and np.array_equal(r.options[0], [0.0, 1.0]) + assert r.occ[0].shape == (2, 5) + # huge bandwidth -> equal NW weights -> OCC = marginal proportions + flat = ksirt_analysis(x, nevalpoints=5, bandwidth=np.array([1e9])) + assert np.allclose(flat.occ[0], 0.5, atol=1e-12) + assert np.allclose(flat.expected[0], 0.5, atol=1e-12) + assert np.allclose(flat.expected_total, 0.5, atol=1e-12) + + +def test_ksirt_analysis_recovery_and_rejects_malformed_input(): + """Monotone-recovery smoke on a simulated 2PL item (crate expected_total + must increase over the grid and crate theta must correlate with true + ability; kills over-collapse/ordering mutants) plus wrapper guard + checks (each raise exercises a validation branch in front of the + crate).""" + import numpy as np + import pytest + from fast_mlsirm import ksirt_analysis + from fast_mlsirm.fitstats import _core_module + + core = _core_module() + if core is None or not hasattr(core, "ksirt_occ"): + pytest.skip("compiled core built without ksirt_occ") + + rng = np.random.default_rng(1991) + n = 800 + t = rng.standard_normal(n) + bs = np.array([-1.0, -0.3, 0.4, 1.1]) + x = (rng.random((n, 4)) < 1.0 / (1.0 + np.exp(-1.5 * (t[:, None] - bs)))).astype(float) + r = ksirt_analysis(x, kernel="gaussian") + d = np.diff(r.expected_total) + assert np.all(d > -1e-9) # non-decreasing everywhere + assert np.all(d[10:-10] > 0.0) # strictly increasing in the interior + assert r.expected_total[-1] - r.expected_total[0] > 2.0 + assert np.corrcoef(r.theta, t)[0, 1] > 0.7 + + ok = x[:8] + with pytest.raises(ValueError, match="2-D"): + ksirt_analysis(ok.reshape(-1)) + with pytest.raises(ValueError, match="complete"): + bad = ok.copy() + bad[0, 0] = np.nan + ksirt_analysis(bad) + with pytest.raises(ValueError, match="kernel"): + ksirt_analysis(ok, kernel="triangular") + with pytest.raises(ValueError, match="nevalpoints"): + ksirt_analysis(ok, nevalpoints=1) + with pytest.raises(ValueError, match="nevalpoints"): + ksirt_analysis(ok, nevalpoints=10**11) # allocation bound + with pytest.raises(ValueError, match="one value per item"): + ksirt_analysis(ok, bandwidth=np.array([0.5])) + with pytest.raises(ValueError, match="positive"): + ksirt_analysis(ok, bandwidth=np.array([0.5, -0.1, 0.5, 0.5])) + + +def test_subscore_analysis_matches_independent_reference(): + """Haberman (2008, as cited in Sinharay, 2010) subscore added-value + analysis: every assert reads the SubscoreResult returned by the crate via + the PyO3 binding, pinned against literals from an independent NumPy + transcription of the CRAN subscore R semantics (ddof=1). The fixture is + asymmetric across subscales (added_value_sx = [False, True]), so + decision-rule and PRMSE mutants (rowsum including the total column, + missing true-variance diagonal, dropped +0.01 margin) all fail here.""" + import numpy as np + import pytest + from fast_mlsirm import SubscoreResult, subscore_analysis + from fast_mlsirm.fitstats import _core_module + + core = _core_module() + if core is None or not hasattr(core, "subscore_analysis"): + pytest.skip("compiled core built without subscore_analysis") + + x = np.array( + [ + [0, 0, 0, 1], [0, 1, 0, 0], [1, 1, 0, 1], [0, 0, 1, 0], + [1, 1, 1, 1], [1, 0, 1, 1], [0, 0, 0, 0], [1, 1, 1, 0], + [0, 1, 1, 1], [1, 1, 0, 0], + ], + dtype=float, + ) + res = subscore_analysis(x, [0, 0, 1, 1]) + assert isinstance(res, SubscoreResult) + np.testing.assert_allclose( + res.alpha, [0.5797101449275359, 0.33333333333333326], atol=1e-10 + ) + assert abs(res.alpha_total - 0.44742729306487705) < 1e-12 + np.testing.assert_allclose( + res.prmse_x, [0.41946308724832204, 0.3020134228187919], atol=1e-10 + ) + np.testing.assert_allclose( + res.prmse_sx, [0.5872524752475244, 0.3663366336633662], atol=1e-10 + ) + np.testing.assert_allclose( + res.tau, [0.1385414635898375, 0.27024443691002326], atol=1e-10 + ) + np.testing.assert_allclose( + res.beta, [0.49752475247524747, 0.21782178217821788], atol=1e-10 + ) + np.testing.assert_allclose( + res.gamma, [0.07178217821782161, 0.09900990099009896], atol=1e-10 + ) + assert res.added_value_s.tolist() == [True, True] + assert res.added_value_sx.tolist() == [False, True] + assert res.corr.shape == (3, 3) + assert abs(res.corr[0, 2] - 0.7791290756515445) < 1e-10 + assert np.isnan(res.disattenuated_corr[0, 0]) + assert abs(res.disattenuated_corr[0, 1] - 0.35355339059327384) < 1e-10 + # person-0 estimates from all three estimators + np.testing.assert_allclose( + res.subscore_s[0], [0.4623188405797105, 1.0], atol=1e-10 + ) + np.testing.assert_allclose( + res.subscore_x[0], [0.7308724832214768, 0.778523489932886], atol=1e-10 + ) + np.testing.assert_allclose( + res.subscore_sx[0], [0.47376237623762407, 0.8910891089108911], atol=1e-10 + ) + assert res.observed.shape == (10, 2) and res.total.shape == (10,) + + +def test_subscore_analysis_rejects_degenerate_inputs(): + """Guard behavior (spec-review mandated): duplicated subscales make a + subscore collinear with the total (tau divides by 1 - r^2 = 0) and must + raise, as must negative within-subscale alpha, bad partitions, and + incomplete data. Every assert reads the ValueError raised by the crate + through the binding; removing any Rust-side guard turns these into NaN + results or panics and fails the test.""" + import numpy as np + import pytest + from fast_mlsirm import subscore_analysis + + from fast_mlsirm.fitstats import _core_module + + core = _core_module() + if core is None or not hasattr(core, "subscore_analysis"): + pytest.skip("compiled core built without subscore_analysis") + + base = np.array( + [ + [0, 0, 0, 1], [0, 1, 0, 0], [1, 1, 0, 1], [0, 0, 1, 0], + [1, 1, 1, 1], [1, 0, 1, 1], [0, 0, 0, 0], [1, 1, 1, 0], + [0, 1, 1, 1], [1, 1, 0, 0], + ], + dtype=float, + ) + dup = base[:, [0, 1, 0, 1]] + with pytest.raises(ValueError, match="collinear"): + subscore_analysis(dup, [0, 0, 1, 1]) + rev = base.copy() + rev[:, 1] = 1.0 - rev[:, 0] # reversed item -> negative alpha + with pytest.raises(ValueError, match="outside"): + subscore_analysis(rev, [0, 0, 1, 1]) + with pytest.raises(ValueError): + subscore_analysis(base, [0, 0, 1]) # groups length mismatch + with pytest.raises(ValueError): + subscore_analysis(base, [0, 0, 0, 0]) # single subscale + with pytest.raises(ValueError): + subscore_analysis(base, [0, 1, 1, 1]) # singleton subscale + nanx = base.copy() + nanx[0, 0] = np.nan + with pytest.raises(ValueError, match="complete"): + subscore_analysis(nanx, [0, 0, 1, 1]) diff --git a/tests/unit/ksirt_tests.rs b/tests/unit/ksirt_tests.rs new file mode 100644 index 000000000..a761ad962 --- /dev/null +++ b/tests/unit/ksirt_tests.rs @@ -0,0 +1,367 @@ +//! Tests for kernel-smoothing nonparametric IRT (`mlsirm_core::ksirt`). +//! +//! Every assert reads values returned by the crate (`KsirtResult` fields). +//! Each test names the crate value it reads and a mutant it kills. +//! +//! Known unkillable mutations (documented per the test discipline): all +//! three supported kernels are symmetric, so flipping the sign of the +//! Nadaraya-Watson argument `(grid - theta)/h` is an identity; likewise a +//! pure multiplicative kernel constant cancels in the NW normalization. +//! No anchor can exist for either — this is a property of the model, not a +//! test gap. + +use super::{ksirt, KsirtKernel}; + +/// qnorm anchors (R `qnorm`, 16 digits) used as independent fixture +/// constants so expected values never route through the crate's quantile. +const QN_1_5: f64 = -0.8416212335729143; // Phi^-1(0.2) +const QN_2_5: f64 = -0.2533471031357997; // Phi^-1(0.4) +const QN_3_5: f64 = 0.2533471031357997; // Phi^-1(0.6) +const QN_4_5: f64 = 0.8416212335729143; // Phi^-1(0.8) + +/// n=4, k=1 binary fixture: subjects 1 and 3 score 1, subjects 0 and 2 +/// score 0. Totals [0,1,0,1] -> first-occurrence ranks [1,3,2,4] -> +/// theta [Phi^-1(.2), Phi^-1(.6), Phi^-1(.4), Phi^-1(.8)]. +fn fixture_theta() -> [f64; 4] { + [QN_1_5, QN_3_5, QN_2_5, QN_4_5] +} + +fn fixture_x() -> Vec> { + vec![vec![0.0], vec![1.0], vec![0.0], vec![1.0]] +} + +/// Reads: `items[0].occ` (both option rows, all grid points). +/// Expected values recomputed here from hard-coded qnorm literals and the +/// Gaussian NW formula (never from crate outputs), so the assert compares +/// the crate curve to an independent hand derivation. +/// Kills: subject/theta index misalignment (swapping two subjects' theta), +/// wrong weight denominator, gaussian<->quadratic kernel swap, and rank +/// mis-ordering (theta enters the expected values). +#[test] +fn hand_fixture_occ_exact() { + let res = ksirt(&fixture_x(), KsirtKernel::Gaussian, 3, Some(&[0.5])).unwrap(); + let theta = fixture_theta(); + let grid = [QN_1_5, 0.0, QN_4_5]; + for (s, &g) in grid.iter().enumerate() { + let kw: Vec = theta + .iter() + .map(|&t| (-0.5 * ((g - t) / 0.5).powi(2)).exp()) + .collect(); + let denom: f64 = kw.iter().sum(); + let exp_p1 = (kw[1] + kw[3]) / denom; // subjects with score 1 + assert!( + (res.items[0].occ[1][s] - exp_p1).abs() < 1e-9, + "occ[1][{s}] = {} expected {exp_p1}", + res.items[0].occ[1][s] + ); + assert!((res.items[0].occ[0][s] - (1.0 - exp_p1)).abs() < 1e-9); + } + // asymmetry anchor: the curve is not flat, low grid point is low + assert!(res.items[0].occ[1][0] < 0.10); + assert!(res.items[0].occ[1][2] > 0.90); +} + +/// Reads: `result.theta`. +/// Kills: `ties.method="first"` violations (tied totals must rank in +/// subject order), and the n+1 -> n denominator mutation (rank 4 of 4 +/// would hit Phi^-1(1.0) = infinity instead of Phi^-1(0.8)). +#[test] +fn theta_rank_ties_first() { + // totals: [2, 5, 2, 7] -> ties between subjects 0 and 2 broken by + // original order: ranks [1, 3, 2, 4] + let x = vec![ + vec![1.0, 1.0, 0.0], + vec![2.0, 2.0, 1.0], + vec![0.0, 1.0, 1.0], + vec![3.0, 3.0, 1.0], + ]; + let res = ksirt(&x, KsirtKernel::Gaussian, 5, None).unwrap(); + let expected = [QN_1_5, QN_3_5, QN_2_5, QN_4_5]; + for i in 0..4 { + assert!( + (res.theta[i] - expected[i]).abs() < 1e-8, + "theta[{i}] = {} expected {}", + res.theta[i], + expected[i] + ); + } +} + +/// Reads: `result.grid`. +/// Kills: wrong endpoints (1/(n+1) vs 1/n), wrong point count, and +/// non-uniform spacing mutations. +#[test] +fn grid_endpoints_and_spacing() { + let x = fixture_x(); + let res = ksirt(&x, KsirtKernel::Gaussian, 5, None).unwrap(); + assert_eq!(res.grid.len(), 5); + assert!((res.grid[0] - QN_1_5).abs() < 1e-8); + assert!((res.grid[4] - QN_4_5).abs() < 1e-8); + let step = (QN_4_5 - QN_1_5) / 4.0; + for s in 1..5 { + // 1e-8: the crate grid endpoints come from the Acklam quantile + // (|error| < 1.15e-9), which propagates into the spacing. + assert!((res.grid[s] - res.grid[s - 1] - step).abs() < 1e-8); + } +} + +/// Reads: `result.bandwidth`. +/// Kills: mutations of the Silverman constant (1.06) or exponent (-1/5), +/// verified against 1.06 * 100^(-0.2) = 0.4219936007867... computed by hand. +#[test] +fn silverman_bandwidth_value() { + let mut x = Vec::new(); + for i in 0..100 { + x.push(vec![(i % 2) as f64, ((i / 2) % 2) as f64]); + } + let res = ksirt(&x, KsirtKernel::Gaussian, 11, None).unwrap(); + let expected = 1.06 * 100f64.powf(-0.2); + assert!((expected - 0.4219936007867).abs() < 1e-9); // pin the hand value + assert_eq!(res.bandwidth.len(), 2); + for &h in &res.bandwidth { + assert!((h - expected).abs() < 1e-12); + } +} + +/// Reads: `items[*].occ` column sums. +/// Weak near-identity on its own (documented); paired with the exact-value +/// tests above. Still kills: per-option denominators (normalizing each +/// option row separately would break cross-option coherence when combined +/// with `hand_fixture_occ_exact`), and dropped options (a missing row makes +/// sums fall short of 1). +#[test] +fn occ_rows_sum_to_one() { + let x = vec![ + vec![0.0, 2.0], + vec![1.0, 0.0], + vec![2.0, 1.0], + vec![1.0, 2.0], + vec![2.0, 2.0], + vec![0.0, 0.0], + ]; + let res = ksirt(&x, KsirtKernel::Gaussian, 7, None).unwrap(); + for item in &res.items { + assert_eq!(item.occ.len(), 3); + for s in 0..7 { + let sum: f64 = item.occ.iter().map(|row| row[s]).sum(); + assert!((sum - 1.0).abs() < 1e-12, "column {s} sums to {sum}"); + } + } +} + +/// 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 + } + 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 binary 2PL responses: P(X=1) = logistic(a*(theta - b)). +fn simulate_2pl(rng: &mut Rng, n: usize, a: &[f64], b: &[f64], skew: bool) -> Vec> { + let k = a.len(); + let mut x = Vec::with_capacity(n); + for _ in 0..n { + let z = rng.next_normal(); + // shifted lognormal with mean approx 0 for the skewed condition + let th = if skew { (0.5 * z).exp() - 1.1331 } else { z }; + let mut row = Vec::with_capacity(k); + for j in 0..k { + let p = 1.0 / (1.0 + (-a[j] * (th - b[j])).exp()); + row.push(if rng.next_f64() < p { 1.0 } else { 0.0 }); + } + x.push(row); + } + x +} + +/// Reads: `result.expected_total`. +/// Kills: theta/rank inversion (descending ranks would flip the curve) and +/// weight-matrix transposition (grid and subject axes swapped destroys the +/// monotone trend). +#[test] +fn expected_score_monotone_recovery() { + let mut rng = Rng::new(20260213); + let a = [1.2, 0.8, 1.5, 1.0, 0.9, 1.3, 1.1, 0.7, 1.4, 1.0]; + let b = [-1.0, -0.5, 0.0, 0.5, 1.0, -0.8, 0.3, 0.8, -0.2, 0.1]; + let x = simulate_2pl(&mut rng, 300, &a, &b, false); + let res = ksirt(&x, KsirtKernel::Gaussian, 51, None).unwrap(); + let q = res.expected_total.len(); + // top decile of grid clearly above bottom decile + assert!( + res.expected_total[q - 3] > res.expected_total[2] + 2.0, + "expected total not increasing: low {} high {}", + res.expected_total[2], + res.expected_total[q - 3] + ); + // and grid itself is ascending (guards a reversed-grid mutant) + assert!(res.grid[q - 1] > res.grid[0]); +} + +/// Reads: `items[0].occ` under a tiny-bandwidth uniform kernel. +/// n=2: theta = [Phi^-1(1/3), Phi^-1(2/3)], grid (q=3) endpoints coincide +/// with the two thetas. With h=0.01 only the co-located subject is in +/// support at each endpoint, and no subject is in support at the middle. +/// Kills: removal of the zero-denominator fallback (0/0 -> NaN would fail +/// the middle-point zero assertions) and support-window widening. +/// Documented limit: the boundary mutation `|u| <= 1` -> `|u| < 1` is NOT +/// killed here (u=0 at the co-located point passes both); no fixture can +/// place a subject exactly on the support edge with irrational thetas. +#[test] +fn uniform_kernel_support() { + let x = vec![vec![0.0], vec![1.0]]; + let res = ksirt(&x, KsirtKernel::Uniform, 3, Some(&[0.01])).unwrap(); + // endpoint 0: only subject 0 (score 0) in support + assert!((res.items[0].occ[0][0] - 1.0).abs() < 1e-12); + assert!(res.items[0].occ[1][0].abs() < 1e-12); + // middle: nobody in support -> all-zero fallback, finite + assert!(res.items[0].occ[0][1].abs() < 1e-12); + assert!(res.items[0].occ[1][1].abs() < 1e-12); + assert!(res.items[0].occ[0][1].is_finite()); + // endpoint 2: only subject 1 (score 1) + assert!((res.items[0].occ[1][2] - 1.0).abs() < 1e-12); +} + +/// Reads: `items[0].occ` under the quadratic kernel, h=1.5. +/// Expected values recomputed from qnorm literals with (1-u^2) truncation; +/// subject 3 falls outside the support at the low endpoint (|u| > 1), so a +/// gaussian<->quadratic swap or a dropped support check changes the value. +/// Kills: kernel dispatch mutations and support-truncation removal. +#[test] +fn quadratic_kernel_exact() { + let res = ksirt(&fixture_x(), KsirtKernel::Quadratic, 3, Some(&[1.5])).unwrap(); + let theta = fixture_theta(); + let g = QN_1_5; // low endpoint + let kw: Vec = theta + .iter() + .map(|&t| { + let u = (g - t) / 1.5; + if u.abs() <= 1.0 { + 1.0 - u * u + } else { + 0.0 + } + }) + .collect(); + assert!(kw[3] == 0.0, "fixture must exercise the out-of-support branch"); + let denom: f64 = kw.iter().sum(); + let expected = (kw[1] + kw[3]) / denom; + assert!( + (res.items[0].occ[1][0] - expected).abs() < 1e-9, + "occ[1][0] = {} expected {expected}", + res.items[0].occ[1][0] + ); +} + +/// Reads: `items[*].expected` and `result.expected_total`. +/// Kills: option-score/probability pairing mutations (expected uses the +/// sorted option scores against their own occ rows) and per-item summation +/// errors in the total. +#[test] +fn expected_score_matches_occ_combination() { + let x = vec![ + vec![0.0, 2.0], + vec![1.0, 0.0], + vec![2.0, 1.0], + vec![1.0, 2.0], + vec![2.0, 2.0], + vec![0.0, 0.0], + ]; + let res = ksirt(&x, KsirtKernel::Gaussian, 7, None).unwrap(); + for s in 0..7 { + let mut total = 0.0; + for item in &res.items { + let mut e = 0.0; + for (l, &opt) in item.options.iter().enumerate() { + e += opt * item.occ[l][s]; + } + assert!( + (item.expected[s] - e).abs() < 1e-12, + "expected[{s}] mismatch" + ); + total += e; + } + assert!((res.expected_total[s] - total).abs() < 1e-12); + } +} + +/// Reads: the `Result` error branch for each documented rejection. +/// Kills: removal of any input validation guard. +#[test] +fn input_rejection() { + let ok = fixture_x(); + assert!(ksirt(&ok[..1], KsirtKernel::Gaussian, 3, None).is_err()); // n < 2 + assert!(ksirt(&[vec![], vec![]], KsirtKernel::Gaussian, 3, None).is_err()); // k = 0 + assert!( + ksirt(&[vec![1.0], vec![1.0, 2.0]], KsirtKernel::Gaussian, 3, None).is_err(), + "ragged" + ); + assert!( + ksirt(&[vec![f64::NAN], vec![1.0]], KsirtKernel::Gaussian, 3, None).is_err(), + "NaN" + ); + assert!(ksirt(&ok, KsirtKernel::Gaussian, 1, None).is_err()); // q < 2 + assert!(ksirt(&ok, KsirtKernel::Gaussian, 3, Some(&[0.5, 0.5])).is_err()); // len + assert!(ksirt(&ok, KsirtKernel::Gaussian, 3, Some(&[0.0])).is_err()); // h <= 0 + assert!(ksirt(&ok, KsirtKernel::Gaussian, 3, Some(&[f64::NAN])).is_err()); +} + +/// Monte Carlo ICC recovery, 500 replications, normal and skewed abilities. +/// Reads: `items[*].occ` (the score-1 row) against the true 2PL ICC on the +/// central grid, averaged over items and replications. Because theta is +/// rank-based on the normal metric, the skewed condition's oracle is the +/// ICC composed with the monotone generating map t(g) = exp(g/2) - 1.1331 +/// (rank invariance: the estimate depends on abilities only through their +/// ranks, so it recovers P(X=1 | z = g) with z the normal driver). +/// Kills: gross formula errors (wrong theta metric, broken smoothing) that +/// the deterministic fixtures could miss at scale. +#[test] +#[ignore = "500-replication Monte Carlo; run explicitly"] +fn mc_2pl_recovery_500() { + let a: Vec = (0..20).map(|j| 0.7 + 0.05 * j as f64).collect(); + let b: Vec = (0..20).map(|j| -1.5 + 0.15 * j as f64).collect(); + for &skew in &[false, true] { + let mut sum_rmse = 0.0; + let mut reps = 0usize; + for rep in 0..500u64 { + let mut rng = Rng::new(7_000_003 * (rep + 1) + skew as u64); + let x = simulate_2pl(&mut rng, 500, &a, &b, skew); + let res = ksirt(&x, KsirtKernel::Gaussian, 51, None).unwrap(); + let mut se_sum = 0.0; + let mut cnt = 0usize; + for (j, item) in res.items.iter().enumerate() { + // score-1 row (options sorted ascending: [0,1]) + let row = &item.occ[item.options.len() - 1]; + for (s, &g) in res.grid.iter().enumerate() { + if g.abs() <= 1.5 { + // effective ability at normal-metric grid point g + let t = if skew { (0.5 * g).exp() - 1.1331 } else { g }; + let p = 1.0 / (1.0 + (-a[j] * (t - b[j])).exp()); + se_sum += (row[s] - p).powi(2); + cnt += 1; + } + } + } + sum_rmse += (se_sum / cnt as f64).sqrt(); + reps += 1; + } + let avg_rmse = sum_rmse / reps as f64; + assert!( + avg_rmse < 0.06, + "avg ICC RMSE {avg_rmse} (skew={skew}) exceeds 0.06" + ); + } +} diff --git a/tests/unit/subscores_tests.rs b/tests/unit/subscores_tests.rs new file mode 100644 index 000000000..0ea703ede --- /dev/null +++ b/tests/unit/subscores_tests.rs @@ -0,0 +1,278 @@ +//! Tests for Haberman (2008) subscore added-value analysis. +//! +//! Fixture literals were generated by an independent Python/NumPy script +//! transcribing the CRAN `subscore` R semantics (ddof=1 throughout), never +//! calling this crate. Each test documents the crate value its asserts read +//! and at least one mutation it kills. + +use super::*; + +const TOL: f64 = 1e-12; + +/// 10 persons x 4 binary items, subscales {0,1} and {2,3}. Chosen so all +/// alphas are in (0,1), correlations are non-degenerate, and the two +/// subscales give DIFFERENT added_value_sx outcomes (asymmetric fixture). +fn fixture() -> (Vec>, Vec) { + let x = vec![ + vec![0., 0., 0., 1.], + vec![0., 1., 0., 0.], + vec![1., 1., 0., 1.], + vec![0., 0., 1., 0.], + vec![1., 1., 1., 1.], + vec![1., 0., 1., 1.], + vec![0., 0., 0., 0.], + vec![1., 1., 1., 0.], + vec![0., 1., 1., 1.], + vec![1., 1., 0., 0.], + ]; + (x, vec![0, 0, 1, 1]) +} + +/// Reads: alpha, alpha_total. Kills: dropping the m/(m-1) factor in +/// cronbach_alpha (alpha would become [0.4348, 0.2500], alpha_x 0.3356); +/// using biased (n) variance in only one of numerator/denominator. +#[test] +fn alpha_matches_independent_reference() { + let (x, g) = fixture(); + let r = subscores(&x, &g).unwrap(); + assert!((r.alpha[0] - 0.5797101449275359).abs() < TOL); + assert!((r.alpha[1] - 0.33333333333333326).abs() < 1e-10); + assert!((r.alpha_total - 0.44742729306487705).abs() < TOL); +} + +/// Reads: prmse_x. Kills the two documented rowsum mutants: +/// including the total column in cov_k gives r_stxt from +/// cov_rowsum = [1.4333, 1.1111]; skipping the true-variance diagonal +/// replacement gives cov_rowsum = [0.8778, 0.7778]. Both move PRMSE_x far +/// outside the asserted tolerance (reference cov_rowsum = [0.5556, 0.3333], +/// r_stxt = [0.9375, 0.675]). +#[test] +fn prmse_x_matches_independent_reference() { + let (x, g) = fixture(); + let r = subscores(&x, &g).unwrap(); + assert!((r.prmse_x[0] - 0.41946308724832204).abs() < 1e-10); + assert!((r.prmse_x[1] - 0.3020134228187919).abs() < 1e-10); +} + +/// Reads: tau, beta, gamma, prmse_sx. Kills: sign flip in the tau numerator +/// (tau[0] would be negative); swapping beta and gamma (they differ by an +/// order of magnitude here); using r instead of r^2 in the tau denominator. +#[test] +fn augmented_regression_matches_independent_reference() { + let (x, g) = fixture(); + let r = subscores(&x, &g).unwrap(); + assert!((r.tau[0] - 0.1385414635898375).abs() < 1e-10); + assert!((r.tau[1] - 0.27024443691002326).abs() < 1e-10); + assert!((r.beta[0] - 0.49752475247524747).abs() < 1e-10); + assert!((r.beta[1] - 0.21782178217821788).abs() < 1e-10); + assert!((r.gamma[0] - 0.07178217821782161).abs() < 1e-10); + assert!((r.gamma[1] - 0.09900990099009896).abs() < 1e-10); + assert!((r.prmse_sx[0] - 0.5872524752475244).abs() < 1e-10); + assert!((r.prmse_sx[1] - 0.3663366336633662).abs() < 1e-10); +} + +/// Reads: corr (K+1 x K+1, total last), prmse_s. Kills: variance/covariance +/// index transposition; normalizing by variance instead of SD. +#[test] +fn correlations_match_independent_reference() { + let (x, g) = fixture(); + let r = subscores(&x, &g).unwrap(); + assert!((r.corr[0][2] - 0.7791290756515445).abs() < 1e-10); + assert!((r.corr[1][2] - 0.7403367031320777).abs() < 1e-10); + assert!((r.corr[2][0] - 0.7791290756515445).abs() < 1e-10); + assert!((r.prmse_s[0] - r.alpha[0]).abs() < TOL); + assert!((r.prmse_s[1] - r.alpha[1]).abs() < TOL); + for k in 0..3 { + assert!((r.corr[k][k] - 1.0).abs() < TOL); + } +} + +/// Reads: subscore_s, subscore_x, subscore_sx for person 0 (s = [0, 1], +/// x = 1). Kills: using observed instead of true SD in the s_hat_x slope +/// (0.7309 vs 0.6152 for subscale 0); centering on x instead of s in +/// s_hat_s; swapping mean_s and mean_x. +#[test] +fn person_estimates_match_independent_reference() { + let (x, g) = fixture(); + let r = subscores(&x, &g).unwrap(); + assert!((r.subscore_s[0][0] - 0.4623188405797105).abs() < 1e-10); + assert!((r.subscore_s[0][1] - 1.0).abs() < 1e-10); + assert!((r.subscore_x[0][0] - 0.7308724832214768).abs() < 1e-10); + assert!((r.subscore_x[0][1] - 0.778523489932886).abs() < 1e-10); + assert!((r.subscore_sx[0][0] - 0.47376237623762407).abs() < 1e-10); + assert!((r.subscore_sx[0][1] - 0.8910891089108911).abs() < 1e-10); + assert!((r.observed[0][0] - 0.0).abs() < TOL); + assert!((r.observed[0][1] - 1.0).abs() < TOL); + assert!((r.total[0] - 1.0).abs() < TOL); +} + +/// Reads: added_value_s, added_value_sx. The fixture is deliberately +/// asymmetric across subscales in added_value_sx ([false, true]), so any +/// mutant collapsing the decision to a constant, comparing against the wrong +/// PRMSE, or dropping the +0.01 margin (subscale 0 has +/// PRMSE_sx - PRMSE_s = 0.00754 < 0.01, which flips to true without the +/// margin) fails here. +#[test] +fn added_value_decisions_match_independent_reference() { + let (x, g) = fixture(); + let r = subscores(&x, &g).unwrap(); + assert_eq!(r.added_value_s, vec![true, true]); + assert_eq!(r.added_value_sx, vec![false, true]); +} + +/// Reads: Err from negatively-keyed within-subscale items (alpha < 0). +/// Kills: removing the alpha-range guard (analysis would emit NaN PRMSEs). +#[test] +fn rejects_negative_alpha() { + // Item 1 is the reverse of item 0 within subscale 0 -> alpha_0 < 0. + let x: Vec> = (0..10) + .map(|p| { + let a = f64::from((p % 2 == 0) as u8); + let c = f64::from((p % 3 == 0) as u8); + let d = f64::from((p % 4 == 0) as u8); + vec![a, 1.0 - a, c, d] + }) + .collect(); + let err = subscores(&x, &[0, 0, 1, 1]).unwrap_err(); + assert!(err.contains("outside (0, 1]"), "{err}"); +} + +/// Reads: Err from a subscore collinear with the total. Duplicated subscales +/// make corr(s_k, x) = 1 exactly; the review mandated this be a REJECTION +/// (tau divides by 1 - r^2 = 0), not a parallel-forms anchor. Kills: +/// removing the collinearity guard (tau/PRMSE_sx would be inf/NaN). +#[test] +fn rejects_collinear_subscore() { + let (x, _) = fixture(); + // Both subscales are the SAME two items duplicated: x = 2 s_1. + let dup: Vec> = x + .iter() + .map(|row| vec![row[0], row[1], row[0], row[1]]) + .collect(); + let err = subscores(&dup, &[0, 0, 1, 1]).unwrap_err(); + assert!(err.contains("collinear"), "{err}"); +} + +/// Reads: Err on structural validation. Kills: skipping partition checks +/// (out-of-bounds panic or silent nonsense instead of Err). +#[test] +fn rejects_bad_inputs() { + let (x, _) = fixture(); + assert!(subscores(&x[..2], &[0, 0, 1, 1]).is_err()); // n < 3 + assert!(subscores(&x, &[0, 0, 1]).is_err()); // groups length mismatch + assert!(subscores(&x, &[0, 0, 0, 0]).is_err()); // K = 1 + assert!(subscores(&x, &[0, 1, 1, 1]).is_err()); // singleton subscale + assert!(subscores(&x, &[0, 0, 1_000_000_000, 1_000_000_000]).is_err()); // sparse index DoS + let mut nanx = x.clone(); + nanx[0][0] = f64::NAN; + assert!(subscores(&nanx, &[0, 0, 1, 1]).is_err()); // incomplete data + let zero = vec![vec![0., 0., 1., 0.]; 5]; + assert!(subscores(&zero, &[0, 0, 1, 1]).is_err()); // zero variance +} + +/// Reads: prmse_sx vs prmse_s/prmse_x across many random guard-passing +/// datasets. The dominance PRMSE_sx >= max(PRMSE_s, PRMSE_x) is CONDITIONAL +/// on the guards (spec-review finding); asserted only on accepted samples. +/// Note: PRMSE_sx >= PRMSE_s is the identity alpha + tau^2 (1 - r^2) >= +/// alpha, which no mutant of the comparison inputs can violate on its own; +/// the discriminating half is PRMSE_sx >= PRMSE_x, which fails under the +/// rowsum/diagonal mutants documented in `prmse_x_matches_independent_reference`. +#[test] +fn dominance_holds_on_guard_passing_data() { + // Deterministic LCG so the test is reproducible without rand dep. + let mut state = 0x2545F4914F6CDD1Du64; + let mut next = move || { + state = state.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407); + ((state >> 33) as f64) / ((1u64 << 31) as f64) + }; + let mut accepted = 0; + for _ in 0..200 { + // Common-factor items so alphas are usually positive. + let n = 40; + let x: Vec> = (0..n) + .map(|_| { + let f = next(); + (0..6) + .map(|_| f64::from((0.6 * f + 0.4 * next() > 0.5) as u8)) + .collect() + }) + .collect(); + if let Ok(r) = subscores(&x, &[0, 0, 0, 1, 1, 1]) { + accepted += 1; + for k in 0..2 { + assert!( + r.prmse_sx[k] + 1e-9 >= r.prmse_s[k].max(r.prmse_x[k]), + "dominance violated: sx={} s={} x={}", + r.prmse_sx[k], + r.prmse_s[k], + r.prmse_x[k] + ); + } + } + } + assert!(accepted >= 50, "too few guard-passing samples: {accepted}"); +} + +/// Reads: disattenuated_corr. Kills: dividing by alpha instead of +/// sqrt(alpha_a alpha_b); forgetting the NaN diagonal. +#[test] +fn disattenuated_correlation() { + let (x, g) = fixture(); + let r = subscores(&x, &g).unwrap(); + let expect = 0.35355339059327384; + assert!((r.disattenuated_corr[0][1] - expect).abs() < 1e-10); + assert!((r.corr[0][1] - 0.15541746804005227).abs() < 1e-10); + assert!(r.disattenuated_corr[0][0].is_nan()); +} + +/// 500-replication Monte Carlo: with a strong common factor plus weak +/// specific factors, the augmented estimator's MSE against the true subscore +/// should not exceed the observed-subscore estimator's MSE on average. +/// Reads: subscore_sx, subscore_s. Kills gross beta/gamma miscalibration +/// that the single-fixture literals could compensate for by luck. +#[test] +#[ignore = "500-rep Monte Carlo; run with -- --ignored"] +fn monte_carlo_augmented_mse() { + let mut state = 0x9E3779B97F4A7C15u64; + let mut next = move || { + state = state.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407); + ((state >> 33) as f64) / ((1u64 << 31) as f64) + }; + let mut mse_s = 0.0f64; + let mut mse_sx = 0.0f64; + let mut reps = 0; + for _ in 0..500 { + let n = 200; + let mut x = Vec::with_capacity(n); + let mut true_s0 = Vec::with_capacity(n); + for _ in 0..n { + let f = next() * 2.0 - 1.0; // common factor + let u0 = next() * 0.6 - 0.3; // weak specific factors + let u1 = next() * 0.6 - 0.3; + let mut row = Vec::with_capacity(8); + let mut t0 = 0.0; + for j in 0..8 { + let load = if j < 4 { f + u0 } else { f + u1 }; + let p = 1.0 / (1.0 + (-2.0 * load).exp()); + if j < 4 { + t0 += p; + } + row.push(f64::from((next() < p) as u8)); + } + true_s0.push(t0); + x.push(row); + } + if let Ok(r) = subscores(&x, &[0, 0, 0, 0, 1, 1, 1, 1]) { + reps += 1; + for p in 0..n { + mse_s += (r.subscore_s[p][0] - true_s0[p]).powi(2); + mse_sx += (r.subscore_sx[p][0] - true_s0[p]).powi(2); + } + } + } + assert!(reps >= 400, "too few valid replications: {reps}"); + assert!( + mse_sx <= mse_s * 1.001, + "augmented MSE {mse_sx} exceeds observed-subscore MSE {mse_s}" + ); +}