feat(facets): many-facet Rasch model rater-severity calibration (Linacre, 1989)

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

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<Py<pyo3::types::PyDict>> {
+    let yy: Vec<usize> = 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::<PyResult<_>>()?;
+    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)?)?;