diff --git a/CHANGELOG.md b/CHANGELOG.md index 494fea9b1..bf870b754 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -114,6 +114,25 @@ ### Added +- **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 + (Andrich, 1978) with a rater facet — to a `persons x items x raters` array with + NaN-missing sparse judging plans. For LLM-as-a-Judge calibration this puts each + judge's severity `c_j` on a common logit scale adjusted for item difficulty and + respondent ability. Estimation is marginal-ML EM on a Gauss-Hermite grid + (Bock & Aitkin, 1981), NOT Linacre's JMLE, and the docs say so: estimates match + the Facets program only up to the JMLE-vs-MMLE difference. Identification: + `theta ~ N(0,1)`, severities and thresholds centered to sum 0 + (`n_parameters = I + (J-1) + (K-2)`). Reports Linacre's connectedness + diagnostic via union-find over the person-mediated item∪rater co-observation + graph; `connected=False` means cross-component severity comparisons rest + solely on the shared trait prior, not the rating design. Rust-only numerics; + the Python wrapper validates and marshals. Tests include FD gradient anchors, + the J=1 RSM-reduction identity, asymmetric-severity recovery, sparse and + disconnected designs, and an `#[ignore]` 500-replicate Monte Carlo + (normal + skew-normal traits) bounding severity bias and RMSE; a gradient + sign-flip mutant was verified to fail 4 tests. - **Warm's weighted-likelihood ability estimation for POLYTOMOUS items** (`fast_mlsirm.score_wle_poly`; new `score_wle_poly` in `mlsirm_core::scoring`; Warm, 1989). The library already had the full polytomous model family and polytomous EAP scoring, but its only bias-reduced ML ability estimator was diff --git a/crates/fast-mlsirm-py/src/lib.rs b/crates/fast-mlsirm-py/src/lib.rs index 13681fe7a..437e33ea8 100644 --- a/crates/fast-mlsirm-py/src/lib.rs +++ b/crates/fast-mlsirm-py/src/lib.rs @@ -59,6 +59,7 @@ use mlsirm_core::poly_marginal::fit_poly_lsirm as core_fit_poly_lsirm; 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::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, @@ -1407,6 +1408,68 @@ fn fit_rsm( Ok(out.into()) } +/// Many-Facet Rasch Model fit (Linacre, 1989; `mlsirm_core::facets::fit_facets`). +/// `y`/`observed` are row-major `n_persons * n_items * n_raters` (rater fastest) +/// with categories `0..n_cat-1`. Adjacent-category log-odds: +/// `ln[P(k)/P(k-1)] = theta - item_difficulty_i - rater_severity_j - threshold_k`, +/// `theta ~ N(0,1)`; severities and thresholds are centered to sum 0. Returns a +/// dict with `item_difficulty` (`n_items`), `rater_severity` (`n_raters`), +/// `thresholds` (`n_cat-1`), `theta` (per-person EAP), `loglik_trace`, `n_iter`, +/// `converged`, `connected` (design-linking flag), `n_parameters`. +#[pyfunction] +#[allow(clippy::too_many_arguments)] +#[pyo3(signature = (y, observed, n_persons, n_items, n_raters, n_cat, q_theta = 41, max_iter = 500, tol = 1e-6))] +fn fit_facets( + py: Python<'_>, + y: PyReadonlyArray1<'_, i64>, + observed: PyReadonlyArray1<'_, bool>, + n_persons: usize, + n_items: usize, + n_raters: usize, + n_cat: usize, + q_theta: usize, + max_iter: usize, + tol: f64, +) -> PyResult> { + let yy: Vec = y + .as_slice()? + .iter() + .map(|&v| { + if v >= 0 { + Ok(v as usize) + } else { + Err(PyValueError::new_err( + "y must be non-negative category indices", + )) + } + }) + .collect::>()?; + let obs = observed.as_slice()?; + let res = core_fit_facets( + &yy, + Some(obs), + n_persons, + n_items, + n_raters, + n_cat, + q_theta, + max_iter, + tol, + ) + .map_err(PyValueError::new_err)?; + let out = pyo3::types::PyDict::new(py); + out.set_item("item_difficulty", res.item_difficulty)?; + out.set_item("rater_severity", res.rater_severity)?; + out.set_item("thresholds", res.thresholds)?; + out.set_item("theta", res.theta)?; + out.set_item("loglik_trace", res.loglik_trace)?; + out.set_item("n_iter", res.n_iter)?; + out.set_item("converged", res.converged)?; + out.set_item("connected", res.connected)?; + out.set_item("n_parameters", res.n_parameters)?; + 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 /// "2pl". `n_classes` latent classes each get their own item parameters. Returns a dict @@ -5002,6 +5065,7 @@ fn fast_mlsirm_core(m: &Bound<'_, PyModule>) -> PyResult<()> { m.add_function(wrap_pyfunction!(fit_gpcm, m)?)?; 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!(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/facets.rs b/crates/mlsirm-core/src/facets.rs new file mode 100644 index 000000000..b3876efa0 --- /dev/null +++ b/crates/mlsirm-core/src/facets.rs @@ -0,0 +1,627 @@ +//! Many-Facet Rasch Model (MFRM; Linacre, 1989) by marginal-ML EM. +//! +//! The MFRM extends the rating scale model with a rater facet: each rating of +//! person `p` on item `i` by rater `j` follows the adjacent-category log-odds +//! +//! ```text +//! ln[ P(Y_pij = k | theta) / P(Y_pij = k-1 | theta) ] +//! = theta_p - d_i - c_j - f_k, k = 1..K-1, +//! ``` +//! +//! with item difficulty `d_i`, rater severity `c_j`, and category thresholds +//! `f_k` shared across items and raters (the rating-scale form of Linacre's +//! model). The cumulative predictor is `psi_k = k*theta - k*(d_i + c_j) - T_k` +//! with `T_k = sum_{m<=k} f_m`, `psi_0 = 0`, `P(Y=k) = softmax_k(psi)` — exactly +//! the RSM cell ([`crate::rsm::rsm_logprobs`]) with location `d_i + c_j`, which +//! this module reuses. At `n_raters = 1` (severity centered to 0) the model +//! reduces to the RSM. +//! +//! Verified formulation: the adjacent-category identity +//! `psi_k - psi_{k-1} = theta - d_i - c_j - f_k` was re-derived here and +//! adversarially checked; it matches the published Linacre (1989) rating-scale +//! MFRM form as documented by Eckes (2015). We did **not** reproduce Linacre's +//! JMLE estimator: Facets-style JMLE is replaced by marginal ML (Bock & Aitkin, +//! 1981) with `theta ~ N(0,1)` on a Gauss-Hermite grid, matching this crate's +//! estimation contract (see `rsm.rs`, `mixed.rs`). Parameter estimates are +//! therefore comparable to Facets output only up to the JMLE-vs-MMLE difference. +//! +//! Identification: the probabilities are invariant under +//! `f_m -> f_m - c, d_i -> d_i + c` and under `c_j -> c_j - c, d_i -> d_i + c`; +//! the trait scale is fixed by `theta ~ N(0,1)`. Both shift redundancies are +//! removed by centering `sum_m f_m = 0` and `sum_j c_j = 0`, leaving +//! `n_items + (n_raters - 1) + (n_cat - 2)` free parameters. +//! +//! Connectedness: Linacre's connectedness requirement concerns *design* +//! linking. We report a `connected` flag from a union-find over facet elements +//! (items and raters), joining every element that appears in the same person's +//! observed cells. When `connected == false`, severity/difficulty comparisons +//! across components are anchored only by the shared `theta ~ N(0,1)` +//! assumption (model-prior linking), not by the rating design itself. +//! +//! # References (APA 7th ed.) +//! Bock, R. D., & Aitkin, M. (1981). Marginal maximum likelihood estimation of +//! item parameters: Application of an EM algorithm. *Psychometrika, 46*(4), +//! 443-459. https://doi.org/10.1007/BF02293801 +//! Eckes, T. (2015). *Introduction to many-facet Rasch measurement* (2nd ed.). +//! Peter Lang. https://doi.org/10.3726/978-3-653-04844-5 +//! Linacre, J. M. (1989). *Many-facet Rasch measurement*. MESA Press. +//! Andrich, D. (1978). A rating formulation for ordered response categories. +//! *Psychometrika, 43*(4), 561-573. https://doi.org/10.1007/BF02293814 + +use crate::poly::solve_small; +use crate::rsm::rsm_logprobs; + +const FACETS_MAX_CAT: usize = 64; +const FACETS_MAX_ITER: usize = 100_000; +const FACETS_MAX_CELLS: usize = 60_000_000; + +/// Fitted many-facet Rasch model (Linacre, 1989). `item_difficulty` is `d_i`; +/// `rater_severity` the centered `c_j` (`sum = 0`, higher = harsher); +/// `thresholds` the `K-1` common category thresholds (centered, `sum = 0`); +/// `theta` the per-person EAP trait. `connected` is the design-linking flag +/// (see module docs). +#[derive(Clone, Debug)] +pub struct FacetsResult { + pub item_difficulty: Vec, + pub rater_severity: Vec, + pub thresholds: Vec, + pub theta: Vec, + pub loglik_trace: Vec, + pub n_iter: usize, + pub converged: bool, + pub connected: bool, + /// `n_items + (n_raters - 1) + (n_cat - 2)`. + pub n_parameters: usize, +} + +/// Fit the many-facet Rasch model (Linacre, 1989) by marginal-ML EM. `y` is +/// `n_persons * n_items * n_raters` row-major (rater fastest) categories +/// `0..n_cat-1`; `observed` marks scored cells (sparse judging plans allowed; +/// `None` = fully crossed). Ability `theta ~ N(0,1)` on the `q_theta`-node +/// Gauss-Hermite grid. +#[allow(clippy::too_many_arguments)] +pub fn fit_facets( + y: &[usize], + observed: Option<&[bool]>, + n_persons: usize, + n_items: usize, + n_raters: usize, + n_cat: usize, + q_theta: usize, + max_iter: usize, + tol: f64, +) -> Result { + if !(2..=FACETS_MAX_CAT).contains(&n_cat) { + return Err(format!("n_cat must be in 2..={FACETS_MAX_CAT}")); + } + if n_persons < 1 || n_items < 1 || n_raters < 1 { + return Err("n_persons, n_items and n_raters must be >= 1".into()); + } + if !(1..=FACETS_MAX_ITER).contains(&max_iter) { + return Err(format!("max_iter must be in 1..={FACETS_MAX_ITER}")); + } + if !tol.is_finite() || tol <= 0.0 { + return Err("tol must be finite and > 0".into()); + } + let n_pairs = + crate::checked_mul_usize(n_items, n_raters, "n_items * n_raters overflows usize")?; + let n_cells = + crate::checked_mul_usize(n_persons, n_pairs, "n_persons * n_items * n_raters overflows")?; + if y.len() != n_cells { + return Err("y must have length n_persons * n_items * n_raters".into()); + } + if let Some(o) = observed { + if o.len() != n_cells { + return Err("observed must have length n_persons * n_items * n_raters".into()); + } + } + for (idx, &cat) in y.iter().enumerate() { + if observed.map_or(true, |o| o[idx]) && cat >= n_cat { + return Err("response category out of range 0..n_cat-1".into()); + } + } + let is_obs = |p: usize, pair: usize| observed.map_or(true, |o| o[p * n_pairs + pair]); + for i in 0..n_items { + if !(0..n_persons) + .any(|p| (0..n_raters).any(|j| is_obs(p, i * n_raters + j))) + { + return Err(format!("item {i} has no observed responses")); + } + } + for j in 0..n_raters { + if !(0..n_persons).any(|p| (0..n_items).any(|i| is_obs(p, i * n_raters + j))) { + return Err(format!("rater {j} has no observed responses")); + } + } + let (nodes, weights) = crate::quadrature::gh_rule(q_theta) + .ok_or_else(|| format!("unsupported q_theta {q_theta}"))?; + let qn = nodes.len(); + let count_cells = n_pairs + .checked_mul(qn) + .and_then(|c| c.checked_mul(n_cat)) + .ok_or_else(|| "pair * node * category table size overflows usize".to_string())?; + if count_cells > FACETS_MAX_CELLS { + return Err(format!( + "count table {count_cells} cells exceeds the cap {FACETS_MAX_CELLS}" + )); + } + let log_w: Vec = weights.iter().map(|w| w.ln()).collect(); + let kb = n_cat - 1; + + let connected = design_connected(n_persons, n_items, n_raters, &is_obs); + + // Init: item difficulty from the item mean category (as in rsm.rs), rater + // severity and thresholds at 0. + let mut d = vec![0.0f64; n_items]; + let mut c = vec![0.0f64; n_raters]; + let mut f = vec![0.0f64; kb]; + for i in 0..n_items { + let (mut s, mut cnt) = (0.0f64, 0.0f64); + for p in 0..n_persons { + for j in 0..n_raters { + if is_obs(p, i * n_raters + j) { + s += y[p * n_pairs + i * n_raters + j] as f64; + cnt += 1.0; + } + } + } + if cnt > 0.0 { + let mean = s / cnt / kb as f64; + d[i] = ((1.0 - mean).clamp(0.02, 0.98) / mean.clamp(0.02, 0.98)).ln(); + } + } + + let mut it = 0usize; + let mut converged = false; + let mut loglik_trace: Vec = Vec::new(); + + while it < max_iter { + // Per-pair cell log-probs at each node (RSM cell, location d_i + c_j). + let item_lp = pair_logprobs(&d, &c, &f, nodes, n_items, n_raters, n_cat); + // E-step: posteriors -> expected counts r[pair][node][k]. + let mut r = vec![vec![0.0f64; qn * n_cat]; n_pairs]; + let mut ll = 0.0f64; + let mut log_node = vec![0.0f64; qn]; + for p in 0..n_persons { + log_node[..qn].copy_from_slice(&log_w[..qn]); + for pair in 0..n_pairs { + if !is_obs(p, pair) { + continue; + } + let yc = y[p * n_pairs + pair]; + for nd in 0..qn { + log_node[nd] += item_lp[pair][nd * n_cat + yc]; + } + } + let mx = log_node.iter().cloned().fold(f64::NEG_INFINITY, f64::max); + let mut denom = 0.0f64; + for nd in 0..qn { + denom += (log_node[nd] - mx).exp(); + } + ll += mx + denom.ln(); + for pair in 0..n_pairs { + if !is_obs(p, pair) { + continue; + } + let yc = y[p * n_pairs + pair]; + for nd in 0..qn { + r[pair][nd * n_cat + yc] += (log_node[nd] - mx).exp() / denom; + } + } + } + + loglik_trace.push(ll); + it += 1; + if loglik_trace.len() > 1 { + let nn = loglik_trace.len(); + if (loglik_trace[nn - 1] - loglik_trace[nn - 2]).abs() + < tol * (1.0 + loglik_trace[nn - 2].abs()) + { + converged = true; + break; + } + } + + // CM-1: per-item Newton on d_i (c, f fixed), aggregated over raters. + // g = -sum_{j,nd,k} k*(r - n*P); h = -sum n*Var(score) < 0 (score = k). + for i in 0..n_items { + for _ in 0..25 { + let (mut g, mut h) = (0.0f64, 0.0f64); + for j in 0..n_raters { + let pair = i * n_raters + j; + location_score_terms( + d[i] + c[j], + &f, + &r[pair], + nodes, + n_cat, + &mut g, + &mut h, + ); + } + if h >= 0.0 { + break; + } + let step = g / h; + let cur = item_ell_d(i, &d, &c, &f, &r, nodes, n_raters, n_cat); + let mut al = 1.0f64; + let mut accepted = false; + for _ in 0..24 { + let cand = d[i] - al * step; + let mut dc = d.clone(); + dc[i] = cand; + if item_ell_d(i, &dc, &c, &f, &r, nodes, n_raters, n_cat) >= cur - 1e-12 { + d[i] = cand; + accepted = true; + break; + } + al *= 0.5; + } + if !accepted || (al * step).abs() < 1e-9 { + break; + } + } + } + + // CM-2: per-rater Newton on c_j (d, f fixed), aggregated over items — + // same algebra as CM-1 by the d<->c symmetry of the location d_i + c_j. + for j in 0..n_raters { + for _ in 0..25 { + let (mut g, mut h) = (0.0f64, 0.0f64); + for i in 0..n_items { + let pair = i * n_raters + j; + location_score_terms( + d[i] + c[j], + &f, + &r[pair], + nodes, + n_cat, + &mut g, + &mut h, + ); + } + if h >= 0.0 { + break; + } + let step = g / h; + let cur = rater_ell_c(j, &d, &c, &f, &r, nodes, n_items, n_raters, n_cat); + let mut al = 1.0f64; + let mut accepted = false; + for _ in 0..24 { + let cand = c[j] - al * step; + let mut cc = c.clone(); + cc[j] = cand; + if rater_ell_c(j, &d, &cc, &f, &r, nodes, n_items, n_raters, n_cat) + >= cur - 1e-12 + { + c[j] = cand; + accepted = true; + break; + } + al *= 0.5; + } + if !accepted || (al * step).abs() < 1e-9 { + break; + } + } + } + + // CM-3: joint Newton on the common thresholds f (d, c fixed), + // aggregated over all (item, rater) pairs; FD Hessian of the gradient. + for _ in 0..25 { + let g = f_gradient(&f, &d, &c, &r, nodes, n_items, n_raters, n_cat); + let mut hess = vec![vec![0.0f64; kb]; kb]; + let eps = 1e-5; + for jj in 0..kb { + let mut fp = f.clone(); + fp[jj] += eps; + let gj = f_gradient(&fp, &d, &c, &r, nodes, n_items, n_raters, n_cat); + for a in 0..kb { + hess[a][jj] = (gj[a] - g[a]) / eps; + } + } + for a in 0..kb { + for b in 0..kb { + hess[a][b] = 0.5 * (hess[a][b] + hess[b][a]); + } + hess[a][a] -= 1e-8; + } + let step = solve_small(hess, g.clone()); + let cur = total_ell(&d, &c, &f, &r, nodes, n_items, n_raters, n_cat); + let mut al = 1.0f64; + let mut accepted = false; + let mut max_step = 0.0f64; + for _ in 0..24 { + let cand: Vec = (0..kb).map(|m| f[m] - al * step[m]).collect(); + if total_ell(&d, &c, &cand, &r, nodes, n_items, n_raters, n_cat) >= cur - 1e-12 { + max_step = (0..kb).map(|m| (al * step[m]).abs()).fold(0.0, f64::max); + f = cand; + accepted = true; + break; + } + al *= 0.5; + } + if !accepted || max_step < 1e-9 { + break; + } + } + + // Re-center: f_m -> f_m - cf shifts T_k by -k*cf, compensated by + // d_i -> d_i + cf (psi_k regains -k*cf through the k*(d+c) term). + let cf = f.iter().sum::() / kb as f64; + for fm in f.iter_mut() { + *fm -= cf; + } + for di in d.iter_mut() { + *di += cf; + } + // c_j -> c_j - cc, d_i -> d_i + cc leaves every location d_i + c_j fixed. + let cc = c.iter().sum::() / n_raters as f64; + for cj in c.iter_mut() { + *cj -= cc; + } + for di in d.iter_mut() { + *di += cc; + } + } + + // Final person EAP pass at the returned parameters. + let item_lp = pair_logprobs(&d, &c, &f, nodes, n_items, n_raters, n_cat); + let mut theta = vec![0.0f64; n_persons]; + let mut final_ll = 0.0f64; + let mut log_node = vec![0.0f64; qn]; + for p in 0..n_persons { + log_node[..qn].copy_from_slice(&log_w[..qn]); + for pair in 0..n_pairs { + if !is_obs(p, pair) { + continue; + } + let yc = y[p * n_pairs + pair]; + for nd in 0..qn { + log_node[nd] += item_lp[pair][nd * n_cat + yc]; + } + } + let mx = log_node.iter().cloned().fold(f64::NEG_INFINITY, f64::max); + let mut denom = 0.0f64; + for nd in 0..qn { + denom += (log_node[nd] - mx).exp(); + } + final_ll += mx + denom.ln(); + let mut m = 0.0f64; + for (nd, &node) in nodes.iter().enumerate() { + m += (log_node[nd] - mx).exp() / denom * node; + } + theta[p] = m; + } + if !converged { + loglik_trace.push(final_ll); + } + + Ok(FacetsResult { + item_difficulty: d, + rater_severity: c, + thresholds: f, + theta, + loglik_trace, + n_iter: it, + converged, + connected, + n_parameters: n_items + (n_raters - 1) + (n_cat - 2), + }) +} + +/// Per-pair RSM cell log-prob tables: `out[i*n_raters + j][nd*n_cat + k]`. +fn pair_logprobs( + d: &[f64], + c: &[f64], + f: &[f64], + nodes: &[f64], + n_items: usize, + n_raters: usize, + n_cat: usize, +) -> Vec> { + let qn = nodes.len(); + let mut out = vec![vec![0.0f64; qn * n_cat]; n_items * n_raters]; + for i in 0..n_items { + for j in 0..n_raters { + let pair = i * n_raters + j; + for (nd, &theta) in nodes.iter().enumerate() { + let lp = rsm_logprobs(theta, d[i] + c[j], f); + out[pair][nd * n_cat..(nd + 1) * n_cat].copy_from_slice(&lp); + } + } + } + out +} + +/// Accumulate the location gradient/Hessian terms of one (pair) count block: +/// `g += -sum_{nd,k} k*(r - n*P)`, `h += -sum_nd n*Var(score)`. Shared by the +/// `d_i` and `c_j` Newton steps (`d ln P_k / d location = -k + E[score]`). +fn location_score_terms( + location: f64, + f: &[f64], + r_pair: &[f64], + nodes: &[f64], + n_cat: usize, + g: &mut f64, + h: &mut f64, +) { + for (nd, &theta) in nodes.iter().enumerate() { + let lp = rsm_logprobs(theta, location, f); + let mut n = 0.0f64; + for k in 0..n_cat { + n += r_pair[nd * n_cat + k]; + } + if n <= 0.0 { + continue; + } + let (mut e1, mut e2) = (0.0f64, 0.0f64); + for k in 0..n_cat { + let pk = lp[k].exp(); + let kf = k as f64; + e1 += kf * pk; + e2 += kf * kf * pk; + *g += -kf * (r_pair[nd * n_cat + k] - n * pk); + } + *h += -n * (e2 - e1 * e1); + } +} + +/// Expected complete-data log-lik of the cells involving item `i` (its row of +/// rater pairs) — the objective ascended by the `d_i` line search. +#[allow(clippy::too_many_arguments)] +fn item_ell_d( + i: usize, + d: &[f64], + c: &[f64], + f: &[f64], + r: &[Vec], + nodes: &[f64], + n_raters: usize, + n_cat: usize, +) -> f64 { + (0..n_raters) + .map(|j| pair_ell(d[i] + c[j], f, &r[i * n_raters + j], nodes, n_cat)) + .sum() +} + +/// Expected complete-data log-lik of the cells involving rater `j` — the +/// objective ascended by the `c_j` line search. +#[allow(clippy::too_many_arguments)] +fn rater_ell_c( + j: usize, + d: &[f64], + c: &[f64], + f: &[f64], + r: &[Vec], + nodes: &[f64], + n_items: usize, + n_raters: usize, + n_cat: usize, +) -> f64 { + (0..n_items) + .map(|i| pair_ell(d[i] + c[j], f, &r[i * n_raters + j], nodes, n_cat)) + .sum() +} + +/// `sum_nd sum_k r[nd][k] * log P(k | theta_nd; location, f)` for one pair. +fn pair_ell(location: f64, f: &[f64], r_pair: &[f64], nodes: &[f64], n_cat: usize) -> f64 { + let mut acc = 0.0f64; + for (nd, &theta) in nodes.iter().enumerate() { + let lp = rsm_logprobs(theta, location, f); + for k in 0..n_cat { + let rc = r_pair[nd * n_cat + k]; + if rc != 0.0 { + acc += rc * lp[k]; + } + } + } + acc +} + +/// Total expected complete-data log-lik over all pairs (for the shared-`f` +/// line search). +#[allow(clippy::too_many_arguments)] +fn total_ell( + d: &[f64], + c: &[f64], + f: &[f64], + r: &[Vec], + nodes: &[f64], + n_items: usize, + n_raters: usize, + n_cat: usize, +) -> f64 { + let mut acc = 0.0f64; + for i in 0..n_items { + for j in 0..n_raters { + acc += pair_ell(d[i] + c[j], f, &r[i * n_raters + j], nodes, n_cat); + } + } + acc +} + +/// Gradient of the expected complete-data objective w.r.t. the common +/// thresholds: `g_m = -sum_{i,j,nd} sum_{k>=m} (r - n*P)` (0-indexed `m` for +/// `f_{m+1}`); suffix-residual form as in `rsm::tau_gradient`. +#[allow(clippy::too_many_arguments)] +fn f_gradient( + f: &[f64], + d: &[f64], + c: &[f64], + r: &[Vec], + nodes: &[f64], + n_items: usize, + n_raters: usize, + n_cat: usize, +) -> Vec { + let kb = f.len(); + let mut g = vec![0.0f64; kb]; + for i in 0..n_items { + for j in 0..n_raters { + let pair = i * n_raters + j; + for (nd, &theta) in nodes.iter().enumerate() { + let lp = rsm_logprobs(theta, d[i] + c[j], f); + let mut n = 0.0f64; + for k in 0..n_cat { + n += r[pair][nd * n_cat + k]; + } + if n <= 0.0 { + continue; + } + let mut suffix = 0.0f64; + for k in (1..n_cat).rev() { + suffix += r[pair][nd * n_cat + k] - n * lp[k].exp(); + g[k - 1] += -suffix; + } + } + } + } + g +} + +/// Design-linking flag: union-find over facet elements (`n_items` item nodes, +/// then `n_raters` rater nodes), joining every element observed for the same +/// person. `true` iff all items and raters form one component. Persons anchor +/// components to the trait scale only through `theta ~ N(0,1)` (module docs). +fn design_connected( + n_persons: usize, + n_items: usize, + n_raters: usize, + is_obs: &dyn Fn(usize, usize) -> bool, +) -> bool { + let n = n_items + n_raters; + let mut parent: Vec = (0..n).collect(); + fn find(parent: &mut [usize], mut x: usize) -> usize { + while parent[x] != x { + parent[x] = parent[parent[x]]; + x = parent[x]; + } + x + } + for p in 0..n_persons { + let mut first: Option = None; + for i in 0..n_items { + for j in 0..n_raters { + if !is_obs(p, i * n_raters + j) { + continue; + } + for node in [i, n_items + j] { + match first { + None => first = Some(node), + Some(anchor) => { + let (ra, rb) = (find(&mut parent, anchor), find(&mut parent, node)); + parent[rb] = ra; + } + } + } + } + } + } + let root = find(&mut parent, 0); + (1..n).all(|x| find(&mut parent, x) == root) +} + +#[cfg(test)] +#[path = "../../../tests/unit/facets_tests.rs"] +mod tests; diff --git a/crates/mlsirm-core/src/lib.rs b/crates/mlsirm-core/src/lib.rs index 36fd85aa2..e55b47ff6 100644 --- a/crates/mlsirm-core/src/lib.rs +++ b/crates/mlsirm-core/src/lib.rs @@ -3,6 +3,7 @@ pub mod cdm; pub mod crm; pub mod dif; pub mod equating; +pub mod facets; pub mod fitstats; pub mod gpcm; pub mod grm; diff --git a/python/fast_mlsirm/__init__.py b/python/fast_mlsirm/__init__.py index 69cbba3ec..d3121b1e9 100644 --- a/python/fast_mlsirm/__init__.py +++ b/python/fast_mlsirm/__init__.py @@ -32,6 +32,7 @@ from .nominal import fit_nominal as fit_nominal, NominalResponseFit as NominalResponseFit 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 .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 @@ -143,6 +144,8 @@ "GpcmFit", "fit_rsm", "RsmFit", + "fit_facets", + "FacetsFit", "fit_mixed_items", "MixedFormatFit", "MixedItemParameters", diff --git a/python/fast_mlsirm/facets.py b/python/fast_mlsirm/facets.py new file mode 100644 index 000000000..57266c588 --- /dev/null +++ b/python/fast_mlsirm/facets.py @@ -0,0 +1,157 @@ +"""Many-Facet Rasch Model (Linacre, 1989): the rating-scale Rasch model with a +rater-severity facet, estimated by marginal-ML EM in the Rust core. All numeric +work happens in Rust; this module only validates and marshals arrays.""" + +from __future__ import annotations + +from dataclasses import dataclass + +import numpy as np + +from .config import MAX_MAX_ITER, MAX_POLYTOMOUS_CATEGORIES + + +@dataclass +class FacetsFit: + """Fitted many-facet Rasch model (Linacre, 1989). + + ``item_difficulty`` is the per-item ``d_i``; ``rater_severity`` the per-rater + ``c_j`` (centered to sum 0; higher = harsher); ``thresholds`` the ``n_cat-1`` + common category thresholds (centered to sum 0); ``theta`` the per-person EAP + trait. The adjacent-category log-odds are + ``ln[P(k)/P(k-1)] = theta - d_i - c_j - f_k``. ``connected`` is False when the + item-rater co-observation design splits into disconnected components — then + severity/difficulty comparisons across components rest solely on the shared + ``theta ~ N(0,1)`` assumption rather than on the rating design (Linacre's + connectedness requirement).""" + + item_difficulty: np.ndarray + rater_severity: np.ndarray + thresholds: np.ndarray + theta: np.ndarray + loglik_trace: np.ndarray + n_iter: int + converged: bool + connected: bool + n_parameters: int + + +def fit_facets( + responses: np.ndarray, + n_cat: int | None = None, + q_theta: int = 41, + max_iter: int = 500, + tol: float = 1e-6, +) -> FacetsFit: + """Fit the many-facet Rasch model (compute in Rust; Linacre, 1989). + + The MFRM extends the rating scale model (Andrich, 1978) with a rater facet: + the rating of person ``p`` on item ``i`` by rater ``j`` follows the + adjacent-category log-odds + ``ln[P(Y=k)/P(Y=k-1)] = theta_p - d_i - c_j - f_k``, where ``d_i`` is item + difficulty, ``c_j`` rater severity, and ``f_k`` the category thresholds + shared across items and raters. ``theta ~ N(0,1)`` fixes the scale; + severities and thresholds are centered to sum to zero. Estimation is + marginal-ML EM (Bock & Aitkin, 1981) on a Gauss-Hermite trait grid — not + Linacre's JMLE, so estimates match Facets output only up to the JMLE-vs-MMLE + difference. + + In LLM-as-a-Judge calibration, raters are judges: ``rater_severity`` + estimates each judge's harshness on a common logit scale, adjusted for item + difficulty and respondent ability. + + ``responses`` is a ``persons x items x raters`` array of integer category + indices ``0..n_cat-1``; ``NaN`` or negative sentinels mark unscored cells + (sparse judging plans), dropped under a missing-at-random assumption. + ``n_cat`` defaults to ``max(responses) + 1``. Every item and every rater + needs at least one observed rating. + + References (APA 7th ed.): + Linacre, J. M. (1989). *Many-facet Rasch measurement*. MESA Press. + Eckes, T. (2015). *Introduction to many-facet Rasch measurement* + (2nd ed.). Peter Lang. https://doi.org/10.3726/978-3-653-04844-5 + Bock, R. D., & Aitkin, M. (1981). Marginal maximum likelihood estimation + of item parameters: Application of an EM algorithm. *Psychometrika, + 46*(4), 443-459. https://doi.org/10.1007/BF02293801 + Andrich, D. (1978). A rating formulation for ordered response + categories. *Psychometrika, 43*(4), 561-573. + https://doi.org/10.1007/BF02293814 + """ + from .fitstats import _core_module + + core = _core_module() + if core is None or not hasattr(core, "fit_facets"): + raise RuntimeError("fit_facets requires the compiled Rust core") + + if not isinstance(n_cat, (int, type(None))) or isinstance(n_cat, bool): + raise ValueError("n_cat must be an integer >= 2") + if n_cat is not None and not (2 <= n_cat <= MAX_POLYTOMOUS_CATEGORIES): + raise ValueError(f"n_cat must be an integer in 2..{MAX_POLYTOMOUS_CATEGORIES}") + if q_theta not in {7, 11, 15, 21, 31, 41}: + raise ValueError("q_theta must be one of 7, 11, 15, 21, 31, 41") + if ( + not isinstance(max_iter, int) + or isinstance(max_iter, bool) + or not (1 <= max_iter <= MAX_MAX_ITER) + ): + raise ValueError(f"max_iter must be an integer in 1..{MAX_MAX_ITER}") + if not np.isfinite(tol) or tol <= 0: + raise ValueError("tol must be finite and > 0") + + y = np.asarray(responses, dtype=np.float64) + if y.ndim != 3: + raise ValueError("responses must be a 3-D persons x items x raters array") + n_persons, n_items, n_raters = y.shape + if n_persons < 1 or n_items < 1 or n_raters < 1: + raise ValueError( + "responses must contain at least one person, one item and one rater" + ) + if np.any(np.isinf(y)): + raise ValueError("observed responses must be finite integer categories") + observed = np.isfinite(y) & (y >= 0) + obs_values = y[observed] + if obs_values.size and np.any(obs_values != np.floor(obs_values)): + raise ValueError("observed responses must be integer categories") + if n_cat is None: + if obs_values.size == 0: + raise ValueError("responses has no observed values") + n_cat = int(obs_values.max()) + 1 + if n_cat < 2: + raise ValueError("responses must contain at least two categories") + if n_cat > MAX_POLYTOMOUS_CATEGORIES: + raise ValueError( + f"responses imply more than {MAX_POLYTOMOUS_CATEGORIES} categories" + ) + if obs_values.size and np.any(obs_values >= n_cat): + raise ValueError( + f"observed responses must be integer categories in 0..{n_cat - 1}" + ) + missing_items = np.flatnonzero(~observed.any(axis=(0, 2))) + if missing_items.size: + raise ValueError(f"item {int(missing_items[0])} has no observed responses") + missing_raters = np.flatnonzero(~observed.any(axis=(0, 1))) + if missing_raters.size: + raise ValueError(f"rater {int(missing_raters[0])} has no observed responses") + yy = np.where(observed, y, 0.0).astype(np.int64).reshape(-1) + res = core.fit_facets( + yy, + observed.reshape(-1), + int(n_persons), + int(n_items), + int(n_raters), + int(n_cat), + int(q_theta), + int(max_iter), + float(tol), + ) + return FacetsFit( + item_difficulty=np.asarray(res["item_difficulty"], dtype=np.float64), + rater_severity=np.asarray(res["rater_severity"], dtype=np.float64), + thresholds=np.asarray(res["thresholds"], dtype=np.float64), + theta=np.asarray(res["theta"], dtype=np.float64), + loglik_trace=np.asarray(res["loglik_trace"], dtype=np.float64), + n_iter=int(res["n_iter"]), + converged=bool(res["converged"]), + connected=bool(res["connected"]), + n_parameters=int(res["n_parameters"]), + ) diff --git a/tests/test_paper_features.py b/tests/test_paper_features.py index 6436d2ac4..c0a6b2964 100644 --- a/tests/test_paper_features.py +++ b/tests/test_paper_features.py @@ -4883,3 +4883,128 @@ def test_fit_testlet_recovers_local_dependence(): with pytest.raises(RuntimeError, match="max_iter_reached"): fit_testlet(y[:40], tid, model="rasch", max_iter=1, require_convergence=True) + + +def test_fit_facets_recovers_rater_severity(): + """Many-facet Rasch model (Linacre, 1989): recover asymmetric rater + severities, item difficulties, and shared thresholds from a sparse judging + plan; single-rater case must agree with fit_rsm (RSM reduction). + + Asserts read crate outputs (res.rater_severity / item_difficulty / + thresholds / theta / loglik_trace / n_parameters / connected). A severity + sign-flip or a d/c dimension-map swap in the Rust core fails the recovery + and reduction checks.""" + import numpy as np + import pytest + from fast_mlsirm import fit_facets, fit_rsm, FacetsFit + from fast_mlsirm.fitstats import _core_module + + core = _core_module() + if core is None or not hasattr(core, "fit_facets"): + pytest.skip("compiled core built without fit_facets") + + rng = np.random.default_rng(1989) + n, n_items, n_raters, n_cat = 1500, 8, 4, 4 + d_true = -0.9 + 0.25 * np.arange(n_items) + c_true = np.array([1.1, -0.2, -0.4, -0.5]) # asymmetric, sums to 0 + f_true = np.array([0.7, 0.1, -0.8]) # sums to 0 + theta = rng.standard_normal(n) + tk = np.concatenate([[0.0], np.cumsum(f_true)]) + ks = np.arange(n_cat) + + y = np.full((n, n_items, n_raters), np.nan) + for p in range(n): + for i in range(n_items): + for j in range(n_raters): + if rng.random() < 0.4: # sparse plan + continue + psi = ks * theta[p] - ks * (d_true[i] + c_true[j]) - tk + pr = np.exp(psi - psi.max()) + y[p, i, j] = rng.choice(n_cat, p=pr / pr.sum()) + + res = fit_facets(y, n_cat=n_cat) + assert isinstance(res, FacetsFit) and res.converged and res.connected + assert np.all(np.diff(res.loglik_trace) >= -1e-6) + assert res.n_parameters == n_items + (n_raters - 1) + (n_cat - 2) + assert abs(res.rater_severity.sum()) < 1e-6 + assert abs(res.thresholds.sum()) < 1e-6 + assert np.sqrt(np.mean((res.rater_severity - c_true) ** 2)) < 0.12 + assert np.sqrt(np.mean((res.item_difficulty - d_true) ** 2)) < 0.15 + assert np.sqrt(np.mean((res.thresholds - f_true) ** 2)) < 0.12 + assert np.corrcoef(res.theta, theta)[0, 1] > 0.85 + + # single-rater reduction: MFRM with J=1 must match RSM (severity absorbed) + y1 = y[:400, :, :1] + keep = ~np.isnan(y1).all(axis=(1, 2)) + y1 = y1[keep] + r_f = fit_facets(y1, n_cat=n_cat) + r_r = fit_rsm(y1[:, :, 0], n_cat=n_cat) + assert np.allclose(r_f.rater_severity, [0.0]) + assert np.allclose(r_f.item_difficulty, r_r.item_location, atol=5e-3) + assert np.allclose(r_f.thresholds, r_r.thresholds, atol=5e-3) + + +def test_fit_facets_rejects_malformed_and_flags_disconnected(): + """MFRM input validation plus Linacre's connectedness diagnostic: a judging + plan whose item-rater graph splits into components must set connected=False + (asserts read res.connected from the crate). The False assert kills a + deleted union-find pass; the True assert on the bridged plan kills the + join-only-item-item-edges mutant, which would leave rater nodes isolated + and report every design as disconnected.""" + import numpy as np + import pytest + from fast_mlsirm import fit_facets + from fast_mlsirm.fitstats import _core_module + + core = _core_module() + if core is None or not hasattr(core, "fit_facets"): + pytest.skip("compiled core built without fit_facets") + + valid = np.zeros((4, 2, 2)) + valid[::2] = 1.0 + with pytest.raises(ValueError, match="3-D"): + fit_facets(np.zeros((4, 2)), n_cat=2) + with pytest.raises(ValueError, match="at least one person"): + fit_facets(np.empty((0, 2, 2)), n_cat=2) + with pytest.raises(ValueError, match="integer categories"): + fit_facets(valid + 0.5, n_cat=2) + with pytest.raises(ValueError, match="in 0..1"): + fit_facets(valid * 3, n_cat=2) + nan_missing = valid.copy() + nan_missing[0, 0, 0] = np.nan + neg_missing = valid.copy() + neg_missing[0, 0, 0] = -1 + res_nan = fit_facets(nan_missing, n_cat=2, q_theta=7, max_iter=5) + res_neg = fit_facets(neg_missing, n_cat=2, q_theta=7, max_iter=5) + assert np.allclose(res_neg.item_difficulty, res_nan.item_difficulty) + assert np.allclose(res_neg.rater_severity, res_nan.rater_severity) + assert np.allclose(res_neg.thresholds, res_nan.thresholds) + assert np.allclose(res_neg.theta, res_nan.theta) + with pytest.raises(ValueError, match="rater 1 has no observed"): + bad = valid.copy() + bad[:, :, 1] = np.nan + fit_facets(bad, n_cat=2) + with pytest.raises(ValueError, match="q_theta"): + fit_facets(valid, n_cat=2, q_theta=10) + with pytest.raises(ValueError, match="max_iter"): + fit_facets(valid, n_cat=2, max_iter=0) + with pytest.raises(ValueError, match="tol"): + fit_facets(valid, n_cat=2, tol=0.0) + + # disconnected plan: persons 0-19 see (item0, rater0), persons 20-39 see + # (item1, rater1) -- no shared element links the two components + rng = np.random.default_rng(7) + y = np.full((40, 2, 2), np.nan) + y[:20, 0, 0] = rng.integers(0, 2, 20) + y[20:, 1, 1] = rng.integers(0, 2, 20) + y[0, 0, 0], y[1, 0, 0] = 0.0, 1.0 + y[20, 1, 1], y[21, 1, 1] = 0.0, 1.0 + res = fit_facets(y, n_cat=2, max_iter=50) + assert res.connected is False + + # bridged plan: person 5 also sees (item1, rater1), joining the components + yb = y.copy() + yb[5, 1, 1] = 1.0 + yb[6, 1, 1] = 0.0 + resb = fit_facets(yb, n_cat=2, max_iter=50) + assert resb.connected is True diff --git a/tests/unit/facets_tests.rs b/tests/unit/facets_tests.rs new file mode 100644 index 000000000..e1a81a657 --- /dev/null +++ b/tests/unit/facets_tests.rs @@ -0,0 +1,365 @@ +use super::*; + +struct Lcg(u64); +impl Lcg { + fn f64(&mut self) -> f64 { + self.0 = self + .0 + .wrapping_mul(6364136223846793005) + .wrapping_add(1442695040888963407); + ((self.0 >> 11) as f64) / ((1u64 << 53) as f64) + } + fn normal(&mut self) -> f64 { + let u1 = self.f64().max(1e-12); + let u2 = self.f64(); + (-2.0 * u1.ln()).sqrt() * (2.0 * std::f64::consts::PI * u2).cos() + } +} + +fn rmse(a: &[f64], b: &[f64]) -> f64 { + (a.iter().zip(b).map(|(x, y)| (x - y).powi(2)).sum::() / a.len() as f64).sqrt() +} +fn corr(x: &[f64], y: &[f64]) -> f64 { + let n = x.len() as f64; + let mx = x.iter().sum::() / n; + let my = y.iter().sum::() / n; + let (mut sxy, mut sxx, mut syy) = (0.0, 0.0, 0.0); + for i in 0..x.len() { + sxy += (x[i] - mx) * (y[i] - my); + sxx += (x[i] - mx).powi(2); + syy += (y[i] - my).powi(2); + } + sxy / (sxx.sqrt() * syy.sqrt()) +} + +/// Draw an MFRM category for ability `theta`, item `d`, rater `c`, thresholds `f`. +fn draw_facets(theta: f64, d: f64, c: f64, f: &[f64], u: f64) -> usize { + let lp = crate::rsm::rsm_logprobs(theta, d + c, f); + let mut cum = 0.0; + for (k, l) in lp.iter().enumerate() { + cum += l.exp(); + if u < cum { + return k; + } + } + lp.len() - 1 +} + +/// Simulate a fully crossed design. Returns row-major `P*I*J` categories. +fn simulate( + seed: u64, + n_persons: usize, + d: &[f64], + c: &[f64], + f: &[f64], +) -> Vec { + let mut rng = Lcg(seed); + let (ni, nj) = (d.len(), c.len()); + let mut y = vec![0usize; n_persons * ni * nj]; + for p in 0..n_persons { + let theta = rng.normal(); + for i in 0..ni { + for j in 0..nj { + y[p * ni * nj + i * nj + j] = draw_facets(theta, d[i], c[j], f, rng.f64()); + } + } + } + y +} + +// --------------------------------------------------------------------------- +// FD anchors on the M-step objective. These asserts read the crate's +// `location_score_terms` / `f_gradient` outputs and compare them against +// central finite differences of the crate's `pair_ell`/`total_ell` at an +// ASYMMETRIC point. Mutations killed: sign flips in the gradients, d<->c index +// transposition (the location derivative would hit the wrong count block), +// suffix-sum off-by-one in `f_gradient` (shifts which residuals feed g_m). +// --------------------------------------------------------------------------- + +#[test] +fn location_gradient_matches_fd() { + // Asymmetric params and asymmetric fake counts (not a fitted state). + let f = [0.9f64, -0.4, -0.5]; + let n_cat = 4usize; + let nodes = [-1.3f64, 0.2, 1.7]; + let mut r = vec![0.0f64; nodes.len() * n_cat]; + let mut rng = Lcg(7); + for v in r.iter_mut() { + *v = 0.05 + rng.f64() * 2.0; + } + let loc0 = 0.37f64; + let (mut g, mut h) = (0.0f64, 0.0f64); + location_score_terms(loc0, &f, &r, &nodes, n_cat, &mut g, &mut h); + let eps = 1e-6; + let fd = (pair_ell(loc0 + eps, &f, &r, &nodes, n_cat) + - pair_ell(loc0 - eps, &f, &r, &nodes, n_cat)) + / (2.0 * eps); + // location_score_terms accumulates -d ell/d location... verify sign + // convention explicitly: Newton uses step = g/h with update loc - step, + // and the code defines g = -sum k (r - nP) = -d ell/d loc? No: the + // derivative d ell/d loc = -sum_k k (r - n P) exactly, so g == fd. + assert!( + (g - fd).abs() < 1e-5, + "analytic {g} vs FD {fd}" + ); + assert!(h < 0.0, "location Hessian must be negative, got {h}"); + // Hessian FD check too (kills Var-of-score sign/formula mutations). + let (mut gp, mut hp) = (0.0f64, 0.0f64); + location_score_terms(loc0 + eps, &f, &r, &nodes, n_cat, &mut gp, &mut hp); + let (mut gm, mut hm) = (0.0f64, 0.0f64); + location_score_terms(loc0 - eps, &f, &r, &nodes, n_cat, &mut gm, &mut hm); + let fd_h = (gp - gm) / (2.0 * eps); + assert!((h - fd_h).abs() < 1e-4, "analytic H {h} vs FD {fd_h}"); +} + +#[test] +fn threshold_gradient_matches_fd() { + let d = [0.3f64, -0.8]; + let c = [0.5f64, -0.1, -0.4]; + let f = [0.7f64, -0.2, -0.5]; + let n_cat = 4usize; + let nodes = [-1.1f64, 0.4, 2.0]; + let n_pairs = d.len() * c.len(); + let mut rng = Lcg(11); + let mut r = vec![vec![0.0f64; nodes.len() * n_cat]; n_pairs]; + for blk in r.iter_mut() { + for v in blk.iter_mut() { + *v = 0.05 + rng.f64(); + } + } + let g = f_gradient(&f, &d, &c, &r, &nodes, d.len(), c.len(), n_cat); + let eps = 1e-6; + for m in 0..f.len() { + let mut fp = f.to_vec(); + fp[m] += eps; + let mut fm = f.to_vec(); + fm[m] -= eps; + let fd = (total_ell(&d, &c, &fp, &r, &nodes, d.len(), c.len(), n_cat) + - total_ell(&d, &c, &fm, &r, &nodes, d.len(), c.len(), n_cat)) + / (2.0 * eps); + assert!( + (g[m] - fd).abs() < 1e-5, + "f[{m}]: analytic {} vs FD {fd}", + g[m] + ); + } +} + +// --------------------------------------------------------------------------- +// J=1 reduction anchor: with one rater the MFRM must reproduce fit_rsm. +// Asserts read fit_facets' item_difficulty/thresholds/loglik and fit_rsm's +// outputs. Mutations killed: wrong aggregation over the rater axis, pair +// indexing bugs (i*n_raters+j vs j*n_items+i), severity leaking into the fit. +// --------------------------------------------------------------------------- +#[test] +fn single_rater_reduces_to_rsm() { + let d_true = [-1.2f64, -0.3, 0.4, 1.1]; + let f_true = [0.8f64, -0.8]; + let y = simulate(42, 400, &d_true, &[0.0], &f_true); + let res = fit_facets(&y, None, 400, 4, 1, 3, 21, 300, 1e-8).unwrap(); + let rsm = crate::rsm::fit_rsm(&y, None, 400, 4, 3, 21, 300, 1e-8).unwrap(); + // sum(c)=0 with one rater forces c_1 = 0 exactly. + assert!(res.rater_severity[0].abs() < 1e-12); + for i in 0..4 { + assert!( + (res.item_difficulty[i] - rsm.item_location[i]).abs() < 1e-4, + "item {i}: facets {} vs rsm {}", + res.item_difficulty[i], + rsm.item_location[i] + ); + } + for m in 0..2 { + assert!((res.thresholds[m] - rsm.thresholds[m]).abs() < 1e-4); + } + let lf = *res.loglik_trace.last().unwrap(); + let lr = *rsm.loglik_trace.last().unwrap(); + assert!((lf - lr).abs() < 1e-4, "loglik facets {lf} vs rsm {lr}"); + assert!(res.connected); + assert_eq!(res.n_parameters, 4 + 0 + 1); +} + +// --------------------------------------------------------------------------- +// Severity recovery with an asymmetric severity vector. Asserts read +// res.rater_severity (crate output). Mutations killed: over-collapse (all +// severities shrink to ~0 -> corr undefined/rmse large), sign flip in the c +// update (corr ~ -1), rater/item dimension-map swap (J=5 != I=6 so shapes +// diverge and recovery fails). +// --------------------------------------------------------------------------- +#[test] +fn recovers_asymmetric_rater_severity() { + let d_true = [-1.5f64, -0.9, -0.2, 0.3, 0.9, 1.6]; + let c_true = [-1.0f64, -0.3, 0.1, 0.4, 0.8]; // deliberately not centered + let f_true = [1.0f64, 0.1, -1.1]; + let y = simulate(2024, 800, &d_true, &c_true, &f_true); + let res = fit_facets(&y, None, 800, 6, 5, 4, 21, 500, 1e-8).unwrap(); + assert!(res.converged); + assert!(res.connected); + // Compare against the centered generating severities (model identifies c + // only up to the sum-zero constraint; the mean shift moves into d). + let mean_c = c_true.iter().sum::() / c_true.len() as f64; + let c_centered: Vec = c_true.iter().map(|v| v - mean_c).collect(); + let r = corr(&res.rater_severity, &c_centered); + let e = rmse(&res.rater_severity, &c_centered); + assert!(r > 0.95, "severity corr {r}"); + assert!(e < 0.15, "severity rmse {e}"); + // Item difficulty absorbs the shift: d_hat ~ d_true + mean_c (+ f-centering + // shift, which is 0 here up to sampling because f_true sums to 0). + let d_shifted: Vec = d_true.iter().map(|v| v + mean_c).collect(); + let rd = corr(&res.item_difficulty, &d_shifted); + assert!(rd > 0.95, "difficulty corr {rd}"); + // Structural invariants of the returned parameters (not test-local math): + // both centerings hold on the crate output. + let sum_c: f64 = res.rater_severity.iter().sum(); + let sum_f: f64 = res.thresholds.iter().sum(); + assert!(sum_c.abs() < 1e-9, "sum(c) = {sum_c}"); + assert!(sum_f.abs() < 1e-9, "sum(f) = {sum_f}"); + assert_eq!(res.n_parameters, 6 + 4 + 2); + // Known limitation: a constant-shift mutation applied jointly to d and -c + // is a model invariance and cannot be detected by any data-based test; + // the discriminating anchors are the centering asserts above. +} + +// --------------------------------------------------------------------------- +// Sparse judging plan: each person is scored by 2 of 5 raters on a rotating +// (non-contiguous) schedule. Asserts read crate outputs. Mutations killed: +// dense-only indexing (missing cells would feed category 0 counts), observed- +// mask offset bugs. +// --------------------------------------------------------------------------- +#[test] +fn sparse_design_recovers_severity_order() { + let d_true = [-0.8f64, 0.0, 0.8]; + let c_true = [-0.9f64, -0.2, 0.0, 0.3, 0.8]; + let f_true = [0.6f64, -0.6]; + let (np, ni, nj) = (1500usize, 3usize, 5usize); + let y = simulate(99, np, &d_true, &c_true, &f_true); + // Rotating pairs (p, p+2 mod 5): non-contiguous rater unions, connected. + let mut obs = vec![false; np * ni * nj]; + for p in 0..np { + let (a, b) = (p % nj, (p + 2) % nj); + for i in 0..ni { + obs[p * ni * nj + i * nj + a] = true; + obs[p * ni * nj + i * nj + b] = true; + } + } + let res = fit_facets(&y, Some(&obs), np, ni, nj, 3, 21, 500, 1e-8).unwrap(); + assert!(res.connected); + let mean_c = c_true.iter().sum::() / nj as f64; + let c_centered: Vec = c_true.iter().map(|v| v - mean_c).collect(); + let r = corr(&res.rater_severity, &c_centered); + assert!(r > 0.9, "sparse severity corr {r}"); + // The recovered severity ORDER must match (kills permutation/off-by-one + // in the rater axis under a sparse mask). + let mut idx: Vec = (0..nj).collect(); + idx.sort_by(|&a, &b| res.rater_severity[a].partial_cmp(&res.rater_severity[b]).unwrap()); + assert_eq!(idx, vec![0, 1, 2, 3, 4]); +} + +// --------------------------------------------------------------------------- +// Connectivity flag. Asserts read res.connected. Mutations killed: joining +// only items to items (not raters), skipping the person-mediated union, or +// hardcoding true. +// --------------------------------------------------------------------------- +#[test] +fn disconnected_design_is_flagged() { + // Two islands: persons 0..P/2 x item 0 x rater 0; persons P/2.. x item 1 x rater 1. + let (np, ni, nj) = (60usize, 2usize, 2usize); + let d = [0.0f64, 0.0]; + let c = [0.5f64, -0.5]; + let f = [0.0f64]; + let y = simulate(5, np, &d, &c, &f); + let mut obs = vec![false; np * ni * nj]; + for p in 0..np { + let island = usize::from(p >= np / 2); + obs[p * ni * nj + island * nj + island] = true; + } + let res = fit_facets(&y, Some(&obs), np, ni, nj, 2, 7, 50, 1e-6).unwrap(); + assert!(!res.connected, "two islands must be flagged disconnected"); + + // Bridging rater: rater 0 also scores item 1 for one person -> connected. + let mut obs2 = obs.clone(); + obs2[0 * ni * nj + 1 * nj + 0] = true; // person 0, item 1, rater 0 + let res2 = fit_facets(&y, Some(&obs2), np, ni, nj, 2, 7, 50, 1e-6).unwrap(); + assert!(res2.connected, "bridge must connect the design"); +} + +// --------------------------------------------------------------------------- +// Validation errors. +// --------------------------------------------------------------------------- +#[test] +fn rejects_bad_inputs() { + let y = vec![0usize; 4]; + assert!(fit_facets(&y, None, 2, 2, 1, 1, 7, 50, 1e-6).is_err()); // n_cat < 2 + assert!(fit_facets(&y, None, 2, 2, 1, 2, 8, 50, 1e-6).is_err()); // bad q + assert!(fit_facets(&y, None, 2, 2, 1, 2, 7, 0, 1e-6).is_err()); // max_iter 0 + assert!(fit_facets(&y, None, 2, 2, 1, 2, 7, 50, f64::NAN).is_err()); + assert!(fit_facets(&y, None, 3, 2, 1, 2, 7, 50, 1e-6).is_err()); // len mismatch + let y2 = vec![0usize, 5, 0, 0]; + assert!(fit_facets(&y2, None, 2, 2, 1, 3, 7, 50, 1e-6).is_err()); // cat >= n_cat + // rater with no observations + let y3 = vec![0usize; 2 * 1 * 2]; + let obs = vec![true, false, true, false]; + assert!(fit_facets(&y3, Some(&obs), 2, 1, 2, 2, 7, 50, 1e-6) + .unwrap_err() + .contains("rater 1")); +} + +#[test] +fn loglik_trace_is_nondecreasing() { + let y = simulate(3, 200, &[-0.5, 0.5], &[-0.4, 0.4], &[0.5, -0.5]); + let res = fit_facets(&y, None, 200, 2, 2, 3, 21, 200, 1e-10).unwrap(); + for w in res.loglik_trace.windows(2) { + assert!(w[1] >= w[0] - 1e-8, "EM must be monotone: {} -> {}", w[0], w[1]); + } +} + +// --------------------------------------------------------------------------- +// Monte-Carlo recovery, 500 replications (heavy; run with --ignored). +// Half the replications generate theta from a skewed distribution (mixture +// shift) to probe prior-misspecification robustness: bias is reported with a +// loose bound rather than asserted tightly. +// --------------------------------------------------------------------------- +#[test] +#[ignore = "500-replication Monte-Carlo; run with --ignored"] +fn monte_carlo_severity_bias_and_rmse() { + let d_true = [-1.0f64, 0.0, 1.0]; + let c_true = [-0.7f64, 0.0, 0.7]; // centered + let f_true = [0.9f64, -0.9]; + let (np, ni, nj) = (300usize, 3usize, 3usize); + let reps = 500usize; + let mut bias = vec![0.0f64; nj]; + let mut mse = vec![0.0f64; nj]; + for rep in 0..reps { + let skewed = rep % 2 == 1; + let mut rng = Lcg(10_000 + rep as u64); + let mut y = vec![0usize; np * ni * nj]; + for p in 0..np { + let theta = if skewed { + // Standardized two-component location mixture (negatively + // skewed), mean 0 / var ~1 by construction below. + let z = rng.normal(); + let comp = if rng.f64() < 0.75 { 0.35 } else { -1.05 }; + (z * 0.8 + comp) / (0.8f64.powi(2) + 0.42f64).sqrt() + } else { + rng.normal() + }; + for i in 0..ni { + for j in 0..nj { + y[p * ni * nj + i * nj + j] = + draw_facets(theta, d_true[i], c_true[j], &f_true, rng.f64()); + } + } + } + let res = fit_facets(&y, None, np, ni, nj, 3, 21, 500, 1e-8).unwrap(); + for j in 0..nj { + let e = res.rater_severity[j] - c_true[j]; + bias[j] += e / reps as f64; + mse[j] += e * e / reps as f64; + } + } + for j in 0..nj { + let rm = mse[j].sqrt(); + // Loose bounds: severity is a fixed effect over 300*3 ratings/rep. + assert!(bias[j].abs() < 0.05, "rater {j} bias {}", bias[j]); + assert!(rm < 0.2, "rater {j} rmse {rm}"); + } +}