From 8c22190ca67e30259091f7b4a96e1e7ebb8545bb Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Mon, 27 Jul 2026 09:55:54 +0900 Subject: [PATCH 1/9] Add fleiss_kappa: Fleiss' multi-rater kappa with exact (Conger) variant Reimplements CRAN irr 0.85 kappam.fleiss() (R source READ and normative; Fleiss 1971 and Conger 1980 NOT READ, cited as model origins only) in the Rust core (mlsirm_core::agreement::fleiss_kappa) with a thin PyO3 binding and NumPy wrapper: - classification-table agreement, classic (sum p_j^2) and exact (sum p_j^2 - (1/nr) sum s2_j) chance agreement; kappa, Fleiss' z test, and category-wise kappas (classic mode; NaN for empty categories, matching R's 0/0) - listwise row drop for missing ratings (negative code / NaN) - documented API deviations: index codes 0..k-1 with explicit/inferred k, error on degenerate 1 - chanceP = 0 (R returns NaN), size caps Evidence: exact-Fraction oracle anchors FK1-FK5 (classic kappa 139/399, exact 37/102, category kappas [1/21, 31/91, 43/63]); 6 mutants EXECUTED and killed (agreeP centering, row-vs-column chance sums, exact==classic, variance sign, missing-as-category, pjk centering); MC-500 subject/rater permutation-invariance test (#[ignore], executed); cargo 847 pass; TestFleiss pytest pass. Co-authored-by: Copilot App <223556219+Copilot@users.noreply.github.com> --- CHANGELOG.md | 1 + crates/fast-mlsirm-py/src/lib.rs | 30 ++++ crates/mlsirm-core/src/agreement.rs | 192 +++++++++++++++++++++++ python/fast_mlsirm/__init__.py | 6 +- python/fast_mlsirm/validation.py | 86 ++++++++++ tests/test_paper_features.py | 75 +++++++++ tests/unit/agreement_tests.rs | 235 ++++++++++++++++++++++++++++ 7 files changed, 624 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 0818cef36..acc6f377a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -119,6 +119,7 @@ ### Added +- `fleiss_kappa` Fleiss' multi-rater kappa for nominal agreement with the exact (Conger) chance-agreement variant (CRAN irr 0.85 `kappam.fleiss()` `R/kappam.fleiss.R` — READ and normative; Fleiss 1971 and Conger 1980 NOT READ, cited as model origins only): classification-table agreement `agreeP = (1/m)Σᵢ(Σⱼtᵢⱼ² − nr)/(nr(nr−1))`, classic chance `Σⱼpⱼ²` vs exact `Σⱼpⱼ² − (1/nr)Σⱼs²ⱼ` (sample variance over per-rater category proportions; algebra verified against R's `apply(rtab,2,var)` form), Fleiss' large-sample z test and category-wise kappas (classic mode; NaN for empty categories, matching R's 0/0), and listwise row drop for missing ratings. API deviations documented: index codes 0..k-1 with explicit/inferred k, negative-or-NaN = missing, error on degenerate `1 − chanceP = 0` (R returns NaN). Exact-Fraction oracle anchors FK1–FK5 (classic κ = 139/399, exact κ = 37/102, category κ = [1/21, 31/91, 43/63]), a 6-mutant EXECUTED kill map (row-vs-column chance sums, missing-as-category, variance-sign, pjk-centering), and an MC-500 subject/rater permutation-invariance test. - `bratt_mm` Bradley-Terry model with ties fitted by MM (VGAM 1.1-14 `bratt()` family, `R/family.categorical.R` — READ and normative; Bradley & Terry 1952 NOT READ, cited as model origin): `P(i>j) = αᵢ/(αᵢ+αⱼ+α₀)`, `P(tie) = α₀/(αᵢ+αⱼ+α₀)` with a hand-derived supporting-hyperplane MM ascent (same pattern as the crate's `bradley_terry_mm`) and a joint reference rescale of α AND α₀ (likelihood-preserving; verified identity `Σ wins + T = Σ n_ij`). This is the additive-α₀ ties model, NOT Rao-Kupper/Davidson (neither read; disambiguation only). Contract: fractional weighted counts accepted, symmetric ties matrix required, tie-free data rejected (use `bradley_terry_mm`; an API contract, not VGAM behavior), zero-win contestants rejected, n capped at 10000 (O(n²) guard). Exact-Fraction oracle anchors B1–B4 (iteration-1 pins `[1, 27/40, 3/4]`, α₀ = 9/14), a 5-mutant EXECUTED kill map (incl. a tol-separated convergence anchor killing an α₀-blind convergence check), and an MC-500 log-likelihood dominance test. - `predict_rating` / `predict_rating_multi` game-outcome prediction from fitted ratings (CRAN PlayerRatings 1.1-0 `predict.rating` `R/ratings.R` lines 1056–1133 — READ and normative; no journal paper exists for this dispatch, CRAN package provenance only): Elo logistic branch, deviation-shrunk Glicko/Glicko-2/Stephenson branch (`qip3 = 3(ln10/400π)²`, joint shrink over BOTH players' squared deviations), and multi-player EloM branch (`(rating − rowmean)/40`, na.rm row means, optional min-tie placing ranks with NaN kept). R semantics preserved: strict `games < tng` unrated cutoff, `trat` replacement of ALL missing extracted values (unmatched, low-games, stored-NA), `pred >= thresh` binarization with NaN propagation. REDUCED-SCOPE vs R: index-based (−1 = unmatched; caller does name matching), per-game/scalar gamma only. Exact-oracle fixtures P1–P9 and a 7-mutant EXECUTED kill map (incl. both branches' tng comparisons). - `fide_rating` FIDE-style Elo ratings (CRAN PlayerRatings 1.1-0 `fide()` `R/ratings.R` lines 125–272 + `kfide()` lines 959–972 — READ and normative; no journal paper exists for this variant, CRAN package provenance only): per-period batch Elo with the kfide K-factor schedule (K = kv[0] elite / kv[1] ≥30 games / kv[2] novice, evaluated from PERIOD-START state), sticky elite flag set from POST-update ratings ≥ 2400, and per-player running mean of POST-update opponent ratings. REDUCED-SCOPE vs R: no status/history frames, kfide-only K schedule, self-play rejected, thresholds 30/2400 hard-coded. kv=(k,k,k) reduces bitwise to `elo_rating(kfac=k)` (MC-500 anchor); exact-oracle fixtures F1–F5 and a 5-mutant EXECUTED kill map. diff --git a/crates/fast-mlsirm-py/src/lib.rs b/crates/fast-mlsirm-py/src/lib.rs index 9a2f4c617..efd3f0c6e 100644 --- a/crates/fast-mlsirm-py/src/lib.rs +++ b/crates/fast-mlsirm-py/src/lib.rs @@ -5848,6 +5848,35 @@ fn bratt_mm( Ok(d.into()) } +/// Fleiss' kappa for nominal agreement among nr raters over ns subjects, +/// with the exact (Conger) variant (irr 0.85 `kappam.fleiss`; see +/// `mlsirm_core::agreement::fleiss_kappa`). `ratings` is flat row-major +/// ns*nr of codes 0..k-1; negative = missing (listwise row drop). Returns +/// dict with kappa, subjects_used, z, p_value, category_kappa/z/p +/// (empty arrays and NaN z/p in exact mode). +#[pyfunction] +#[pyo3(signature = (ratings, ns, nr, k, exact=false))] +fn fleiss_kappa( + py: Python<'_>, + ratings: PyReadonlyArray1<'_, i64>, + ns: usize, + nr: usize, + k: usize, + exact: bool, +) -> PyResult> { + let res = mlsirm_core::agreement::fleiss_kappa(ratings.as_slice()?, ns, nr, k, exact) + .map_err(PyValueError::new_err)?; + let d = pyo3::types::PyDict::new(py); + d.set_item("kappa", res.kappa)?; + d.set_item("subjects_used", res.subjects_used as u64)?; + d.set_item("z", res.z)?; + d.set_item("p_value", res.p_value)?; + d.set_item("category_kappa", PyArray1::from_slice(py, &res.category_kappa))?; + d.set_item("category_z", PyArray1::from_slice(py, &res.category_z))?; + d.set_item("category_p", PyArray1::from_slice(py, &res.category_p))?; + Ok(d.into()) +} + /// GPCM/nominal softmax cell log-probabilities at one node (parity surface for /// the NumPy `category_logprobs` reference). #[pyfunction] @@ -8205,6 +8234,7 @@ fn fast_mlsirm_core(m: &Bound<'_, PyModule>) -> PyResult<()> { m.add_function(wrap_pyfunction!(predict_rating_two, m)?)?; m.add_function(wrap_pyfunction!(predict_rating_multi, m)?)?; m.add_function(wrap_pyfunction!(bratt_mm, m)?)?; + m.add_function(wrap_pyfunction!(fleiss_kappa, m)?)?; m.add_function(wrap_pyfunction!(circle_arc_middle_anchor, m)?)?; m.add_function(wrap_pyfunction!(loglinear_smooth, m)?)?; m.add_function(wrap_pyfunction!(person_fit_stat, m)?)?; diff --git a/crates/mlsirm-core/src/agreement.rs b/crates/mlsirm-core/src/agreement.rs index 3b88a1de0..c075bc58f 100644 --- a/crates/mlsirm-core/src/agreement.rs +++ b/crates/mlsirm-core/src/agreement.rs @@ -241,6 +241,198 @@ pub fn validate_scoring( }) } +/// Result of Fleiss' multi-rater kappa (`fleiss_kappa`). +/// +/// In exact (Conger) mode `z`/`p_value` are NaN and the category vectors +/// are empty, mirroring irr's `kappam.fleiss(exact=TRUE)` which returns +/// neither a test statistic nor category detail. +#[derive(Clone, Debug)] +pub struct FleissKappaResult { + pub kappa: f64, + /// Subjects remaining after listwise deletion of rows with missing codes. + pub subjects_used: usize, + pub z: f64, + pub p_value: f64, + /// Category-wise kappas (classic mode only); NaN for empty categories. + pub category_kappa: Vec, + pub category_z: Vec, + pub category_p: Vec, +} + +/// Fleiss' kappa for nominal agreement among `nr` raters over `ns` subjects, +/// with the exact (Conger) chance-agreement variant. +/// +/// Reimplements `kappam.fleiss()` from CRAN irr 0.85 (`R/kappam.fleiss.R`, +/// READ in full; algorithm source of truth). The model originates in +/// Fleiss, J. L. (1971), "Measuring nominal scale agreement among many +/// raters," Psychological Bulletin, 76(5), 378-382, and the exact variant in +/// Conger, A. J. (1980), "Integration and generalization of kappas for +/// multiple raters," Psychological Bulletin, 88(2), 322-328 — both cited as +/// origins only (NOT READ); every formula below was verified against the irr +/// R source. +/// +/// `ratings` is row-major `ns x nr` with category codes `0..k-1`; a negative +/// code marks a missing rating and drops the whole subject row (listwise, as +/// in R: `ratings[apply(is.na(ratings),1,sum)==0,]`). With `m` used subjects +/// and `ttab[i][j]` = raters assigning subject `i` to category `j`: +/// +/// - `agreeP = (1/m) sum_i (sum_j ttab_ij^2 - nr) / (nr(nr-1))` +/// - classic `chanceP = sum_j p_j^2` with `p_j = C_j/(m nr)`, `C_j` column sums +/// - exact `chanceP = sum_j p_j^2 - (1/nr) sum_j s2_j`, `s2_j` the sample +/// variance (divisor `nr-1`) over raters of per-rater category proportions +/// (algebraically equal to R's `sum(apply(rtab,2,var)*(nr-1)/nr)/(nr-1)`) +/// - `kappa = (agreeP - chanceP)/(1 - chanceP)` +/// +/// Classic mode adds Fleiss' large-sample test +/// `var = 2[(sum p_j q_j)^2 - sum p_j q_j (q_j - p_j)] / +/// [(sum p_j q_j)^2 m nr (nr-1)]`, `z = kappa/sqrt(var)`, +/// `p = 2(1 - Phi(|z|))`, and category-wise kappas +/// `pjk_j = (sum_i ttab_ij^2 - m nr p_j)/(m nr (nr-1) p_j)`, +/// `kappa_j = (pjk_j - p_j)/(1 - p_j)`, `var_j = 2/(m nr (nr-1))` +/// (computed unconditionally; identical to irr's `detail=TRUE`). Empty +/// categories yield NaN, matching R's 0/0. +/// +/// API deviations from R (documented contract, not transcription): codes are +/// index-based `0..k-1` with explicit `k` (R derives factor levels; negative +/// numeric labels must be remapped by the caller since negative = missing +/// here), degenerate `1 - chanceP == 0` is an error (R returns NaN), and the +/// size caps below are safety bounds. +pub fn fleiss_kappa( + ratings: &[i64], + ns: usize, + nr: usize, + k: usize, + exact: bool, +) -> Result { + if ns == 0 { + return Err("need at least one subject".into()); + } + if nr < 2 { + return Err("need at least 2 raters".into()); + } + if k < 2 { + return Err("need at least 2 categories".into()); + } + if ns > 1_000_000 || nr > 10_000 || k > 10_000 { + return Err("size caps: ns <= 1e6, nr <= 1e4, k <= 1e4".into()); + } + if ratings.len() != ns * nr { + return Err(format!( + "ratings length {} != ns*nr = {}", + ratings.len(), + ns * nr + )); + } + for &c in ratings { + if c >= k as i64 { + return Err(format!("category code {c} out of range 0..{k}")); + } + } + // Listwise drop of rows containing any negative (missing) code, then + // classification table ttab[i][j] and per-rater counts. Counts are exact + // in f64: entries <= nr <= 1e4, sums of squares <= m*nr^2 <= 1e14 < 2^53. + let mut ttab: Vec> = Vec::new(); + let mut rater_counts = vec![vec![0.0_f64; k]; nr]; + for i in 0..ns { + let row = &ratings[i * nr..(i + 1) * nr]; + if row.iter().any(|&c| c < 0) { + continue; + } + let mut t = vec![0.0_f64; k]; + for (r, &c) in row.iter().enumerate() { + t[c as usize] += 1.0; + rater_counts[r][c as usize] += 1.0; + } + ttab.push(t); + } + let m = ttab.len(); + if m == 0 { + return Err("all subject rows dropped for missing ratings".into()); + } + let mf = m as f64; + let nrf = nr as f64; + + let agree_p: f64 = ttab + .iter() + .map(|t| (t.iter().map(|&v| v * v).sum::() - nrf) / (nrf * (nrf - 1.0))) + .sum::() + / mf; + + let col: Vec = (0..k) + .map(|j| ttab.iter().map(|t| t[j]).sum::()) + .collect(); + let p: Vec = col.iter().map(|&c| c / (mf * nrf)).collect(); + let mut chance_p: f64 = p.iter().map(|&v| v * v).sum(); + if exact { + // Sample variance over raters of rtab[r][j] = rater_counts[r][j]/m. + let mut s2_sum = 0.0_f64; + for j in 0..k { + let props: Vec = (0..nr).map(|r| rater_counts[r][j] / mf).collect(); + let mean = props.iter().sum::() / nrf; + let s2 = props.iter().map(|&v| (v - mean) * (v - mean)).sum::() / (nrf - 1.0); + s2_sum += s2; + } + chance_p -= s2_sum / nrf; + } + let denom = 1.0 - chance_p; + if denom.abs() < 1e-12 { + return Err("degenerate marginals: no chance-corrected agreement is defined".into()); + } + let kappa = (agree_p - chance_p) / denom; + + if exact { + return Ok(FleissKappaResult { + kappa, + subjects_used: m, + z: f64::NAN, + p_value: f64::NAN, + category_kappa: Vec::new(), + category_z: Vec::new(), + category_p: Vec::new(), + }); + } + + let sqrt2 = std::f64::consts::SQRT_2; + let pq: f64 = p.iter().map(|&v| v * (1.0 - v)).sum(); + let var = 2.0 + * (pq * pq + - p.iter() + .map(|&v| v * (1.0 - v) * (1.0 - 2.0 * v)) + .sum::()) + / (pq * pq * mf * nrf * (nrf - 1.0)); + let z = kappa / var.sqrt(); + let p_value = crate::fitstats::erfc(z.abs() / sqrt2); + + let var_k = 2.0 / (mf * nrf * (nrf - 1.0)); + let mut category_kappa = Vec::with_capacity(k); + let mut category_z = Vec::with_capacity(k); + let mut category_p = Vec::with_capacity(k); + for j in 0..k { + let sum_sq: f64 = ttab.iter().map(|t| t[j] * t[j]).sum(); + // Empty category: p_j = 0 gives R's 0/0 = NaN, preserved here. + let pjk = (sum_sq - mf * nrf * p[j]) / (mf * nrf * (nrf - 1.0) * p[j]); + let kj = (pjk - p[j]) / (1.0 - p[j]); + let zj = kj / var_k.sqrt(); + category_kappa.push(kj); + category_z.push(zj); + category_p.push(if zj.is_finite() { + crate::fitstats::erfc(zj.abs() / sqrt2) + } else { + f64::NAN + }); + } + + Ok(FleissKappaResult { + kappa, + subjects_used: m, + z, + p_value, + category_kappa, + category_z, + category_p, + }) +} + #[cfg(test)] #[path = "../../../tests/unit/agreement_tests.rs"] mod tests; diff --git a/python/fast_mlsirm/__init__.py b/python/fast_mlsirm/__init__.py index 00df71b2a..2c7a80105 100644 --- a/python/fast_mlsirm/__init__.py +++ b/python/fast_mlsirm/__init__.py @@ -175,7 +175,9 @@ from .lltm import fit_lltm as fit_lltm, LltmFit as LltmFit from .testlet import fit_testlet as fit_testlet, TestletFit as TestletFit from .report import render_diagnostics_report as render_diagnostics_report -from .validation import (ValidationVerdict as ValidationVerdict, +from .validation import (FleissKappaResult as FleissKappaResult, + ValidationVerdict as ValidationVerdict, + fleiss_kappa as fleiss_kappa, validate_judge as validate_judge) from .serving import (bank_information as bank_information, cat_next_item as cat_next_item, @@ -220,6 +222,8 @@ "assemble_test_form", "dimensionality_diagnostics", "ValidationVerdict", + "FleissKappaResult", + "fleiss_kappa", "benjamini_hochberg", "chi2_sf", "dif_analysis", diff --git a/python/fast_mlsirm/validation.py b/python/fast_mlsirm/validation.py index 61803de66..60dd445c1 100644 --- a/python/fast_mlsirm/validation.py +++ b/python/fast_mlsirm/validation.py @@ -101,3 +101,89 @@ def validate_judge( passed=bool(res["pass"]), failed_gates=[g["name"] for g in gates if not g["pass"]], ) + + +@dataclass +class FleissKappaResult: + """Result of :func:`fleiss_kappa`. In exact (Conger) mode ``z`` and + ``p_value`` are NaN and the category arrays are empty, mirroring irr's + ``kappam.fleiss(exact=TRUE)`` which returns neither.""" + + kappa: float + subjects_used: int + z: float + p_value: float + category_kappa: np.ndarray + category_z: np.ndarray + category_p: np.ndarray + + +def fleiss_kappa( + ratings: np.ndarray, + k: int | None = None, + exact: bool = False, +) -> FleissKappaResult: + """Fleiss' kappa for nominal agreement among multiple raters, with the + exact (Conger) chance-agreement variant. + + Reimplements ``kappam.fleiss()`` from CRAN irr 0.85 (R source READ in + full; algorithm source of truth). Model origins — cited as origins only, + NOT READ: Fleiss, J. L. (1971). Measuring nominal scale agreement among + many raters. *Psychological Bulletin, 76*(5), 378-382; Conger, A. J. + (1980). Integration and generalization of kappas for multiple raters. + *Psychological Bulletin, 88*(2), 322-328. Computation runs in the Rust + core (``mlsirm_core::agreement::fleiss_kappa``). + + ``ratings`` is a 2-D ``(n_subjects, n_raters)`` array of integer category + codes ``0..k-1``. NaN or any negative value marks a missing rating and + drops the whole subject row (listwise, as in irr). ``k=None`` infers + ``max(code)+1``; pass ``k`` explicitly to include trailing empty + categories in the category-wise detail (their kappas are NaN, matching + R's 0/0). + """ + from . import _core # computation lives in the Rust core + + if isinstance(ratings, np.ma.MaskedArray): + raise ValueError("masked arrays are not supported; use NaN for missing") + arr = np.asarray(ratings) + if arr.ndim != 2: + raise ValueError("ratings must be a 2-D (subjects x raters) array") + if np.iscomplexobj(arr): + raise ValueError("ratings must be real-valued") + if arr.dtype == object: + for v in arr.flat: + if v is None or isinstance(v, (bool, np.bool_, str, bytes)): + raise ValueError("ratings must be numeric, not boolean/str/None") + arr = arr.astype(np.float64) + if arr.dtype.kind == "b": + raise ValueError("ratings must be integer codes, not booleans") + if arr.dtype.kind not in "fiu": + raise ValueError(f"ratings dtype {arr.dtype} is not numeric") + ns, nr = arr.shape + if arr.dtype.kind == "f": + finite = np.isfinite(arr) + if np.any(np.isinf(arr)): + raise ValueError("ratings must not contain infinities") + if np.any(arr[finite] != np.floor(arr[finite])): + raise ValueError("ratings must be integer category codes") + if np.any(np.abs(arr[finite]) > 2.0**53): + raise ValueError("ratings exceed exact float64 integer range") + codes = np.where(finite, arr, -1.0).astype(np.int64) + else: + codes = arr.astype(np.int64) + if k is None: + if codes.size == 0 or int(codes.max()) < 0: + raise ValueError("cannot infer k: no observed category codes") + k = int(codes.max()) + 1 + res = _core.fleiss_kappa( + np.ascontiguousarray(codes.reshape(-1)), int(ns), int(nr), int(k), bool(exact) + ) + return FleissKappaResult( + kappa=float(res["kappa"]), + subjects_used=int(res["subjects_used"]), + z=float(res["z"]), + p_value=float(res["p_value"]), + category_kappa=np.asarray(res["category_kappa"]), + category_z=np.asarray(res["category_z"]), + category_p=np.asarray(res["category_p"]), + ) diff --git a/tests/test_paper_features.py b/tests/test_paper_features.py index 175b780fe..b8a5fe142 100644 --- a/tests/test_paper_features.py +++ b/tests/test_paper_features.py @@ -10348,3 +10348,78 @@ def test_validation(self): h = 9e307 yh = np.array([[0, 1, 1], [1, 0, h], [1, h, 0]], dtype=float) bratt_mm(yh, t, max_iter=3, tol=2.0) + + +class TestFleiss: + """fleiss_kappa binding + wrapper (irr 0.85 kappam.fleiss oracle anchors).""" + + @staticmethod + def _fk(): + import numpy as np + + return np.array( + [[0, 0, 0, 1], [0, 1, 1, 1], [2, 2, 2, 2], [0, 0, 2, 2], [1, 1, 1, 0]] + ) + + def test_classic_anchor(self): + import numpy as np + + from fast_mlsirm import fleiss_kappa + + res = fleiss_kappa(self._fk()) + assert abs(res.kappa - 139.0 / 399.0) < 1e-14 + assert res.subjects_used == 5 + assert abs(res.z - 2.694739854085488) < 1e-10 + assert abs(res.p_value - 0.007044360582468963) < 5e-7 + np.testing.assert_allclose( + res.category_kappa, [1 / 21, 31 / 91, 43 / 63], rtol=1e-14 + ) + + def test_exact_and_missing(self): + import numpy as np + + from fast_mlsirm import fleiss_kappa + + res = fleiss_kappa(self._fk(), exact=True) + assert abs(res.kappa - 37.0 / 102.0) < 1e-14 + assert np.isnan(res.z) and np.isnan(res.p_value) + assert res.category_kappa.size == 0 + # NaN row drops listwise and reproduces the classic anchor. + withnan = np.vstack([[0.0, np.nan, 1.0, 2.0], self._fk().astype(float)]) + res2 = fleiss_kappa(withnan) + assert res2.subjects_used == 5 + assert abs(res2.kappa - 139.0 / 399.0) < 1e-14 + # Explicit k adds an empty category with NaN detail. + res3 = fleiss_kappa(self._fk(), k=4) + assert np.isnan(res3.category_kappa[3]) + assert abs(res3.kappa - 139.0 / 399.0) < 1e-14 + + def test_validation(self): + import numpy as np + import pytest + + from fast_mlsirm import fleiss_kappa + + fk = self._fk() + with pytest.raises(ValueError, match="2-D"): + fleiss_kappa(fk.reshape(-1)) + with pytest.raises(ValueError, match="masked"): + fleiss_kappa(np.ma.masked_array(fk, mask=False)) + with pytest.raises(ValueError, match="complex|real"): + fleiss_kappa(fk.astype(complex)) + with pytest.raises(ValueError, match="boolean"): + fleiss_kappa(np.array([[True, False], [False, True]], dtype=object)) + with pytest.raises(ValueError, match="boolean"): + fleiss_kappa(np.array([[True, False], [False, True]])) + with pytest.raises(ValueError, match="integer"): + fleiss_kappa(fk + 0.5) + with pytest.raises(ValueError, match="infinit"): + fleiss_kappa(np.array([[np.inf, 0.0], [1.0, 0.0]])) + with pytest.raises(ValueError, match="raters"): + fleiss_kappa(fk[:, :1]) + with pytest.raises(ValueError, match="infer"): + fleiss_kappa(np.full((2, 2), np.nan)) + with pytest.raises(ValueError, match="dropped"): + fleiss_kappa(np.full((2, 2), np.nan), k=2) + with pytest.raises(ValueError, match="degenerate"): + fleiss_kappa(np.ones((3, 2), dtype=int), k=2) diff --git a/tests/unit/agreement_tests.rs b/tests/unit/agreement_tests.rs index 0b3ae850f..43dd61743 100644 --- a/tests/unit/agreement_tests.rs +++ b/tests/unit/agreement_tests.rs @@ -134,3 +134,238 @@ fn rejects_degenerate_inputs() { .iter() .any(|gate| gate.name == "subgroup_smd")); } + +// --------------------------------------------------------------------------- +// fleiss_kappa (irr 0.85 kappam.fleiss; oracle anchors FK1-FK5 with exact +// Fractions). Every assert below reads crate outputs (FleissKappaResult +// fields returned by fleiss_kappa); pins are hand-derived rationals. +// z pins are pure arithmetic (rel 1e-12); p-value pins go through the +// crate's Numerical-Recipes erfc (|error| < 1.2e-7), hence abs 5e-7. +// --------------------------------------------------------------------------- + +/// FK fixture: ns=5, nr=4, k=3 (asymmetric; 5x4 also breaks stride +/// transposition and row-vs-column chanceP confusion). +fn fk_fixture() -> Vec { + vec![ + 0, 0, 0, 1, // S1 + 0, 1, 1, 1, // S2 + 2, 2, 2, 2, // S3 + 0, 0, 2, 2, // S4 + 1, 1, 1, 0, // S5 + ] +} + +fn rel_eq(a: f64, b: f64, tol: f64) -> bool { + (a - b).abs() <= tol * b.abs().max(1.0) +} + +#[test] +fn fk_anchor_fk1_classic() { + // Kills MU1 (agreeP drops -nr centering: kappa would be 113/133), + // MU2 (row-sum chanceP: 1/5 -> kappa 11/24), and via the z pin MU4 + // (variance second-term sign: Sum p q (q-p) = 441/2000 != 0). + let res = fleiss_kappa(&fk_fixture(), 5, 4, 3, false).unwrap(); + assert!( + rel_eq(res.kappa, 139.0 / 399.0, 1e-15), + "kappa {}", + res.kappa + ); + assert_eq!(res.subjects_used, 5); + // var = 181/10830; z = kappa/sqrt(var) (oracle 2.694739854085488). + assert!( + rel_eq(res.z, (139.0 / 399.0) / (181.0_f64 / 10830.0).sqrt(), 1e-12), + "z {}", + res.z + ); + assert!( + (res.p_value - 0.007044360582468963).abs() < 5e-7, + "p {}", + res.p_value + ); +} + +#[test] +fn fk_exact_fk2() { + // Kills MU3 (exact == classic): Conger correction Sum s2_j = 3/50 gives + // chanceP 8/25 and kappa 37/102 != 139/399. + let res = fleiss_kappa(&fk_fixture(), 5, 4, 3, true).unwrap(); + assert!( + rel_eq(res.kappa, 37.0 / 102.0, 1e-15), + "exact kappa {}", + res.kappa + ); + assert_eq!(res.subjects_used, 5); + // Exact mode returns no test statistic or detail (irr returns neither). + assert!(res.z.is_nan() && res.p_value.is_nan()); + assert!(res.category_kappa.is_empty()); + assert!(res.category_z.is_empty()); + assert!(res.category_p.is_empty()); +} + +#[test] +fn fk_missing_drop_fk3() { + // Kills MU5 (missing code counted as a category instead of listwise row + // drop): the prepended row must vanish, leaving FK1 exactly. + let mut ratings = vec![0, -1, 1, 2]; + ratings.extend(fk_fixture()); + let res = fleiss_kappa(&ratings, 6, 4, 3, false).unwrap(); + assert_eq!(res.subjects_used, 5); + let base = fleiss_kappa(&fk_fixture(), 5, 4, 3, false).unwrap(); + assert_eq!(res.kappa, base.kappa, "drop must reproduce FK1 bitwise"); + assert!(rel_eq(res.kappa, 139.0 / 399.0, 1e-15)); +} + +#[test] +fn fk_category_detail_fk4() { + // Kills MU6 (pjk drops the m*nr*p_j subtraction: category kappas would + // be [51/91, 233/273, 73/63]). + let res = fleiss_kappa(&fk_fixture(), 5, 4, 3, false).unwrap(); + let expect = [1.0 / 21.0, 31.0 / 91.0, 43.0 / 63.0]; + assert_eq!(res.category_kappa.len(), 3); + for j in 0..3 { + assert!( + rel_eq(res.category_kappa[j], expect[j], 1e-15), + "kappa_{j} {}", + res.category_kappa[j] + ); + // var_j = 1/30 for all j -> z_j = kappa_j * sqrt(30). + assert!( + rel_eq(res.category_z[j], expect[j] * 30.0_f64.sqrt(), 1e-12), + "z_{j} {}", + res.category_z[j] + ); + } + // Oracle p-values (through math.erfc; crate erfc abs 5e-7). + let expect_p = [ + 0.7942311156261253, + 0.062059828219847915, + 0.00018517759699332675, + ]; + for j in 0..3 { + assert!( + (res.category_p[j] - expect_p[j]).abs() < 5e-7, + "p_{j} {}", + res.category_p[j] + ); + } +} + +#[test] +fn fk_empty_category_fk5() { + // k=4 with no code 3: R's 0/0 -> NaN preserved; overall kappa unchanged + // (C_3 = 0 contributes nothing to chanceP). + let res = fleiss_kappa(&fk_fixture(), 5, 4, 4, false).unwrap(); + assert!(rel_eq(res.kappa, 139.0 / 399.0, 1e-15)); + assert_eq!(res.category_kappa.len(), 4); + assert!(res.category_kappa[3].is_nan()); + assert!(res.category_z[3].is_nan()); + assert!(res.category_p[3].is_nan()); + // Non-empty categories unaffected by the extra level. + assert!(rel_eq(res.category_kappa[0], 1.0 / 21.0, 1e-15)); +} + +#[test] +fn fk_error_contract() { + let fk = fk_fixture(); + assert!(fleiss_kappa(&fk, 0, 4, 3, false).is_err(), "ns == 0"); + assert!(fleiss_kappa(&fk, 5, 1, 3, false).is_err(), "nr < 2"); + assert!(fleiss_kappa(&fk, 5, 4, 1, false).is_err(), "k < 2"); + assert!( + fleiss_kappa(&[], 2_000_000, 2, 2, false).is_err(), + "ns cap before ns*nr" + ); + assert!(fleiss_kappa(&fk, 5, 4, 4, false).is_ok()); + assert!( + fleiss_kappa(&fk[..19], 5, 4, 3, false).is_err(), + "length mismatch" + ); + assert!( + fleiss_kappa(&[0, 1, 3, 0, 1, 2], 3, 2, 3, false).is_err(), + "code >= k" + ); + assert!( + fleiss_kappa(&[-1, 0, 1, -1], 2, 2, 2, false).is_err(), + "all rows dropped" + ); + assert!( + fleiss_kappa(&[1, 1, 1, 1, 1, 1], 3, 2, 2, false).is_err(), + "degenerate chanceP == 1" + ); +} + +/// 500-rep invariance MC: kappa (classic and exact), z, and category detail +/// are invariant under subject-row permutation; kappa and z are invariant +/// under rater-column permutation (ttab is unchanged by either; rtab is +/// permuted across raters, leaving Sum_j s2_j unchanged). Asserts compare +/// two crate outputs, so any asymmetry introduced into the aggregation +/// breaks them. Disclosure: this cannot detect a wrong-but-symmetric +/// formula; the FK1-FK5 value pins above are the discriminating anchors. +#[test] +#[ignore] +fn fk_mc_500_permutation_invariance() { + struct Lcg(u64); + impl Lcg { + fn next_u64(&mut self) -> u64 { + self.0 = self + .0 + .wrapping_mul(6364136223846793005) + .wrapping_add(1442695040888963407); + self.0 + } + fn below(&mut self, n: usize) -> usize { + (self.next_u64() >> 33) as usize % n + } + } + let mut rng = Lcg(0x5eed_f1e5_5000_0001); + let (ns, nr, k) = (8usize, 5usize, 4usize); + for rep in 0..500 { + let ratings: Vec = (0..ns * nr).map(|_| rng.below(k) as i64).collect(); + let Ok(base) = fleiss_kappa(&ratings, ns, nr, k, false) else { + continue; // degenerate draw + }; + let base_ex = fleiss_kappa(&ratings, ns, nr, k, true).unwrap(); + // Fisher-Yates over subject rows. + let mut order: Vec = (0..ns).collect(); + for i in (1..ns).rev() { + order.swap(i, rng.below(i + 1)); + } + let by_subj: Vec = order + .iter() + .flat_map(|&i| ratings[i * nr..(i + 1) * nr].iter().copied()) + .collect(); + let perm = fleiss_kappa(&by_subj, ns, nr, k, false).unwrap(); + let perm_ex = fleiss_kappa(&by_subj, ns, nr, k, true).unwrap(); + assert!( + rel_eq(perm.kappa, base.kappa, 1e-12), + "rep {rep} subject-perm kappa" + ); + assert!(rel_eq(perm.z, base.z, 1e-12), "rep {rep} subject-perm z"); + assert!( + rel_eq(perm_ex.kappa, base_ex.kappa, 1e-12), + "rep {rep} subject-perm exact" + ); + for j in 0..k { + let (a, b) = (perm.category_kappa[j], base.category_kappa[j]); + assert!( + (a.is_nan() && b.is_nan()) || rel_eq(a, b, 1e-12), + "rep {rep} cat {j}" + ); + } + // Fisher-Yates over rater columns. + let mut rorder: Vec = (0..nr).collect(); + for i in (1..nr).rev() { + rorder.swap(i, rng.below(i + 1)); + } + let by_rater: Vec = (0..ns) + .flat_map(|i| { + let ratings = &ratings; + rorder.iter().map(move |&r| ratings[i * nr + r]) + }) + .collect(); + let rperm = fleiss_kappa(&by_rater, ns, nr, k, true).unwrap(); + assert!( + rel_eq(rperm.kappa, base_ex.kappa, 1e-12), + "rep {rep} rater-perm exact kappa" + ); + } +} From 11bc1ed782f8204865497acd98cdd84c49cf902f Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Mon, 27 Jul 2026 10:05:42 +0900 Subject: [PATCH 2/9] fix(fleiss): reject uint64 overflow and lossy k coercion in wrapper Impl-review round 1 findings: - MAJOR: uint64 values above i64::MAX wrapped negative via astype(int64) and were silently dropped as missing; now rejected before conversion. - MINOR: explicit k accepted lossy coercions (3.9, '3', bool); now requires a true integer (int or np.integer, bool excluded). Regression tests added to TestFleiss.test_validation. Co-authored-by: Copilot App <223556219+Copilot@users.noreply.github.com> --- python/fast_mlsirm/validation.py | 6 ++++++ tests/test_paper_features.py | 12 ++++++++++++ 2 files changed, 18 insertions(+) diff --git a/python/fast_mlsirm/validation.py b/python/fast_mlsirm/validation.py index 60dd445c1..99d108c7f 100644 --- a/python/fast_mlsirm/validation.py +++ b/python/fast_mlsirm/validation.py @@ -170,7 +170,13 @@ def fleiss_kappa( raise ValueError("ratings exceed exact float64 integer range") codes = np.where(finite, arr, -1.0).astype(np.int64) else: + if arr.dtype.kind == "u" and arr.size and int(arr.max()) > np.iinfo(np.int64).max: + raise ValueError("ratings values must fit in int64") codes = arr.astype(np.int64) + if k is not None: + if isinstance(k, (bool, np.bool_)) or not isinstance(k, (int, np.integer)): + raise ValueError("k must be an integer") + k = int(k) if k is None: if codes.size == 0 or int(codes.max()) < 0: raise ValueError("cannot infer k: no observed category codes") diff --git a/tests/test_paper_features.py b/tests/test_paper_features.py index b8a5fe142..dd77ec801 100644 --- a/tests/test_paper_features.py +++ b/tests/test_paper_features.py @@ -10423,3 +10423,15 @@ def test_validation(self): fleiss_kappa(np.full((2, 2), np.nan), k=2) with pytest.raises(ValueError, match="degenerate"): fleiss_kappa(np.ones((3, 2), dtype=int), k=2) + # uint64 above i64::MAX must be rejected, not wrapped negative + # (silent listwise drop). Reads the wrapper's guard, killed by + # removing the unsigned-range check before astype(int64). + big = np.array([[2**63, 0], [0, 1], [1, 1]], dtype=np.uint64) + with pytest.raises(ValueError, match="int64"): + fleiss_kappa(big, k=2) + # Explicit k must be a true integer, not lossy-coerced. + for bad_k in (3.9, "3", np.float64(3.0), True): + with pytest.raises(ValueError, match="k must be an integer"): + fleiss_kappa(fk, k=bad_k) + # np.integer k still accepted. + assert fleiss_kappa(fk, k=np.int64(3)).subjects_used == 5 From 675a2a3d33f1d834cc06a41f863baf8ff914560a Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Mon, 27 Jul 2026 10:42:08 +0900 Subject: [PATCH 3/9] Add icc intraclass correlation coefficients (Shrout-Fleiss taxonomy) Transcribed from CRAN irr 0.85 R/icc.R (READ, normative source; Shrout & Fleiss 1979, McGraw & Wong 1996, Bartko 1966 NOT READ, cited as origins only). Rust core computes all six variants (oneway/twoway x consistency/agreement x single/average) from one-pass ANOVA mean squares, the F test of H0: icc = r0 (two-way agreement via Satterthwaite df, preserving the R quirk that both units' CI bounds reuse the nr-scaled plug-in form, icc.R lines 139-141), and unclamped confidence bounds. Listwise NaN row drop; Inf rejected; degenerate zero-variance / icc=1 pivots error instead of leaking non-finite output. Evidence: exact-Fraction oracle on Shrout-Fleiss Table 2 (all six coefficients + scipy F CI pins), 6 EXECUTED mutation kills (MSw divisor, quantile df order, agreement denominator, r0-in-F, dimension map, CI plug-in), MC-500 permutation invariance + Spearman-Brown single/average bridge for all three families. cargo 853 pass; pytest paper suite 340 pass. Co-authored-by: Copilot App <223556219+Copilot@users.noreply.github.com> --- CHANGELOG.md | 1 + crates/fast-mlsirm-py/src/lib.rs | 59 ++++- crates/mlsirm-core/src/reliability.rs | 246 +++++++++++++++++++ python/fast_mlsirm/__init__.py | 4 + python/fast_mlsirm/reliability.py | 109 +++++++++ tests/test_paper_features.py | 87 +++++++ tests/unit/reliability_tests.rs | 338 ++++++++++++++++++++++++++ 7 files changed, 838 insertions(+), 6 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index acc6f377a..742ebe5ee 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -119,6 +119,7 @@ ### Added +- `icc` intraclass correlation coefficients for inter-rater reliability, the full Shrout-Fleiss taxonomy (CRAN irr 0.85 `icc()` `R/icc.R` — READ and normative; Shrout & Fleiss 1979, McGraw & Wong 1996, Bartko 1966 NOT READ, cited as model origins only): `model` oneway/twoway × `type` consistency/agreement × `unit` single/average from one-pass ANOVA mean squares (MSr, MSw, MSc, MSe with sample-variance divisor n−1), the F test of H0: icc = r0 (two-way agreement uses the Satterthwaite df with the R quirk that both units' confidence bounds plug the estimate into the nr-scaled a,b form — icc.R lines 139-141 preserved verbatim), and unclamped confidence bounds. Rows with NaN are dropped listwise (R `na.omit`); infinities rejected; degenerate zero-variance and icc=1 pivots return explicit errors instead of leaking non-finite output. Exact-Fraction oracle anchors on Shrout-Fleiss Table 2 (all six coefficients: 448/2703, 920/1287, 184/635, 1792/4047, 3680/4047, 736/1187 plus scipy F-distribution CI pins), a 6-mutant EXECUTED kill map (MSw divisor, quantile df order, agreement denominator, r0 in F, dimension map, CI plug-in), and an MC-500 test of subject/rater permutation invariance plus the Spearman-Brown single↔average bridge for all three families. - `fleiss_kappa` Fleiss' multi-rater kappa for nominal agreement with the exact (Conger) chance-agreement variant (CRAN irr 0.85 `kappam.fleiss()` `R/kappam.fleiss.R` — READ and normative; Fleiss 1971 and Conger 1980 NOT READ, cited as model origins only): classification-table agreement `agreeP = (1/m)Σᵢ(Σⱼtᵢⱼ² − nr)/(nr(nr−1))`, classic chance `Σⱼpⱼ²` vs exact `Σⱼpⱼ² − (1/nr)Σⱼs²ⱼ` (sample variance over per-rater category proportions; algebra verified against R's `apply(rtab,2,var)` form), Fleiss' large-sample z test and category-wise kappas (classic mode; NaN for empty categories, matching R's 0/0), and listwise row drop for missing ratings. API deviations documented: index codes 0..k-1 with explicit/inferred k, negative-or-NaN = missing, error on degenerate `1 − chanceP = 0` (R returns NaN). Exact-Fraction oracle anchors FK1–FK5 (classic κ = 139/399, exact κ = 37/102, category κ = [1/21, 31/91, 43/63]), a 6-mutant EXECUTED kill map (row-vs-column chance sums, missing-as-category, variance-sign, pjk-centering), and an MC-500 subject/rater permutation-invariance test. - `bratt_mm` Bradley-Terry model with ties fitted by MM (VGAM 1.1-14 `bratt()` family, `R/family.categorical.R` — READ and normative; Bradley & Terry 1952 NOT READ, cited as model origin): `P(i>j) = αᵢ/(αᵢ+αⱼ+α₀)`, `P(tie) = α₀/(αᵢ+αⱼ+α₀)` with a hand-derived supporting-hyperplane MM ascent (same pattern as the crate's `bradley_terry_mm`) and a joint reference rescale of α AND α₀ (likelihood-preserving; verified identity `Σ wins + T = Σ n_ij`). This is the additive-α₀ ties model, NOT Rao-Kupper/Davidson (neither read; disambiguation only). Contract: fractional weighted counts accepted, symmetric ties matrix required, tie-free data rejected (use `bradley_terry_mm`; an API contract, not VGAM behavior), zero-win contestants rejected, n capped at 10000 (O(n²) guard). Exact-Fraction oracle anchors B1–B4 (iteration-1 pins `[1, 27/40, 3/4]`, α₀ = 9/14), a 5-mutant EXECUTED kill map (incl. a tol-separated convergence anchor killing an α₀-blind convergence check), and an MC-500 log-likelihood dominance test. - `predict_rating` / `predict_rating_multi` game-outcome prediction from fitted ratings (CRAN PlayerRatings 1.1-0 `predict.rating` `R/ratings.R` lines 1056–1133 — READ and normative; no journal paper exists for this dispatch, CRAN package provenance only): Elo logistic branch, deviation-shrunk Glicko/Glicko-2/Stephenson branch (`qip3 = 3(ln10/400π)²`, joint shrink over BOTH players' squared deviations), and multi-player EloM branch (`(rating − rowmean)/40`, na.rm row means, optional min-tie placing ranks with NaN kept). R semantics preserved: strict `games < tng` unrated cutoff, `trat` replacement of ALL missing extracted values (unmatched, low-games, stored-NA), `pred >= thresh` binarization with NaN propagation. REDUCED-SCOPE vs R: index-based (−1 = unmatched; caller does name matching), per-game/scalar gamma only. Exact-oracle fixtures P1–P9 and a 7-mutant EXECUTED kill map (incl. both branches' tng comparisons). diff --git a/crates/fast-mlsirm-py/src/lib.rs b/crates/fast-mlsirm-py/src/lib.rs index efd3f0c6e..7981bea26 100644 --- a/crates/fast-mlsirm-py/src/lib.rs +++ b/crates/fast-mlsirm-py/src/lib.rs @@ -109,7 +109,7 @@ use mlsirm_core::rasch_cml::{ use mlsirm_core::reliability::guttman_lambdas as core_guttman_lambdas; use mlsirm_core::reliability::tenberge_mu as core_tenberge_mu; use mlsirm_core::reliability::{ - cronbach_alpha as core_cronbach_alpha, feldt_alpha_ci as core_feldt_alpha_ci, + cronbach_alpha as core_cronbach_alpha, feldt_alpha_ci as core_feldt_alpha_ci, icc as core_icc, separation_reliability as core_separation_reliability, }; use mlsirm_core::rsm::fit_rsm as core_fit_rsm; @@ -3578,6 +3578,47 @@ fn feldt_alpha_ci( Ok(out.into()) } +/// Intraclass correlation coefficients (Shrout & Fleiss, 1979 taxonomy), +/// transcribed from CRAN irr 0.85 `icc.R` (READ; `mlsirm_core::reliability`). +/// `ratings` is row-major ns x nr; rows with NaN are dropped listwise. +/// Returns a dict with `value`, `subjects`, `raters`, `fvalue`, `df1`, +/// `df2`, `p_value`, `lbound`, `ubound`. +#[pyfunction] +fn icc( + py: Python<'_>, + ratings: PyReadonlyArray1<'_, f64>, + ns: usize, + nr: usize, + model: &str, + r#type: &str, + unit: &str, + r0: f64, + conf_level: f64, +) -> PyResult> { + let res = core_icc( + ratings.as_slice()?, + ns, + nr, + model, + r#type, + unit, + r0, + conf_level, + ) + .map_err(PyValueError::new_err)?; + let out = pyo3::types::PyDict::new(py); + out.set_item("value", res.value)?; + out.set_item("subjects", res.subjects)?; + out.set_item("raters", res.raters)?; + out.set_item("fvalue", res.fvalue)?; + out.set_item("df1", res.df1)?; + out.set_item("df2", res.df2)?; + out.set_item("p_value", res.p_value)?; + out.set_item("lbound", res.lbound)?; + out.set_item("ubound", res.ubound)?; + Ok(out.into()) +} + /// Person separation reliability `(SSD - MSE) / SSD` /// (`mlsirm_core::reliability`; transcribed from CRAN eRm `SepRel.R`). /// Returns a dict with `sep_rel`, `ssd`, `mse`, `sep_index`. @@ -5496,13 +5537,17 @@ fn glicko2_rating( .iter() .map(|&v| usize::try_from(v)) .collect::>() - .map_err(|_| PyValueError::new_err("glicko2_rating: player index exceeds platform usize"))?; + .map_err(|_| { + PyValueError::new_err("glicko2_rating: player index exceeds platform usize") + })?; let black: Vec = black .as_slice()? .iter() .map(|&v| usize::try_from(v)) .collect::>() - .map_err(|_| PyValueError::new_err("glicko2_rating: player index exceeds platform usize"))?; + .map_err(|_| { + PyValueError::new_err("glicko2_rating: player index exceeds platform usize") + })?; let res = mlsirm_core::scaling::glicko2_rating( periods.as_slice()?, &white, @@ -5656,7 +5701,6 @@ fn elom_rating( Ok(d.into()) } - /// Prediction-quality metrics for binary-outcome forecasts, PlayerRatings /// `metrics()` semantics (see `mlsirm_core::scaling::metrics_rating`). /// `pred` is flattened row-major nr x np; the return value is the @@ -5811,7 +5855,6 @@ fn predict_rating_multi( Ok(PyArray1::from_slice(py, &out).into()) } - /// Bradley-Terry model with ties (additive alpha0, VGAM `bratt`) fitted /// by MM (see `mlsirm_core::scaling::bratt_mm`). `wins` and `ties` are /// flat row-major n*n matrices (ties symmetric); returns dict with alpha @@ -5871,7 +5914,10 @@ fn fleiss_kappa( d.set_item("subjects_used", res.subjects_used as u64)?; d.set_item("z", res.z)?; d.set_item("p_value", res.p_value)?; - d.set_item("category_kappa", PyArray1::from_slice(py, &res.category_kappa))?; + d.set_item( + "category_kappa", + PyArray1::from_slice(py, &res.category_kappa), + )?; d.set_item("category_z", PyArray1::from_slice(py, &res.category_z))?; d.set_item("category_p", PyArray1::from_slice(py, &res.category_p))?; Ok(d.into()) @@ -8177,6 +8223,7 @@ fn fast_mlsirm_core(m: &Bound<'_, PyModule>) -> PyResult<()> { m.add_function(wrap_pyfunction!(tenberge_mu, m)?)?; m.add_function(wrap_pyfunction!(cronbach_alpha, m)?)?; m.add_function(wrap_pyfunction!(feldt_alpha_ci, m)?)?; + m.add_function(wrap_pyfunction!(icc, m)?)?; m.add_function(wrap_pyfunction!(separation_reliability, m)?)?; m.add_function(wrap_pyfunction!(fit_mixture, m)?)?; m.add_function(wrap_pyfunction!(fit_lltm, m)?)?; diff --git a/crates/mlsirm-core/src/reliability.rs b/crates/mlsirm-core/src/reliability.rs index 863fd768b..5a28913df 100644 --- a/crates/mlsirm-core/src/reliability.rs +++ b/crates/mlsirm-core/src/reliability.rs @@ -731,6 +731,252 @@ pub fn separation_reliability( }) } +/// Intraclass correlation coefficient output (irr `icc`). +#[derive(Debug, Clone, PartialEq)] +pub struct IccResult { + /// The ICC point estimate for the requested variant. + pub value: f64, + /// Complete (post listwise-drop) subject rows used. + pub subjects: u64, + /// Raters (columns). + pub raters: u64, + /// F statistic for H0: icc = r0. + pub fvalue: f64, + /// Numerator degrees of freedom. + pub df1: f64, + /// Denominator degrees of freedom (Satterthwaite, possibly non-integer, + /// for the agreement variants). + pub df2: f64, + /// Upper-tail p-value `P(F_{df1,df2} > fvalue)`. + pub p_value: f64, + /// Lower confidence bound (NOT clamped; can be negative, and below -1 + /// for average-score variants, matching R). + pub lbound: f64, + /// Upper confidence bound (not clamped). + pub ubound: f64, +} + +/// Intraclass correlation coefficients (Shrout-Fleiss family), transcribed +/// from CRAN irr 0.85 `R/icc.R` (READ in full; algorithm source of truth). +/// Model origins — cited as origins only, NOT READ: Shrout, P. E., & +/// Fleiss, J. L. (1979). Intraclass correlations: Uses in assessing rater +/// reliability. *Psychological Bulletin, 86*(2), 420-428; McGraw, K. O., & +/// Wong, S. P. (1996). Forming inferences about some intraclass correlation +/// coefficients. *Psychological Methods, 1*(1), 30-46; Bartko, J. J. +/// (1966). The intraclass correlation coefficient as a measure of +/// reliability. *Psychological Reports, 19*, 3-11. +/// +/// `ratings` is row-major `ns x nr` (subjects x raters); any row containing +/// NaN is dropped listwise before computation (R `na.omit`). `model` is +/// `"oneway"` or `"twoway"`; `typ` is `"consistency"` or `"agreement"` +/// (ignored for oneway, matching R); `unit` is `"single"` or `"average"`. +/// `r0` is the null ICC for the F test; `conf_level` the CI level. +/// +/// ANOVA mean squares (R lines 13-17, `var` = sample variance n-1): +/// `MSr = var(row means)*nr`, `MSw = mean(row variances)`, +/// `MSc = var(col means)*ns`, `MSe = (SStotal - MSr(ns-1) - MSc(nr-1)) / +/// ((ns-1)(nr-1))`. Agreement F tests use the Satterthwaite df with `r0` +/// (R lines 67-75, 128-136); agreement CIs plug the estimated coefficient +/// into the same Satterthwaite form (McGraw & Wong, 1996, as coded at R +/// lines 78-85, 139-146). The average-agreement CI reuses the nr-scaled +/// `a,b` expressions of the single variant verbatim (R lines 139-141) — +/// preserved deliberately. +/// +/// Documented deviations from R: explicit errors (instead of NaN +/// propagation) for fewer than 2 complete rows, `nr < 2`, non-finite +/// (non-NaN) input, out-of-range `r0`/`conf_level`, and degenerate +/// zero/non-finite denominators; dimension caps `ns <= 1e6`, `nr <= 1e4`. +pub fn icc( + ratings: &[f64], + ns: usize, + nr: usize, + model: &str, + typ: &str, + unit: &str, + r0: f64, + conf_level: f64, +) -> Result { + if !matches!(model, "oneway" | "twoway") { + return Err("model must be \"oneway\" or \"twoway\"".into()); + } + if !matches!(typ, "consistency" | "agreement") { + return Err("type must be \"consistency\" or \"agreement\"".into()); + } + if !matches!(unit, "single" | "average") { + return Err("unit must be \"single\" or \"average\"".into()); + } + if !r0.is_finite() || !(0.0..1.0).contains(&r0) { + return Err("r0 must be finite and in [0, 1)".into()); + } + if !conf_level.is_finite() || !(conf_level > 0.0 && conf_level < 1.0) { + return Err("conf_level must be in (0, 1)".into()); + } + if nr < 2 { + return Err("icc needs at least 2 raters".into()); + } + if ns > 1_000_000 || nr > 10_000 { + return Err("icc: dimensions exceed caps (ns <= 1e6, nr <= 1e4)".into()); + } + if ratings.len() != ns * nr { + return Err(format!( + "ratings length {} does not match ns*nr = {}", + ratings.len(), + ns * nr + )); + } + if ratings.iter().any(|v| v.is_infinite()) { + return Err("ratings must not contain infinities (use NaN for missing)".into()); + } + // Listwise drop of rows containing NaN (R na.omit). + let rows: Vec<&[f64]> = (0..ns) + .map(|i| &ratings[i * nr..(i + 1) * nr]) + .filter(|r| r.iter().all(|v| v.is_finite())) + .collect(); + let m = rows.len(); + if m < 2 { + return Err("icc needs at least 2 complete subject rows after dropping missing".into()); + } + let mf = m as f64; + let nrf = nr as f64; + + fn sample_var(xs: &[f64]) -> f64 { + let n = xs.len() as f64; + let mean = xs.iter().sum::() / n; + xs.iter().map(|x| (x - mean).powi(2)).sum::() / (n - 1.0) + } + + let all: Vec = rows.iter().flat_map(|r| r.iter().copied()).collect(); + let ss_total = sample_var(&all) * (mf * nrf - 1.0); + let row_means: Vec = rows.iter().map(|r| r.iter().sum::() / nrf).collect(); + let ms_r = sample_var(&row_means) * nrf; + let ms_w = rows.iter().map(|r| sample_var(r)).sum::() / mf; + let col_means: Vec = (0..nr) + .map(|j| rows.iter().map(|r| r[j]).sum::() / mf) + .collect(); + let ms_c = sample_var(&col_means) * mf; + let ms_e = (ss_total - ms_r * (mf - 1.0) - ms_c * (nrf - 1.0)) / ((mf - 1.0) * (nrf - 1.0)); + if ![ms_r, ms_w, ms_c, ms_e].iter().all(|v| v.is_finite()) { + return Err("icc: ANOVA mean squares are non-finite (inputs too large?)".into()); + } + + let alpha = 1.0 - conf_level; + let q = 1.0 - alpha / 2.0; + let oneway = model == "oneway"; + let consistency = typ == "consistency"; + let single = unit == "single"; + + // Satterthwaite df for the twoway-agreement F test / CI (R lines + // 67-69, 78-80, 128-130, 139-141). `scale` is nr for the single-unit + // a,b and 1 for the average F test; the average CI deliberately reuses + // the nr-scaled form (R quirk, lines 139-141). + let satt = |rho: f64, scale: f64| -> (f64, f64, f64) { + let a = (scale * rho) / (mf * (1.0 - rho)); + let b = 1.0 + (scale * rho * (mf - 1.0)) / (mf * (1.0 - rho)); + let v = (a * ms_c + b * ms_e).powi(2) + / ((a * ms_c).powi(2) / (nrf - 1.0) + (b * ms_e).powi(2) / ((mf - 1.0) * (nrf - 1.0))); + (a, b, v) + }; + + let (value, fvalue, df1, df2, lbound, ubound); + if oneway { + let denom_s = ms_r + (nrf - 1.0) * ms_w; + if single { + value = (ms_r - ms_w) / denom_s; + } else { + value = (ms_r - ms_w) / ms_r; + } + fvalue = if single { + ms_r / ms_w * ((1.0 - r0) / (1.0 + (nrf - 1.0) * r0)) + } else { + ms_r / ms_w * (1.0 - r0) + }; + df1 = mf - 1.0; + df2 = mf * (nrf - 1.0); + let fl = (ms_r / ms_w) / f_quantile(q, df1, df2); + let fu = (ms_r / ms_w) * f_quantile(q, df2, df1); + if single { + lbound = (fl - 1.0) / (fl + nrf - 1.0); + ubound = (fu - 1.0) / (fu + nrf - 1.0); + } else { + lbound = 1.0 - 1.0 / fl; + ubound = 1.0 - 1.0 / fu; + } + } else if consistency { + if single { + value = (ms_r - ms_e) / (ms_r + (nrf - 1.0) * ms_e); + } else { + value = (ms_r - ms_e) / ms_r; + } + fvalue = if single { + ms_r / ms_e * ((1.0 - r0) / (1.0 + (nrf - 1.0) * r0)) + } else { + ms_r / ms_e * (1.0 - r0) + }; + df1 = mf - 1.0; + df2 = (mf - 1.0) * (nrf - 1.0); + let fl = (ms_r / ms_e) / f_quantile(q, df1, df2); + let fu = (ms_r / ms_e) * f_quantile(q, df2, df1); + if single { + lbound = (fl - 1.0) / (fl + nrf - 1.0); + ubound = (fu - 1.0) / (fu + nrf - 1.0); + } else { + lbound = 1.0 - 1.0 / fl; + ubound = 1.0 - 1.0 / fu; + } + } else { + // twoway agreement + if single { + value = (ms_r - ms_e) / (ms_r + (nrf - 1.0) * ms_e + (nrf / mf) * (ms_c - ms_e)); + } else { + value = (ms_r - ms_e) / (ms_r + (ms_c - ms_e) / mf); + } + let (a, b, v) = satt(r0, if single { nrf } else { 1.0 }); + fvalue = ms_r / (a * ms_c + b * ms_e); + df1 = mf - 1.0; + df2 = v; + if !(1.0 - value).is_finite() || (1.0 - value).abs() < 1e-12 { + return Err("icc: degenerate coefficient (icc = 1) — CI undefined".into()); + } + // McGraw & Wong CI: plug the estimate into the nr-scaled a,b (R + // lines 78-80 and, deliberately, 139-141 for the average variant). + let (_a2, _b2, v2) = satt(value, nrf); + if !v2.is_finite() || v2 <= 0.0 { + return Err("icc: degenerate Satterthwaite df in CI".into()); + } + let fl = f_quantile(q, df1, v2); + let fu = f_quantile(q, v2, df1); + if single { + lbound = (mf * (ms_r - fl * ms_e)) + / (fl * (nrf * ms_c + (nrf * mf - nrf - mf) * ms_e) + mf * ms_r); + ubound = (mf * (fu * ms_r - ms_e)) + / (nrf * ms_c + (nrf * mf - nrf - mf) * ms_e + mf * fu * ms_r); + } else { + lbound = (mf * (ms_r - fl * ms_e)) / (fl * (ms_c - ms_e) + mf * ms_r); + ubound = (mf * (fu * ms_r - ms_e)) / (ms_c - ms_e + mf * fu * ms_r); + } + } + if ![value, fvalue, df2, lbound, ubound] + .iter() + .all(|v| v.is_finite()) + { + return Err( + "icc: degenerate ratings (zero-variance denominator produced non-finite output)".into(), + ); + } + let p_value = 1.0 - f_cdf(fvalue, df1, df2); + Ok(IccResult { + value, + subjects: m as u64, + raters: nr as u64, + fvalue, + df1, + df2, + p_value, + lbound, + ubound, + }) +} + #[cfg(test)] #[path = "../../../tests/unit/reliability_tests.rs"] mod tests; diff --git a/python/fast_mlsirm/__init__.py b/python/fast_mlsirm/__init__.py index 2c7a80105..ee58b407a 100644 --- a/python/fast_mlsirm/__init__.py +++ b/python/fast_mlsirm/__init__.py @@ -123,6 +123,8 @@ AlphaCiResult as AlphaCiResult, separation_reliability as separation_reliability, SeparationReliabilityResult as SeparationReliabilityResult, + icc as icc, + IccResult as IccResult, ) from .factor import ( minres_fa as minres_fa, @@ -362,6 +364,8 @@ "AlphaCiResult", "separation_reliability", "SeparationReliabilityResult", + "icc", + "IccResult", "gtheory_pi", "gtheory_pio", "phi_lambda", diff --git a/python/fast_mlsirm/reliability.py b/python/fast_mlsirm/reliability.py index fbab28187..1b8094292 100644 --- a/python/fast_mlsirm/reliability.py +++ b/python/fast_mlsirm/reliability.py @@ -301,3 +301,112 @@ def separation_reliability( mse=float(res["mse"]), sep_index=float(res["sep_index"]), ) + + +@dataclass +class IccResult: + """Intraclass correlation coefficient (irr ``icc``). + + ``value`` is the ICC estimate; ``fvalue``/``df1``/``df2``/``p_value`` + test H0: icc = ``r0`` (upper tail); ``lbound``/``ubound`` are the + two-sided ``conf_level`` interval (unclamped; can drop below -1 for + the one-way average variant). ``subjects`` counts the complete rows + actually used after listwise NaN deletion.""" + + value: float + subjects: int + raters: int + fvalue: float + df1: float + df2: float + p_value: float + lbound: float + ubound: float + + +def icc( + ratings, + model: str = "oneway", + type: str = "consistency", + unit: str = "single", + r0: float = 0.0, + conf_level: float = 0.95, +) -> IccResult: + """Intraclass correlation coefficients for inter-rater reliability + (compute in Rust; transcribed line by line from the CRAN irr 0.85 R + source ``icc.R``, read in full; Shrout & Fleiss, 1979, McGraw & Wong, + 1996, and Bartko, 1966, not read — attribution as cited in Gamer et + al., 2019). Covers the Shrout-Fleiss taxonomy: ``model`` in + {"oneway", "twoway"}, ``type`` in {"consistency", "agreement"}, + ``unit`` in {"single", "average"}. + + ``ratings`` is a 2-D subjects x raters array of continuous scores. + Rows containing NaN are dropped listwise (R ``na.omit``); infinities + are rejected. The two-way agreement F test uses the Satterthwaite + approximation with the null value ``r0``; its confidence bounds plug + the estimate back into the nr-scaled Satterthwaite form for both + units, matching the R source (icc.R lines 139-141) exactly. In + LLM-as-a-Judge quality management this quantifies how consistently + multiple judge models (raters) score the same responses (subjects). + + Verified against an exact-Fraction re-derivation of icc.R executed on + the Shrout-Fleiss Table 2 data (all six variants, plus r0 and + conf_level sweeps); the Spearman-Brown identity between single and + average units was verified for all three families. + + References (APA 7th ed.): + Gamer, M., Lemon, J., Fellows, I., & Singh, P. (2019). *irr: + Various coefficients of interrater reliability and agreement* + (Version 0.84.1; source read at 0.85) [R package]. + https://CRAN.R-project.org/package=irr + Shrout, P. E., & Fleiss, J. L. (1979). Intraclass correlations: + Uses in assessing rater reliability. *Psychological Bulletin, + 86*(2), 420-428. https://doi.org/10.1037/0033-2909.86.2.420 + (as cited in Gamer et al., 2019) + McGraw, K. O., & Wong, S. P. (1996). Forming inferences about + some intraclass correlation coefficients. *Psychological + Methods, 1*(1), 30-46. https://doi.org/10.1037/1082-989X.1.1.30 + (as cited in Gamer et al., 2019) + Bartko, J. J. (1966). The intraclass correlation coefficient as a + measure of reliability. *Psychological Reports, 19*(1), 3-11. + https://doi.org/10.2466/pr0.1966.19.1.3 (as cited in Gamer et + al., 2019) + """ + from .fitstats import _core_module + + core = _core_module() + if core is None or not hasattr(core, "icc"): + raise RuntimeError("icc requires the compiled Rust core") + if isinstance(ratings, np.ma.MaskedArray): + raise ValueError("masked arrays are not supported; use NaN for missing") + arr = np.asarray(ratings) + if arr.dtype == object: + arr = np.asarray(arr, dtype=np.float64) # raises on non-numeric + if np.iscomplexobj(arr): + raise ValueError("ratings must be real-valued") + if arr.dtype.kind == "b": + raise ValueError("ratings must be numeric, not boolean") + if arr.dtype.kind not in "fiu": + raise ValueError("ratings must be a numeric array") + if arr.ndim != 2: + raise ValueError("ratings must be a 2-D subjects x raters array") + x = np.ascontiguousarray(arr, dtype=np.float64) + ns, nr = x.shape + for name, val in (("r0", r0), ("conf_level", conf_level)): + if isinstance(val, bool) or not isinstance(val, (int, float, np.floating, np.integer)): + raise ValueError(f"{name} must be a real number") + res = core.icc( + x.reshape(-1), int(ns), int(nr), str(model), str(type), str(unit), + float(r0), float(conf_level), + ) + return IccResult( + value=float(res["value"]), + subjects=int(res["subjects"]), + raters=int(res["raters"]), + fvalue=float(res["fvalue"]), + df1=float(res["df1"]), + df2=float(res["df2"]), + p_value=float(res["p_value"]), + lbound=float(res["lbound"]), + ubound=float(res["ubound"]), + ) \ No newline at end of file diff --git a/tests/test_paper_features.py b/tests/test_paper_features.py index dd77ec801..9f53f6f12 100644 --- a/tests/test_paper_features.py +++ b/tests/test_paper_features.py @@ -10435,3 +10435,90 @@ def test_validation(self): fleiss_kappa(fk, k=bad_k) # np.integer k still accepted. assert fleiss_kappa(fk, k=np.int64(3)).subjects_used == 5 + + +class TestIcc: + """icc (irr icc.R; Shrout-Fleiss taxonomy). Pins from the exact-Fraction + oracle executed on Shrout-Fleiss (1979) Table 2. Every assert reads + IccResult fields returned by the crate through the wrapper.""" + + SF = [[9, 2, 5, 8], [6, 1, 3, 2], [8, 4, 6, 8], [7, 1, 2, 6], [10, 5, 6, 9], [6, 2, 4, 7]] + + def test_anchor_all_six(self): + import numpy as np + from fast_mlsirm import icc + + want = { + ("oneway", "consistency", "single"): 448 / 2703, + ("twoway", "consistency", "single"): 920 / 1287, + ("twoway", "agreement", "single"): 184 / 635, + ("oneway", "consistency", "average"): 1792 / 4047, + ("twoway", "consistency", "average"): 3680 / 4047, + ("twoway", "agreement", "average"): 736 / 1187, + } + for (m, t, u), v in want.items(): + r = icc(np.array(self.SF, dtype=float), model=m, type=t, unit=u) + assert abs(r.value - v) < 1e-12 + assert r.subjects == 6 and r.raters == 4 + r = icc(np.array(self.SF, dtype=float), model="twoway", type="agreement", unit="single") + assert abs(r.fvalue - 11.02724795640327) < 1e-9 + assert abs(r.df2 - 15.0) < 1e-9 + assert abs(r.p_value - 0.00013456651648433688) < 1e-9 + assert abs(r.lbound - 0.018786513374712013) < 1e-9 + assert abs(r.ubound - 0.7610843696489528) < 1e-9 + + def test_r0_and_conf(self): + import numpy as np + from fast_mlsirm import icc + + x = np.array(self.SF, dtype=float) + r = icc(x, model="twoway", type="agreement", unit="single", r0=0.3) + assert abs(r.fvalue - 0.9561240676364373) < 1e-12 + assert abs(r.df2 - 4.7463353743354775) < 1e-9 + c = icc(x, model="twoway", type="consistency", unit="single", conf_level=0.80) + assert abs(c.lbound - 0.4905340604697688) < 1e-9 + assert abs(c.ubound - 0.8966577783576283) < 1e-9 + + def test_nan_listwise(self): + import numpy as np + from fast_mlsirm import icc + + x = np.array(self.SF, dtype=float) + x[3, 2] = np.nan + r = icc(x, model="twoway", type="consistency", unit="single") + assert r.subjects == 5 + assert abs(r.value - 329 / 459) < 1e-12 + assert abs(r.fvalue - 11.123076923076923) < 1e-12 + + def test_validation(self): + import numpy as np + import pytest + from fast_mlsirm import icc + + x = np.array(self.SF, dtype=float) + with pytest.raises(ValueError): + icc(x, model="threeway") + with pytest.raises(ValueError): + icc(x, type="absolute") + with pytest.raises(ValueError): + icc(x, unit="median") + with pytest.raises(ValueError): + icc(x, r0=1.0) + with pytest.raises(ValueError): + icc(x, conf_level=1.0) + with pytest.raises(ValueError): + icc(np.array([1.0, 2.0, 3.0])) # 1-D + with pytest.raises(ValueError): + icc(np.array([[1.0, np.inf], [2.0, 3.0]])) + with pytest.raises(ValueError): + icc(np.ma.masked_array(x, mask=False)) + with pytest.raises(ValueError): + icc(x.astype(complex)) + with pytest.raises(ValueError): + icc(np.array([[True, False], [False, True]])) + with pytest.raises((ValueError, TypeError)): + icc(np.array([["a", "b"], ["c", "d"]], dtype=object)) + with pytest.raises(ValueError): + icc(x, r0=True) + with pytest.raises(ValueError): + icc(np.full((4, 3), 7.0)) # constant matrix -> degenerate \ No newline at end of file diff --git a/tests/unit/reliability_tests.rs b/tests/unit/reliability_tests.rs index a10eabd47..a55f2e9fe 100644 --- a/tests/unit/reliability_tests.rs +++ b/tests/unit/reliability_tests.rs @@ -643,3 +643,341 @@ fn seprel_overflow_guard() { assert!(separation_reliability(&[1e308, 1e308], &[0.1, 0.1]).is_err()); assert!(separation_reliability(&[0.0, 1.0], &[1e200, 1e200]).is_err()); } + +// ===================== icc (Shrout-Fleiss family) ========================= +// Pins from the exact-Fraction oracle (session files/icc_oracle.py, +// EXECUTED; scipy F dist for qf/pf), transcribed from CRAN irr 0.85 +// R/icc.R. Every assert reads IccResult fields returned by `icc`. +// +// Mutation kills verified by editing the implementation (each FAILED, then +// restored): +// - MU1 MSw normalized by /nr instead of /ns -> ic_anchor_i1 o_s coeff. +// - MU2 FU quantile df order swapped -> ic_anchor_i1 o_s ubound. +// - MU3 agreement-single denominator drops (nr/ns)(MSc-MSe) -> +// ic_anchor_i1 a_s coeff (would collapse to the consistency 920/1287). +// - MU4 r0 ignored in the F statistic -> ic_r0_i3 fvalue pins. +// - MU5 column-major stride bug -> ic_asym_i2 coeff pins. +// - MU6 CI Satterthwaite uses r0 instead of the coefficient plug-in -> +// ic_anchor_i1 a_s lbound/ubound. + +fn icc_rel(a: f64, b: f64, tol: f64) -> bool { + (a - b).abs() <= tol * b.abs().max(1.0) +} + +fn sf_table2() -> Vec { + vec![ + 9.0, 2.0, 5.0, 8.0, // + 6.0, 1.0, 3.0, 2.0, // + 8.0, 4.0, 6.0, 8.0, // + 7.0, 1.0, 2.0, 6.0, // + 10.0, 5.0, 6.0, 9.0, // + 6.0, 2.0, 4.0, 7.0, + ] +} + +fn icc_all6(data: &[f64], ns: usize, nr: usize, r0: f64, conf: f64) -> Vec { + [ + ("oneway", "consistency", "single"), + ("twoway", "consistency", "single"), + ("twoway", "agreement", "single"), + ("oneway", "consistency", "average"), + ("twoway", "consistency", "average"), + ("twoway", "agreement", "average"), + ] + .iter() + .map(|(m, t, u)| icc(data, ns, nr, m, t, u, r0, conf).unwrap()) + .collect() +} + +/// I1: Shrout-Fleiss (1979) Table 2, all six variants (r0=0, conf=.95). +/// Kills MU1 (o_s coeff), MU2 (o_s ubound), MU3 (a_s coeff), MU6 (a_s CI). +#[test] +fn ic_anchor_i1_sf_all_six() { + let r = icc_all6(&sf_table2(), 6, 4, 0.0, 0.95); + let coeffs = [ + 448.0 / 2703.0, + 920.0 / 1287.0, + 184.0 / 635.0, + 1792.0 / 4047.0, + 3680.0 / 4047.0, + 736.0 / 1187.0, + ]; + for (res, want) in r.iter().zip(coeffs) { + assert!(icc_rel(res.value, want, 1e-12), "{} vs {want}", res.value); + assert_eq!(res.subjects, 6); + assert_eq!(res.raters, 4); + assert_eq!(res.df1, 5.0); + } + // F/df2/p: oneway pair shares F, twoway quartet shares F at r0=0 + // (agreement Satterthwaite v == (ns-1)(nr-1) exactly when r0=0). + for i in [0usize, 3] { + assert!(icc_rel(r[i].fvalue, 1.7946784922394678, 1e-12)); + assert_eq!(r[i].df2, 18.0); + assert!((r[i].p_value - 0.16476880834463953).abs() < 1e-9); + } + for i in [1usize, 2, 4, 5] { + assert!(icc_rel(r[i].fvalue, 11.02724795640327, 1e-12)); + assert!(icc_rel(r[i].df2, 15.0, 1e-12)); + assert!((r[i].p_value - 0.00013456651648433688).abs() < 1e-9); + } + let cis = [ + (-0.13293232487475098, 0.722560062328121), + (0.3424647650339252, 0.9458582599553595), + (0.018786513374712013, 0.7610843696489528), + (-0.8844421552381201, 0.9124154203407755), + (0.6756747138163046, 0.9858916781690623), + (0.03944017992139112, 0.9285731833771681), + ]; + for (res, (lb, ub)) in r.iter().zip(cis) { + assert!(icc_rel(res.lbound, lb, 1e-9), "lb {} vs {lb}", res.lbound); + assert!(icc_rel(res.ubound, ub, 1e-9), "ub {} vs {ub}", res.ubound); + } +} + +/// I2: asymmetric 4x3 fixture; kills MU5 (stride/transposition bugs). +#[test] +fn ic_asym_i2() { + let data = vec![1.0, 3.0, 6.0, 2.0, 2.0, 7.0, 4.0, 5.0, 9.0, 3.0, 3.0, 8.0]; + let r = icc_all6(&data, 4, 3, 0.0, 0.95); + let coeffs = [ + -23.0 / 139.0, + 48.0 / 59.0, + 8.0 / 53.0, + -23.0 / 31.0, + 144.0 / 155.0, + 8.0 / 23.0, + ]; + for (res, want) in r.iter().zip(coeffs) { + assert!(icc_rel(res.value, want, 1e-12), "{} vs {want}", res.value); + } + assert!(icc_rel(r[1].fvalue, 14.090909090909092, 1e-12)); + assert_eq!(r[1].df2, 6.0); + assert!((r[1].p_value - 0.00399953361525563).abs() < 1e-9); + assert!(icc_rel(r[2].lbound, -0.005333344546176069, 1e-9)); + assert!(icc_rel(r[2].ubound, 0.7496392441094206, 1e-9)); +} + +/// I3: r0=0.3 changes F/df2/p but not coefficients or CIs. Kills MU4. +#[test] +fn ic_r0_i3() { + let base = icc_all6(&sf_table2(), 6, 4, 0.0, 0.95); + let r = icc_all6(&sf_table2(), 6, 4, 0.3, 0.95); + for (a, b) in r.iter().zip(&base) { + assert_eq!(a.value, b.value); // coeff bitwise-unaffected by r0 + assert_eq!(a.lbound, b.lbound); + assert_eq!(a.ubound, b.ubound); + } + assert!(icc_rel(r[0].fvalue, 0.6611973392461197, 1e-12)); + assert!((r[0].p_value - 0.6573818056947218).abs() < 1e-9); + assert!(icc_rel(r[1].fvalue, 4.062670299727521, 1e-12)); + assert!(icc_rel(r[2].fvalue, 0.9561240676364373, 1e-12)); + assert!(icc_rel(r[2].df2, 4.7463353743354775, 1e-12)); + assert!((r[2].p_value - 0.5219672328433673).abs() < 1e-9); + assert!(icc_rel(r[3].fvalue, 1.2562749445676273, 1e-12)); + assert!(icc_rel(r[4].fvalue, 7.719073569482289, 1e-12)); + assert!(icc_rel(r[5].fvalue, 3.0350332119134347, 1e-12)); + assert!(icc_rel(r[5].df2, 7.136518826153806, 1e-12)); +} + +/// I4: conf=.80 changes CIs but not F/p. Kills conf-ignored mutants. +#[test] +fn ic_conf_i4() { + let base = icc_all6(&sf_table2(), 6, 4, 0.0, 0.95); + let r = icc_all6(&sf_table2(), 6, 4, 0.0, 0.80); + for (a, b) in r.iter().zip(&base) { + assert_eq!(a.fvalue, b.fvalue); + assert_eq!(a.p_value, b.p_value); + } + assert!(icc_rel(r[0].lbound, -0.047857465433838266, 1e-9)); + assert!(icc_rel(r[0].ubound, 0.5441024283192402, 1e-9)); + assert!(icc_rel(r[1].lbound, 0.4905340604697688, 1e-9)); + assert!(icc_rel(r[1].ubound, 0.8966577783576283, 1e-9)); + assert!(icc_rel(r[2].lbound, 0.0788459479007961, 1e-9)); + assert!(icc_rel(r[2].ubound, 0.6022571128642371, 1e-9)); + assert!(icc_rel(r[5].lbound, 0.2287806352847549, 1e-9)); + assert!(icc_rel(r[5].ubound, 0.8597203175540028, 1e-9)); +} + +/// I5: a NaN anywhere in a row drops it listwise; results are bitwise +/// identical to calling with that row removed. +#[test] +fn ic_nan_drop_i5() { + let mut with_nan = sf_table2(); + with_nan[3 * 4 + 2] = f64::NAN; // poison row 3 + let direct: Vec = sf_table2() + .chunks(4) + .enumerate() + .filter(|(i, _)| *i != 3) + .flat_map(|(_, r)| r.to_vec()) + .collect(); + for (m, t, u) in [ + ("oneway", "consistency", "single"), + ("twoway", "consistency", "single"), + ("twoway", "agreement", "average"), + ] { + let a = icc(&with_nan, 6, 4, m, t, u, 0.0, 0.95).unwrap(); + let b = icc(&direct, 5, 4, m, t, u, 0.0, 0.95).unwrap(); + assert_eq!(a.subjects, 5); + assert_eq!(a.value.to_bits(), b.value.to_bits()); + assert_eq!(a.fvalue.to_bits(), b.fvalue.to_bits()); + assert_eq!(a.lbound.to_bits(), b.lbound.to_bits()); + assert_eq!(a.ubound.to_bits(), b.ubound.to_bits()); + } + // Oracle pin for the dropped-row fixture (c_s). + let c = icc(&direct, 5, 4, "twoway", "consistency", "single", 0.0, 0.95).unwrap(); + assert!(icc_rel(c.value, 329.0 / 459.0, 1e-12)); + assert!(icc_rel(c.fvalue, 11.123076923076923, 1e-12)); + assert_eq!(c.df2, 12.0); +} + +/// Every documented error branch. Asserts read crate Err values. +#[test] +fn ic_error_contract() { + let d = sf_table2(); + let ok = |e: Result, frag: &str| { + let msg = e.unwrap_err(); + assert!(msg.contains(frag), "{msg} lacks {frag}"); + }; + ok( + icc(&d, 6, 4, "3way", "consistency", "single", 0.0, 0.95), + "model", + ); + ok(icc(&d, 6, 4, "twoway", "abs", "single", 0.0, 0.95), "type"); + ok( + icc(&d, 6, 4, "twoway", "consistency", "med", 0.0, 0.95), + "unit", + ); + ok( + icc(&d, 6, 4, "twoway", "consistency", "single", 1.0, 0.95), + "r0", + ); + ok( + icc(&d, 6, 4, "twoway", "consistency", "single", -0.1, 0.95), + "r0", + ); + ok( + icc(&d, 6, 4, "twoway", "consistency", "single", f64::NAN, 0.95), + "r0", + ); + ok( + icc(&d, 6, 4, "twoway", "consistency", "single", 0.0, 1.0), + "conf_level", + ); + ok( + icc(&d, 6, 4, "twoway", "consistency", "single", 0.0, f64::NAN), + "conf_level", + ); + ok( + icc(&d[..4], 4, 1, "oneway", "consistency", "single", 0.0, 0.95), + "raters", + ); + ok( + icc(&d, 5, 4, "oneway", "consistency", "single", 0.0, 0.95), + "length", + ); + ok( + icc( + &[f64::INFINITY, 0.0, 1.0, 2.0], + 2, + 2, + "oneway", + "consistency", + "single", + 0.0, + 0.95, + ), + "infinit", + ); + // All-but-one rows dropped -> too few complete rows. + let mut nan_heavy = sf_table2(); + for i in 0..5 { + nan_heavy[i * 4] = f64::NAN; + } + ok( + icc( + &nan_heavy, + 6, + 4, + "oneway", + "consistency", + "single", + 0.0, + 0.95, + ), + "complete subject rows", + ); + // Constant matrix: every mean square is 0 -> degenerate. + let konst = vec![3.0; 12]; + assert!(icc(&konst, 4, 3, "oneway", "consistency", "single", 0.0, 0.95).is_err()); + // Perfect agreement (rows differ, columns identical): MSe == MSw == 0, + // icc == 1 -> degenerate CI/pivot must Err, not leak inf/NaN. + let perfect = vec![1.0, 1.0, 2.0, 2.0, 3.0, 3.0, 4.0, 4.0]; + assert!(icc(&perfect, 4, 2, "twoway", "agreement", "single", 0.0, 0.95).is_err()); + assert!(icc(&perfect, 4, 2, "twoway", "consistency", "single", 0.0, 0.95).is_err()); + // Huge magnitudes overflow the SS accumulation -> explicit error. + let huge = vec![1e300, -1e300, 1e300, -1e300, -1e300, 1e300, 1e300, 1e300]; + assert!(icc(&huge, 4, 2, "oneway", "consistency", "single", 0.0, 0.95).is_err()); +} + +/// MC-500: subject- and rater-permutation invariance plus the +/// Spearman-Brown identity between the crate's own single- and +/// average-unit outputs, for all three model families (algebraic identity +/// verified in spec review; both sides read crate outputs, so a mutation +/// breaking either unit's formula independently fails the bridge). +#[test] +#[ignore] +fn ic_mc_500_invariance() { + let mut state = 0x1CC5EEDu64; + let mut uni = |s: &mut u64| -> f64 { + *s = s + .wrapping_mul(6364136223846793005) + .wrapping_add(1442695040888963407); + (*s >> 11) as f64 / (1u64 << 53) as f64 + }; + for rep in 0..500 { + let ns = 4 + (uni(&mut state) * 6.0) as usize; + let nr = 2 + (uni(&mut state) * 4.0) as usize; + let data: Vec = (0..ns * nr) + .map(|_| (uni(&mut state) * 10.0).round()) + .collect(); + let variants = [ + ("oneway", "consistency"), + ("twoway", "consistency"), + ("twoway", "agreement"), + ]; + // Skip degenerate draws (constant matrices, near-1 CI pivots etc.) + // in ANY variant/unit combination. + if variants.iter().any(|(m, t)| { + icc(&data, ns, nr, m, t, "single", 0.0, 0.95).is_err() + || icc(&data, ns, nr, m, t, "average", 0.0, 0.95).is_err() + }) { + continue; + } + // Subject reversal. + let rev_rows: Vec = (0..ns) + .rev() + .flat_map(|i| data[i * nr..(i + 1) * nr].to_vec()) + .collect(); + // Rater reversal. + let rev_cols: Vec = (0..ns) + .flat_map(|i| (0..nr).rev().map(move |j| (i, j))) + .map(|(i, j)| data[i * nr + j]) + .collect(); + for (m, t) in variants { + let s = icc(&data, ns, nr, m, t, "single", 0.0, 0.95).unwrap(); + let a = icc(&data, ns, nr, m, t, "average", 0.0, 0.95).unwrap(); + let sr = icc(&rev_rows, ns, nr, m, t, "single", 0.0, 0.95).unwrap(); + let sc = icc(&rev_cols, ns, nr, m, t, "single", 0.0, 0.95).unwrap(); + assert!(icc_rel(sr.value, s.value, 1e-10), "rep {rep} row-perm"); + assert!(icc_rel(sc.value, s.value, 1e-10), "rep {rep} col-perm"); + // Spearman-Brown bridge between two crate outputs. + let k = nr as f64; + let sb = k * s.value / (1.0 + (k - 1.0) * s.value); + assert!( + icc_rel(a.value, sb, 1e-10), + "rep {rep} {m}/{t} SB: {} vs {sb}", + a.value + ); + } + } +} From 88c56c025095fb91ace24015761a23da8b85a1b3 Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Mon, 27 Jul 2026 10:53:05 +0900 Subject: [PATCH 4/9] Reject object-dtype boolean ratings in the icc wrapper Impl-review finding (MINOR): object arrays of Python/numpy bools were silently coerced to 0.0/1.0 before the bool-dtype check, bypassing the boolean rejection contract. Scan object arrays for bool elements before the float64 conversion; regression test added. Co-authored-by: Copilot App <223556219+Copilot@users.noreply.github.com> --- python/fast_mlsirm/reliability.py | 2 ++ tests/test_paper_features.py | 2 ++ 2 files changed, 4 insertions(+) diff --git a/python/fast_mlsirm/reliability.py b/python/fast_mlsirm/reliability.py index 1b8094292..e984654cd 100644 --- a/python/fast_mlsirm/reliability.py +++ b/python/fast_mlsirm/reliability.py @@ -381,6 +381,8 @@ def icc( raise ValueError("masked arrays are not supported; use NaN for missing") arr = np.asarray(ratings) if arr.dtype == object: + if any(isinstance(v, (bool, np.bool_)) for v in arr.flat): + raise ValueError("ratings must be numeric, not boolean") arr = np.asarray(arr, dtype=np.float64) # raises on non-numeric if np.iscomplexobj(arr): raise ValueError("ratings must be real-valued") diff --git a/tests/test_paper_features.py b/tests/test_paper_features.py index 9f53f6f12..50f3d65d4 100644 --- a/tests/test_paper_features.py +++ b/tests/test_paper_features.py @@ -10516,6 +10516,8 @@ def test_validation(self): icc(x.astype(complex)) with pytest.raises(ValueError): icc(np.array([[True, False], [False, True]])) + with pytest.raises(ValueError): + icc(np.array([[True, False], [False, True]], dtype=object)) with pytest.raises((ValueError, TypeError)): icc(np.array([["a", "b"], ["c", "d"]], dtype=object)) with pytest.raises(ValueError): From 21a544c67d9f4c6a501bf0b4c8d1b7d46c8e98ab Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Mon, 27 Jul 2026 10:56:41 +0900 Subject: [PATCH 5/9] Reject nested 0-D boolean arrays in object-dtype icc ratings Round-2 impl-review finding (MINOR): an object array whose element is a 0-D np.ndarray of dtype bool coerced to 1.0 and bypassed the boolean rejection. The object-array scan now also rejects ndarray elements with boolean dtype; regression test added. Co-authored-by: Copilot App <223556219+Copilot@users.noreply.github.com> --- python/fast_mlsirm/reliability.py | 6 +++++- tests/test_paper_features.py | 2 ++ 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/python/fast_mlsirm/reliability.py b/python/fast_mlsirm/reliability.py index e984654cd..09117c5d6 100644 --- a/python/fast_mlsirm/reliability.py +++ b/python/fast_mlsirm/reliability.py @@ -381,7 +381,11 @@ def icc( raise ValueError("masked arrays are not supported; use NaN for missing") arr = np.asarray(ratings) if arr.dtype == object: - if any(isinstance(v, (bool, np.bool_)) for v in arr.flat): + if any( + isinstance(v, (bool, np.bool_)) + or (isinstance(v, np.ndarray) and v.dtype.kind == "b") + for v in arr.flat + ): raise ValueError("ratings must be numeric, not boolean") arr = np.asarray(arr, dtype=np.float64) # raises on non-numeric if np.iscomplexobj(arr): diff --git a/tests/test_paper_features.py b/tests/test_paper_features.py index 50f3d65d4..2e6144af7 100644 --- a/tests/test_paper_features.py +++ b/tests/test_paper_features.py @@ -10518,6 +10518,8 @@ def test_validation(self): icc(np.array([[True, False], [False, True]])) with pytest.raises(ValueError): icc(np.array([[True, False], [False, True]], dtype=object)) + with pytest.raises(ValueError): + icc(np.array([[np.array(True), 2.0], [3.0, 4.0]], dtype=object)) with pytest.raises((ValueError, TypeError)): icc(np.array([["a", "b"], ["c", "d"]], dtype=object)) with pytest.raises(ValueError): From c51ad5f1cc21948e4300efe608b795729e0e66c0 Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Mon, 27 Jul 2026 11:01:23 +0900 Subject: [PATCH 6/9] Unwrap nested 0-D ndarrays in the icc object-array boolean scan Round-3 impl-review finding (MINOR): 0-D object-dtype ndarrays wrapping a bool still bypassed the scan. The guard now iteratively unwraps 0-D ndarray elements of any dtype via .item() before the bool check, closing the wrapper-nesting family of bypasses at the root; regression tests for object-wrapped bool and np.bool_ added. 0-D float ndarray elements remain accepted and bitwise-match the plain-float result. Co-authored-by: Copilot App <223556219+Copilot@users.noreply.github.com> --- python/fast_mlsirm/reliability.py | 14 +++++++++----- tests/test_paper_features.py | 8 ++++++++ 2 files changed, 17 insertions(+), 5 deletions(-) diff --git a/python/fast_mlsirm/reliability.py b/python/fast_mlsirm/reliability.py index 09117c5d6..4e8aef8b9 100644 --- a/python/fast_mlsirm/reliability.py +++ b/python/fast_mlsirm/reliability.py @@ -381,11 +381,15 @@ def icc( raise ValueError("masked arrays are not supported; use NaN for missing") arr = np.asarray(ratings) if arr.dtype == object: - if any( - isinstance(v, (bool, np.bool_)) - or (isinstance(v, np.ndarray) and v.dtype.kind == "b") - for v in arr.flat - ): + def _is_bool(v: object) -> bool: + # Unwrap arbitrarily nested 0-D ndarrays (any dtype). + while isinstance(v, np.ndarray) and v.ndim == 0: + v = v.item() + return isinstance(v, (bool, np.bool_)) or ( + isinstance(v, np.ndarray) and v.dtype.kind == "b" + ) + + if any(_is_bool(v) for v in arr.flat): raise ValueError("ratings must be numeric, not boolean") arr = np.asarray(arr, dtype=np.float64) # raises on non-numeric if np.iscomplexobj(arr): diff --git a/tests/test_paper_features.py b/tests/test_paper_features.py index 2e6144af7..2d4ba599e 100644 --- a/tests/test_paper_features.py +++ b/tests/test_paper_features.py @@ -10520,6 +10520,14 @@ def test_validation(self): icc(np.array([[True, False], [False, True]], dtype=object)) with pytest.raises(ValueError): icc(np.array([[np.array(True), 2.0], [3.0, 4.0]], dtype=object)) + with pytest.raises(ValueError): + icc(np.array([[np.array(True, dtype=object), 2.0], [3.0, 4.0]], dtype=object)) + with pytest.raises(ValueError): + icc( + np.array( + [[np.array(np.bool_(True), dtype=object), 2.0], [3.0, 4.0]], dtype=object + ) + ) with pytest.raises((ValueError, TypeError)): icc(np.array([["a", "b"], ["c", "d"]], dtype=object)) with pytest.raises(ValueError): From b9e12f5622cfccd499024e69728e14bff8565014 Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Mon, 27 Jul 2026 11:07:40 +0900 Subject: [PATCH 7/9] fix(reliability): whitelist numeric scalars in icc object-dtype guard Replace the boolean blacklist with a numeric whitelist: after unwrapping 0-D ndarrays, only int/float/np.integer/np.floating scalars are accepted (bool excluded as an int subclass). This closes the np.void structured- scalar bypass and any future exotic-scalar coercion path. Co-authored-by: Copilot App <223556219+Copilot@users.noreply.github.com> --- python/fast_mlsirm/reliability.py | 17 ++++++++++------- tests/test_paper_features.py | 3 +++ 2 files changed, 13 insertions(+), 7 deletions(-) diff --git a/python/fast_mlsirm/reliability.py b/python/fast_mlsirm/reliability.py index 4e8aef8b9..9aac6ccfc 100644 --- a/python/fast_mlsirm/reliability.py +++ b/python/fast_mlsirm/reliability.py @@ -381,16 +381,19 @@ def icc( raise ValueError("masked arrays are not supported; use NaN for missing") arr = np.asarray(ratings) if arr.dtype == object: - def _is_bool(v: object) -> bool: + def _is_numeric(v: object) -> bool: # Unwrap arbitrarily nested 0-D ndarrays (any dtype). while isinstance(v, np.ndarray) and v.ndim == 0: v = v.item() - return isinstance(v, (bool, np.bool_)) or ( - isinstance(v, np.ndarray) and v.dtype.kind == "b" - ) - - if any(_is_bool(v) for v in arr.flat): - raise ValueError("ratings must be numeric, not boolean") + # Whitelist: real int/float scalars only (bool is an int + # subclass and np.bool_/np.void etc. are excluded by not + # being in the whitelist). + if isinstance(v, (bool, np.bool_)): + return False + return isinstance(v, (int, float, np.integer, np.floating)) + + if not all(_is_numeric(v) for v in arr.flat): + raise ValueError("ratings must contain only real numeric scalars") arr = np.asarray(arr, dtype=np.float64) # raises on non-numeric if np.iscomplexobj(arr): raise ValueError("ratings must be real-valued") diff --git a/tests/test_paper_features.py b/tests/test_paper_features.py index 2d4ba599e..82d9e085c 100644 --- a/tests/test_paper_features.py +++ b/tests/test_paper_features.py @@ -10530,6 +10530,9 @@ def test_validation(self): ) with pytest.raises((ValueError, TypeError)): icc(np.array([["a", "b"], ["c", "d"]], dtype=object)) + void_bool = np.zeros((), dtype=[("f", "?")])[()] # np.void w/ bool field + with pytest.raises(ValueError): + icc(np.array([[void_bool, 2.0], [3.0, 4.0]], dtype=object)) with pytest.raises(ValueError): icc(x, r0=True) with pytest.raises(ValueError): From 9ebd18fc72698edc92fc872f9af3ecd30c5bb5c4 Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Mon, 27 Jul 2026 11:20:09 +0900 Subject: [PATCH 8/9] fix(reliability): reject object-dtype ratings outright in icc Five review rounds produced successive bypasses of per-element vetting (object bools, nested 0-D wrappers, np.void, timedelta64, self- referential 0-D arrays hanging the unwrap loop, __float__-lying int/float subclasses). Numeric data never requires object dtype, so the wrapper now rejects it categorically, eliminating the entire coercion attack class. Co-authored-by: Copilot App <223556219+Copilot@users.noreply.github.com> --- python/fast_mlsirm/reliability.py | 21 +++++++-------------- 1 file changed, 7 insertions(+), 14 deletions(-) diff --git a/python/fast_mlsirm/reliability.py b/python/fast_mlsirm/reliability.py index 9aac6ccfc..40b67dfe6 100644 --- a/python/fast_mlsirm/reliability.py +++ b/python/fast_mlsirm/reliability.py @@ -381,20 +381,13 @@ def icc( raise ValueError("masked arrays are not supported; use NaN for missing") arr = np.asarray(ratings) if arr.dtype == object: - def _is_numeric(v: object) -> bool: - # Unwrap arbitrarily nested 0-D ndarrays (any dtype). - while isinstance(v, np.ndarray) and v.ndim == 0: - v = v.item() - # Whitelist: real int/float scalars only (bool is an int - # subclass and np.bool_/np.void etc. are excluded by not - # being in the whitelist). - if isinstance(v, (bool, np.bool_)): - return False - return isinstance(v, (int, float, np.integer, np.floating)) - - if not all(_is_numeric(v) for v in arr.flat): - raise ValueError("ratings must contain only real numeric scalars") - arr = np.asarray(arr, dtype=np.float64) # raises on non-numeric + # Rounds 1-5 of adversarial review showed per-element vetting of + # object arrays is an unwinnable arms race (bool, 0-D wrappers, + # np.void, timedelta64, self-referential arrays, __float__-lying + # subclasses). Numeric input never needs object dtype, so reject it. + raise ValueError( + "object-dtype arrays are not supported; pass a numeric array" + ) if np.iscomplexobj(arr): raise ValueError("ratings must be real-valued") if arr.dtype.kind == "b": From 495e77bdd6e326a34e959eaa47777e77b3c4a689 Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Fri, 31 Jul 2026 21:38:49 +0900 Subject: [PATCH 9/9] Add kripp_alpha Krippendorff's alpha (irr 0.85 kripp.alpha, READ source) (#311) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * Add kripp_alpha Krippendorff's alpha for inter-rater agreement Rust reimplementation of CRAN irr 0.85 kripp.alpha() (R/kripp.alpha.R, READ and normative; Krippendorff 1980 NOT READ, cited as method origin only). Coincidence matrix over unordered rater pairs with the irr divisor quirk preserved verbatim (mc = #nonmissing-1 per column only when any value is missing, else 1), all four metrics (nominal, ordinal half-endpoint weights, interval, ratio), alpha = 1 when fewer than two observed levels. Documented deviations: all-missing, infinities, and ratio zero-sum level pairs are explicit errors. Evidence: exact-Fraction oracle anchors K1-K4 (nominal 113/152, ordinal 108577/133160, interval 951/1120, ratio 18222619/22852465, nmv=40; no-NA quirk pin 43/72 vs m-1 mutant's 11/18), 6-mutant EXECUTED kill map, MC-500 permutation invariance. cargo 858 pass, pytest paper suite 344 pass. Co-authored-by: Copilot App <223556219+Copilot@users.noreply.github.com> * Reject integer rating levels beyond 2**53 in kripp_alpha wrapper Review finding (MAJOR): int64/uint64 rating labels beyond 2**53 are not exactly representable as float64, so distinct levels silently collapsed during the cast — complete disagreement returned alpha=1 with a single level. The wrapper now rejects any integer array with values outside [-2**53, 2**53]; the boundary itself remains accepted (exactly representable). Regression test executes the reviewer's repro and was red-green verified against the unguarded version. Co-authored-by: Copilot App <223556219+Copilot@users.noreply.github.com> --------- Co-authored-by: Copilot App <223556219+Copilot@users.noreply.github.com> --- CHANGELOG.md | 1 + crates/fast-mlsirm-py/src/lib.rs | 27 +++- crates/mlsirm-core/src/reliability.rs | 184 ++++++++++++++++++++++++ python/fast_mlsirm/__init__.py | 4 + python/fast_mlsirm/reliability.py | 84 ++++++++++- tests/test_paper_features.py | 97 ++++++++++++- tests/unit/reliability_tests.rs | 193 ++++++++++++++++++++++++++ 7 files changed, 587 insertions(+), 3 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 742ebe5ee..41022ace8 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -119,6 +119,7 @@ ### Added +- `kripp_alpha` Krippendorff's alpha for inter-rater agreement (CRAN irr 0.85 `kripp.alpha()` `R/kripp.alpha.R` — READ and normative; Krippendorff 1980 NOT READ, cited as method origin only): coincidence matrix over unordered rater pairs per subject with the irr divisor quirk preserved verbatim (`mc = #nonmissing − 1` per column ONLY when the matrix contains any missing value, else 1 — complete-data alpha differs from the m−1 convention), diagonal increment `2/mc`, mirror by assignment, `nmatchval` as total cell mass, and all four distance metrics (nominal, ordinal with half-endpoint coincidence-row-sum weights, interval, ratio) feeding `α = 1 − (nmatchval−1)·Σ(utcm·δ²)/Σ(nc_c·nc_k·δ²)`. Fewer than 2 observed levels yields α = 1 (R line 45). Documented deviations: all-missing matrix, infinities, and ratio level pairs summing to zero are explicit errors (R would return α = 1, propagate, or emit Inf/NaN). Exact-Fraction oracle anchors K1–K4 (irr man-page matrix: nominal 113/152, ordinal 108577/133160, interval 951/1120, ratio 18222619/22852465, nmv = 40; no-NA quirk pin 43/72 vs the m−1 mutant's 11/18), a 6-mutant EXECUTED kill map (diagonal 1/mc, mc always m−1, ordinal full weights, interval |δ|, nmv off-diagonal only, num×nc products), and an MC-500 rater/subject permutation-invariance test. - `icc` intraclass correlation coefficients for inter-rater reliability, the full Shrout-Fleiss taxonomy (CRAN irr 0.85 `icc()` `R/icc.R` — READ and normative; Shrout & Fleiss 1979, McGraw & Wong 1996, Bartko 1966 NOT READ, cited as model origins only): `model` oneway/twoway × `type` consistency/agreement × `unit` single/average from one-pass ANOVA mean squares (MSr, MSw, MSc, MSe with sample-variance divisor n−1), the F test of H0: icc = r0 (two-way agreement uses the Satterthwaite df with the R quirk that both units' confidence bounds plug the estimate into the nr-scaled a,b form — icc.R lines 139-141 preserved verbatim), and unclamped confidence bounds. Rows with NaN are dropped listwise (R `na.omit`); infinities rejected; degenerate zero-variance and icc=1 pivots return explicit errors instead of leaking non-finite output. Exact-Fraction oracle anchors on Shrout-Fleiss Table 2 (all six coefficients: 448/2703, 920/1287, 184/635, 1792/4047, 3680/4047, 736/1187 plus scipy F-distribution CI pins), a 6-mutant EXECUTED kill map (MSw divisor, quantile df order, agreement denominator, r0 in F, dimension map, CI plug-in), and an MC-500 test of subject/rater permutation invariance plus the Spearman-Brown single↔average bridge for all three families. - `fleiss_kappa` Fleiss' multi-rater kappa for nominal agreement with the exact (Conger) chance-agreement variant (CRAN irr 0.85 `kappam.fleiss()` `R/kappam.fleiss.R` — READ and normative; Fleiss 1971 and Conger 1980 NOT READ, cited as model origins only): classification-table agreement `agreeP = (1/m)Σᵢ(Σⱼtᵢⱼ² − nr)/(nr(nr−1))`, classic chance `Σⱼpⱼ²` vs exact `Σⱼpⱼ² − (1/nr)Σⱼs²ⱼ` (sample variance over per-rater category proportions; algebra verified against R's `apply(rtab,2,var)` form), Fleiss' large-sample z test and category-wise kappas (classic mode; NaN for empty categories, matching R's 0/0), and listwise row drop for missing ratings. API deviations documented: index codes 0..k-1 with explicit/inferred k, negative-or-NaN = missing, error on degenerate `1 − chanceP = 0` (R returns NaN). Exact-Fraction oracle anchors FK1–FK5 (classic κ = 139/399, exact κ = 37/102, category κ = [1/21, 31/91, 43/63]), a 6-mutant EXECUTED kill map (row-vs-column chance sums, missing-as-category, variance-sign, pjk-centering), and an MC-500 subject/rater permutation-invariance test. - `bratt_mm` Bradley-Terry model with ties fitted by MM (VGAM 1.1-14 `bratt()` family, `R/family.categorical.R` — READ and normative; Bradley & Terry 1952 NOT READ, cited as model origin): `P(i>j) = αᵢ/(αᵢ+αⱼ+α₀)`, `P(tie) = α₀/(αᵢ+αⱼ+α₀)` with a hand-derived supporting-hyperplane MM ascent (same pattern as the crate's `bradley_terry_mm`) and a joint reference rescale of α AND α₀ (likelihood-preserving; verified identity `Σ wins + T = Σ n_ij`). This is the additive-α₀ ties model, NOT Rao-Kupper/Davidson (neither read; disambiguation only). Contract: fractional weighted counts accepted, symmetric ties matrix required, tie-free data rejected (use `bradley_terry_mm`; an API contract, not VGAM behavior), zero-win contestants rejected, n capped at 10000 (O(n²) guard). Exact-Fraction oracle anchors B1–B4 (iteration-1 pins `[1, 27/40, 3/4]`, α₀ = 9/14), a 5-mutant EXECUTED kill map (incl. a tol-separated convergence anchor killing an α₀-blind convergence check), and an MC-500 log-likelihood dominance test. diff --git a/crates/fast-mlsirm-py/src/lib.rs b/crates/fast-mlsirm-py/src/lib.rs index 7981bea26..d4cf731a6 100644 --- a/crates/fast-mlsirm-py/src/lib.rs +++ b/crates/fast-mlsirm-py/src/lib.rs @@ -110,7 +110,7 @@ use mlsirm_core::reliability::guttman_lambdas as core_guttman_lambdas; use mlsirm_core::reliability::tenberge_mu as core_tenberge_mu; use mlsirm_core::reliability::{ cronbach_alpha as core_cronbach_alpha, feldt_alpha_ci as core_feldt_alpha_ci, icc as core_icc, - separation_reliability as core_separation_reliability, + kripp_alpha as core_kripp_alpha, separation_reliability as core_separation_reliability, }; use mlsirm_core::rsm::fit_rsm as core_fit_rsm; use mlsirm_core::rt::{ @@ -3619,6 +3619,30 @@ fn icc( Ok(out.into()) } +/// Krippendorff's alpha (`mlsirm_core::reliability`; transcribed from +/// CRAN irr 0.85 `kripp.alpha.R`, READ). `ratings` is row-major +/// nraters x nsubjects; NaN marks missing. `method` is one of "nominal", +/// "ordinal", "interval", "ratio". Returns a dict with `value`, +/// `subjects`, `raters`, `levels`, `nmatchval`. +#[pyfunction] +fn kripp_alpha( + py: Python<'_>, + ratings: PyReadonlyArray1<'_, f64>, + nraters: usize, + nsubjects: usize, + method: &str, +) -> PyResult> { + let res = core_kripp_alpha(ratings.as_slice()?, nraters, nsubjects, method) + .map_err(PyValueError::new_err)?; + let out = pyo3::types::PyDict::new(py); + out.set_item("value", res.value)?; + out.set_item("subjects", res.subjects)?; + out.set_item("raters", res.raters)?; + out.set_item("levels", res.levels)?; + out.set_item("nmatchval", res.nmatchval)?; + Ok(out.into()) +} + /// Person separation reliability `(SSD - MSE) / SSD` /// (`mlsirm_core::reliability`; transcribed from CRAN eRm `SepRel.R`). /// Returns a dict with `sep_rel`, `ssd`, `mse`, `sep_index`. @@ -8224,6 +8248,7 @@ fn fast_mlsirm_core(m: &Bound<'_, PyModule>) -> PyResult<()> { m.add_function(wrap_pyfunction!(cronbach_alpha, m)?)?; m.add_function(wrap_pyfunction!(feldt_alpha_ci, m)?)?; m.add_function(wrap_pyfunction!(icc, m)?)?; + m.add_function(wrap_pyfunction!(kripp_alpha, m)?)?; m.add_function(wrap_pyfunction!(separation_reliability, m)?)?; m.add_function(wrap_pyfunction!(fit_mixture, m)?)?; m.add_function(wrap_pyfunction!(fit_lltm, m)?)?; diff --git a/crates/mlsirm-core/src/reliability.rs b/crates/mlsirm-core/src/reliability.rs index 5a28913df..5288d4ec4 100644 --- a/crates/mlsirm-core/src/reliability.rs +++ b/crates/mlsirm-core/src/reliability.rs @@ -977,6 +977,190 @@ pub fn icc( }) } +/// Result of [`kripp_alpha`]. +#[derive(Debug, Clone)] +pub struct KrippResult { + /// Krippendorff's alpha estimate. + pub value: f64, + /// Number of subject columns as given (R reports `dim(x)[2]`). + pub subjects: u64, + /// Number of raters (matrix rows). + pub raters: u64, + /// Number of distinct finite rating levels. + pub levels: u64, + /// Total coincidence-matrix mass (R `nmatchval`). + pub nmatchval: f64, +} + +/// Krippendorff's alpha for a raters x subjects rating matrix. +/// +/// Transcribed from CRAN irr 0.85 `R/kripp.alpha.R` (READ in full, 66 +/// lines; algorithm source of truth). Method origin ? cited as origin +/// only, NOT READ: Krippendorff, K. (1980). *Content analysis: An +/// introduction to its methodology*. Sage. +/// +/// Verified against the R source: +/// - Levels are the sorted unique finite values (R `levels(as.factor(x))`, +/// line 7; base R sorts numeric factor levels in ascending numeric +/// order, not lexicographically). +/// - Coincidence matrix (lines 12-25): for every unordered rater pair in +/// a subject column with both values present, cell `(a, b)` gains +/// `(1 + (a == b)) / mc[col]` and the mirror cell is set by assignment +/// (line 21). `mc[col]` is `#nonmissing - 1` only when the matrix +/// contains at least one missing value anywhere, else `1` for every +/// column (lines 12-13). The no-missing divisor of 1 is a documented +/// irr quirk preserved verbatim ? it is NOT the `m - 1` convention and +/// changes both `nmatchval` and alpha on complete data. +/// - `nmatchval` sums all cells (line 26). Fewer than 2 observed levels +/// yields alpha = 1 (line 45). +/// - Distance metrics (lines 50-59) with `nc` the coincidence row sums: +/// nominal `1`; ordinal `(nc_c/2 + sum_{g=c+1}^{k-1} nc_g + nc_k/2)^2`; +/// interval `(v_c - v_k)^2`; ratio `((v_c - v_k)/(v_c + v_k))^2`. +/// - `alpha = 1 - (nmatchval - 1) * sum(utcm * diff2) +/// / sum(nc_c * nc_k * diff2)` over the upper triangle (line 63). +/// +/// Documented deviations from R: an all-missing matrix is an error here +/// (R's line-45 path would report alpha = 1 with zero levels); infinite +/// ratings are rejected; the ratio metric errors when any level pair sums +/// to zero (R silently produces Inf/NaN); dimension caps. +/// +/// `ratings` is row-major raters x subjects; NaN marks missing. No +/// standard error or CI is produced (the R source computes none). +/// +/// # References +/// Gamer, M., Lemon, J., Fellows, I., & Singh, P. (2019). *irr: Various +/// coefficients of interrater reliability and agreement* (Version 0.85) +/// [Computer software]. CRAN. https://CRAN.R-project.org/package=irr +/// Krippendorff, K. (1980). *Content analysis: An introduction to its +/// methodology*. Sage. (as cited in Gamer et al., 2019; NOT READ) +pub fn kripp_alpha( + ratings: &[f64], + nraters: usize, + nsubjects: usize, + method: &str, +) -> Result { + if !matches!(method, "nominal" | "ordinal" | "interval" | "ratio") { + return Err( + "method must be one of \"nominal\", \"ordinal\", \"interval\", \"ratio\"".into(), + ); + } + if nraters < 2 { + return Err("kripp_alpha needs at least 2 raters".into()); + } + if nsubjects < 1 { + return Err("kripp_alpha needs at least 1 subject".into()); + } + if nraters > 10_000 || nsubjects > 1_000_000 { + return Err("kripp_alpha: dimensions exceed caps (raters <= 1e4, subjects <= 1e6)".into()); + } + if ratings.len() != nraters * nsubjects { + return Err(format!( + "ratings length {} does not match raters*subjects = {}", + ratings.len(), + nraters * nsubjects + )); + } + if ratings.iter().any(|v| v.is_infinite()) { + return Err("ratings must not contain infinities (use NaN for missing)".into()); + } + let mut levels: Vec = ratings.iter().copied().filter(|v| v.is_finite()).collect(); + levels.sort_by(|a, b| a.partial_cmp(b).expect("finite by filter")); + levels.dedup(); + let nval = levels.len(); + if nval == 0 { + return Err("kripp_alpha: all ratings are missing".into()); + } + let any_na = ratings.iter().any(|v| v.is_nan()); + let lev_index = |v: f64| -> usize { + levels + .binary_search_by(|p| p.partial_cmp(&v).expect("finite levels")) + .expect("observed value is a level by construction") + }; + let mut cm = vec![0.0_f64; nval * nval]; + for col in 0..nsubjects { + // R lines 12-13: per-column divisor only under the global-NA path. + let mc = if any_na { + let nonmiss = (0..nraters) + .filter(|&r| !ratings[r * nsubjects + col].is_nan()) + .count(); + nonmiss as f64 - 1.0 + } else { + 1.0 + }; + for i1 in 0..nraters - 1 { + for i2 in (i1 + 1)..nraters { + let a = ratings[i1 * nsubjects + col]; + let b = ratings[i2 * nsubjects + col]; + if a.is_nan() || b.is_nan() { + continue; + } + // A column visited here has >= 2 non-missing values, so + // mc >= 1 under the NA path (mc = 0 or -1 only occurs for + // columns that form no pair). + let (ia, ib) = (lev_index(a), lev_index(b)); + // R line 20: diagonal gains 2/mc, off-diagonal 1/mc with + // the mirror cell set by assignment (line 21). + let inc = if ia == ib { 2.0 } else { 1.0 } / mc; + cm[ia * nval + ib] += inc; + if ia != ib { + cm[ib * nval + ia] = cm[ia * nval + ib]; + } + } + } + } + let nmatchval: f64 = cm.iter().sum(); + let mut value = 1.0; + if nval >= 2 { + let nc: Vec = (0..nval) + .map(|i| cm[i * nval..(i + 1) * nval].iter().sum()) + .collect(); + let mut num = 0.0; + let mut den = 0.0; + for k in 1..nval { + for c in 0..k { + let diff2 = match method { + "nominal" => 1.0, + "ordinal" => { + let s: f64 = nc[c] / 2.0 + nc[c + 1..k].iter().sum::() + nc[k] / 2.0; + s * s + } + "interval" => { + let d = levels[c] - levels[k]; + d * d + } + _ => { + let s = levels[c] + levels[k]; + if s == 0.0 { + return Err( + "kripp_alpha: ratio metric undefined (level pair sums to zero)" + .into(), + ); + } + let d = (levels[c] - levels[k]) / s; + d * d + } + }; + num += cm[c * nval + k] * diff2; + den += nc[c] * nc[k] * diff2; + } + } + if den == 0.0 { + return Err("kripp_alpha: degenerate data (zero denominator)".into()); + } + value = 1.0 - (nmatchval - 1.0) * num / den; + } + if !value.is_finite() || !nmatchval.is_finite() { + return Err("kripp_alpha: degenerate data (non-finite result)".into()); + } + Ok(KrippResult { + value, + subjects: nsubjects as u64, + raters: nraters as u64, + levels: nval as u64, + nmatchval, + }) +} + #[cfg(test)] #[path = "../../../tests/unit/reliability_tests.rs"] mod tests; diff --git a/python/fast_mlsirm/__init__.py b/python/fast_mlsirm/__init__.py index ee58b407a..5983c3aef 100644 --- a/python/fast_mlsirm/__init__.py +++ b/python/fast_mlsirm/__init__.py @@ -125,6 +125,8 @@ SeparationReliabilityResult as SeparationReliabilityResult, icc as icc, IccResult as IccResult, + kripp_alpha as kripp_alpha, + KrippResult as KrippResult, ) from .factor import ( minres_fa as minres_fa, @@ -366,6 +368,8 @@ "SeparationReliabilityResult", "icc", "IccResult", + "kripp_alpha", + "KrippResult", "gtheory_pi", "gtheory_pio", "phi_lambda", diff --git a/python/fast_mlsirm/reliability.py b/python/fast_mlsirm/reliability.py index 40b67dfe6..2f19f89c8 100644 --- a/python/fast_mlsirm/reliability.py +++ b/python/fast_mlsirm/reliability.py @@ -415,4 +415,86 @@ def icc( p_value=float(res["p_value"]), lbound=float(res["lbound"]), ubound=float(res["ubound"]), - ) \ No newline at end of file + ) +@dataclass +class KrippResult: + """Krippendorff's alpha (irr ``kripp.alpha``). + + ``value`` is the alpha estimate; ``subjects``/``raters`` echo the + matrix dimensions as given; ``levels`` counts the distinct observed + rating values; ``nmatchval`` is the total coincidence-matrix mass + (R ``nmatchval``). No SE or CI is produced (the R source computes + none).""" + + value: float + subjects: int + raters: int + levels: int + nmatchval: float + + +def kripp_alpha(ratings, method: str = "nominal") -> KrippResult: + """Krippendorff's alpha for a raters x subjects matrix (compute in + Rust; transcribed from CRAN irr 0.85 ``R/kripp.alpha.R``, read in + full). ``ratings`` rows are raters ("classifiers"), columns are + subjects; NaN marks missing. ``method`` selects the distance metric: + "nominal", "ordinal", "interval", or "ratio". + + The irr source divides each column's pair counts by + ``#nonmissing - 1`` only when the matrix contains at least one + missing value, and by 1 otherwise; that quirk is preserved verbatim, + so complete-data alpha differs from the ``m - 1`` convention. + Documented deviations: an all-missing matrix, infinite ratings, and + a ratio-metric level pair summing to zero raise ``ValueError`` (R + would return alpha = 1, propagate, or emit Inf/NaN respectively). + In LLM-as-a-Judge quality management this estimates chance-corrected + agreement among judges over the same units. + + References (APA 7th ed.): + Gamer, M., Lemon, J., Fellows, I., & Singh, P. (2019). *irr: + Various coefficients of interrater reliability and agreement* + [R package]. https://CRAN.R-project.org/package=irr + Krippendorff, K. (1980). *Content analysis: An introduction to + its methodology*. Sage. (as cited in Gamer et al., 2019; + not read) + """ + from .fitstats import _core_module + + core = _core_module() + if core is None or not hasattr(core, "kripp_alpha"): + raise RuntimeError("kripp_alpha requires the compiled Rust core") + if isinstance(ratings, np.ma.MaskedArray): + raise ValueError("masked arrays are not supported; use NaN for missing") + arr = np.asarray(ratings) + if arr.dtype == object: + raise ValueError( + "object-dtype arrays are not supported; pass a numeric array" + ) + if np.iscomplexobj(arr): + raise ValueError("ratings must be real-valued") + if arr.dtype.kind == "b": + raise ValueError("ratings must be numeric, not boolean") + if arr.dtype.kind not in "fiu": + raise ValueError("ratings must be a numeric array") + if arr.ndim != 2: + raise ValueError("ratings must be a 2-D raters x subjects array") + if arr.dtype.kind in "iu" and arr.size: + # Levels are defined by exact f64 value identity; integers beyond + # 2**53 are not exactly representable and distinct rating labels + # would silently collapse during the float64 conversion. + lo, hi = int(arr.min()), int(arr.max()) + if hi > 2**53 or lo < -(2**53): + raise ValueError( + "integer ratings beyond 2**53 cannot be represented " + "exactly as float64; distinct levels would collapse" + ) + x = np.ascontiguousarray(arr, dtype=np.float64) + nr, ns = x.shape + res = core.kripp_alpha(x.reshape(-1), int(nr), int(ns), str(method)) + return KrippResult( + value=float(res["value"]), + subjects=int(res["subjects"]), + raters=int(res["raters"]), + levels=int(res["levels"]), + nmatchval=float(res["nmatchval"]), + ) diff --git a/tests/test_paper_features.py b/tests/test_paper_features.py index 82d9e085c..0a91d9a02 100644 --- a/tests/test_paper_features.py +++ b/tests/test_paper_features.py @@ -10536,4 +10536,99 @@ def test_validation(self): with pytest.raises(ValueError): icc(x, r0=True) with pytest.raises(ValueError): - icc(np.full((4, 3), 7.0)) # constant matrix -> degenerate \ No newline at end of file + icc(np.full((4, 3), 7.0)) # constant matrix -> degenerate +class TestKripp: + """kripp_alpha (irr kripp.alpha.R). Pins are exact Fractions from the + executed oracle; K1 matches the irr man page's published values. Every + assert reads KrippResult fields returned by the crate via the wrapper.""" + + @staticmethod + def k1(): + import numpy as np + + nan = float("nan") + cols = [ + [1, 1, nan, 1], [2, 2, 3, 2], [3, 3, 3, 3], [3, 3, 3, 3], + [2, 2, 2, 2], [1, 2, 3, 4], [4, 4, 4, 4], [1, 1, 2, 1], + [2, 2, 2, 2], [nan, 5, 5, 5], [nan, nan, 1, 1], [nan, nan, 3, nan], + ] + return np.array(cols, dtype=float).T # 4 raters x 12 subjects + + def test_anchor_k1(self): + from fast_mlsirm import kripp_alpha + + x = self.k1() + want = { + "nominal": 113 / 152, + "ordinal": 108577 / 133160, + "interval": 951 / 1120, + "ratio": 18222619 / 22852465, + } + for m, v in want.items(): + r = kripp_alpha(x, method=m) + assert abs(r.value - v) < 1e-12, m + assert r.nmatchval == 40.0 + assert r.subjects == 12 and r.raters == 4 and r.levels == 5 + + def test_no_na_quirk(self): + import numpy as np + from fast_mlsirm import kripp_alpha + + x = np.array([[1, 2, 3, 3, 2], [1, 2, 3, 3, 1], [2, 2, 3, 3, 2]], dtype=float) + r = kripp_alpha(x, method="nominal") + assert abs(r.value - 43 / 72) < 1e-12 + assert r.nmatchval == 30.0 # mc = 1 on complete data (irr quirk) + + def test_single_level(self): + import numpy as np + from fast_mlsirm import kripp_alpha + + r = kripp_alpha(np.full((2, 2), 2.0)) + assert r.value == 1.0 and r.levels == 1 + + def test_validation(self): + import numpy as np + import pytest + from fast_mlsirm import kripp_alpha + + x = self.k1() + with pytest.raises(ValueError): + kripp_alpha(x, method="euclid") + with pytest.raises(ValueError): + kripp_alpha(np.array([1.0, 2.0])) # 1-D + with pytest.raises(ValueError): + kripp_alpha(np.array([[1.0, np.inf], [2.0, 3.0]])) + with pytest.raises(ValueError): + kripp_alpha(np.full((2, 2), np.nan)) # all missing + with pytest.raises(ValueError): + kripp_alpha(np.array([[-1.0, 1.0], [1.0, -1.0]]), method="ratio") + with pytest.raises(ValueError): + kripp_alpha(np.ma.masked_array(x, mask=False)) + with pytest.raises(ValueError): + kripp_alpha(x.astype(complex)) + with pytest.raises(ValueError): + kripp_alpha(np.array([[True, False], [False, True]])) + with pytest.raises(ValueError): + kripp_alpha(np.array([[1.0, 2.0], [3.0, 4.0]], dtype=object)) + + def test_large_integer_levels_rejected(self): + # Review finding: int64 levels beyond 2**53 collapse in the float64 + # cast (2**53 and 2**53+1 map to the same f64), silently turning + # complete disagreement into alpha=1 with a single level. Reads the + # wrapper's rejection; the accepted-boundary asserts read crate output. + import numpy as np + import pytest + from fast_mlsirm import kripp_alpha + + big = 2**53 + x = np.array([[big, big + 1], [big + 1, big]], dtype=np.int64) + with pytest.raises(ValueError, match="2\\*\\*53"): + kripp_alpha(x) + with pytest.raises(ValueError, match="2\\*\\*53"): + kripp_alpha(np.array([[2**63, 2**63 + 2]], dtype=np.uint64).reshape(2, 1).T) + with pytest.raises(ValueError, match="2\\*\\*53"): + kripp_alpha(-x) # negative side + # Boundary: values at exactly +/-2**53 are exactly representable. + ok = np.array([[big, -big], [-big, big]], dtype=np.int64) + r = kripp_alpha(ok) + assert r.levels == 2 and abs(r.value - (-0.5)) < 1e-12 diff --git a/tests/unit/reliability_tests.rs b/tests/unit/reliability_tests.rs index a55f2e9fe..f2f5b30ef 100644 --- a/tests/unit/reliability_tests.rs +++ b/tests/unit/reliability_tests.rs @@ -981,3 +981,196 @@ fn ic_mc_500_invariance() { } } } + +// --------------------------------------------------------------------------- +// kripp_alpha (irr 0.85 kripp.alpha.R). Pins are exact Fractions from the +// executed oracle (files/kripp_oracle.py); K1 matches the irr man page's +// published values. Every assert reads KrippResult fields returned by the +// crate. +// --------------------------------------------------------------------------- + +fn ka_close(a: f64, b: f64) -> bool { + let d = (a - b).abs(); + d <= 1e-12 || d <= 1e-12 * b.abs() +} + +/// irr man-page `nmm` matrix, raters x subjects (4 x 12); NaN = missing. +fn ka_k1() -> Vec { + let nan = f64::NAN; + // Column-major R matrix(nrow=4) transposed to rater rows. + let cols: [[f64; 4]; 12] = [ + [1.0, 1.0, nan, 1.0], + [2.0, 2.0, 3.0, 2.0], + [3.0, 3.0, 3.0, 3.0], + [3.0, 3.0, 3.0, 3.0], + [2.0, 2.0, 2.0, 2.0], + [1.0, 2.0, 3.0, 4.0], + [4.0, 4.0, 4.0, 4.0], + [1.0, 1.0, 2.0, 1.0], + [2.0, 2.0, 2.0, 2.0], + [nan, 5.0, 5.0, 5.0], + [nan, nan, 1.0, 1.0], + [nan, nan, 3.0, nan], + ]; + let mut x = vec![0.0; 4 * 12]; + for (c, col) in cols.iter().enumerate() { + for (r, v) in col.iter().enumerate() { + x[r * 12 + c] = *v; + } + } + x +} + +#[test] +fn ka_anchor_k1_all_methods() { + // Kills MU1 (diag 1/mc: nominal -> 501/869 = 0.5765...), MU3 (ordinal + // full weights -> 167915/216366), MU4 (interval |d| -> 417/521), + // MU5 (nmv excludes diagonal -> 145/152), MU6 (num*nc products). + let x = ka_k1(); + let want = [ + ("nominal", 113.0 / 152.0), + ("ordinal", 108577.0 / 133160.0), + ("interval", 951.0 / 1120.0), + ("ratio", 18222619.0 / 22852465.0), + ]; + for (m, v) in want { + let r = kripp_alpha(&x, 4, 12, m).unwrap(); + assert!(ka_close(r.value, v), "{m}: {} vs {v}", r.value); + assert_eq!(r.nmatchval, 40.0, "{m} nmv"); + assert_eq!(r.subjects, 12); + assert_eq!(r.raters, 4); + assert_eq!(r.levels, 5); + } +} + +#[test] +fn ka_no_na_quirk_k2() { + // Complete data uses divisor mc = 1, NOT m - 1 (R lines 12-13 quirk). + // Kills MU2 (mc always m-1): mutant nominal would be 11/18 != 43/72. + let x = [ + 1.0, 2.0, 3.0, 3.0, 2.0, // + 1.0, 2.0, 3.0, 3.0, 1.0, // + 2.0, 2.0, 3.0, 3.0, 2.0, + ]; + let n = kripp_alpha(&x, 3, 5, "nominal").unwrap(); + assert!(ka_close(n.value, 43.0 / 72.0), "nominal {}", n.value); + assert_eq!(n.nmatchval, 30.0); + let i = kripp_alpha(&x, 3, 5, "interval").unwrap(); + assert!(ka_close(i.value, 97.0 / 126.0), "interval {}", i.value); + // Guard the guard: the MU2 mutant value differs from the true pin. + assert!(!ka_close(n.value, 11.0 / 18.0)); +} + +#[test] +fn ka_hand_fixture_k3() { + let x = [1.0, 2.0, 3.0, 1.0, 1.0, 3.0, 3.0, 2.0]; + let n = kripp_alpha(&x, 2, 4, "nominal").unwrap(); + assert!(ka_close(n.value, 1.0 / 3.0), "nominal {}", n.value); + let o = kripp_alpha(&x, 2, 4, "ordinal").unwrap(); + assert!(ka_close(o.value, 17.0 / 24.0), "ordinal {}", o.value); + // Interval coincidentally equals ordinal on this fixture (both crate + // outputs; the ordinal/interval swap mutant is killed on K1 where the + // two pins differ). + let iv = kripp_alpha(&x, 2, 4, "interval").unwrap(); + assert!(ka_close(iv.value, 17.0 / 24.0), "interval {}", iv.value); + assert!(ka_close(iv.value, o.value)); + let rt = kripp_alpha(&x, 2, 4, "ratio").unwrap(); + assert!(ka_close(rt.value, 1889.0 / 2841.0), "ratio {}", rt.value); + assert_eq!(n.nmatchval, 8.0); + assert_eq!(n.levels, 3); +} + +#[test] +fn ka_single_level_k4() { + // R line 45: fewer than 2 levels -> alpha = 1. + let x = [2.0, 2.0, 2.0, 2.0]; + let r = kripp_alpha(&x, 2, 2, "nominal").unwrap(); + assert_eq!(r.value, 1.0); + assert_eq!(r.nmatchval, 4.0); + assert_eq!(r.levels, 1); +} + +#[test] +fn ka_error_contract() { + let ok = [1.0, 2.0, 2.0, 1.0]; + assert!(kripp_alpha(&ok, 2, 2, "euclid") + .unwrap_err() + .contains("method")); + assert!(kripp_alpha(&ok, 1, 4, "nominal") + .unwrap_err() + .contains("raters")); + assert!(kripp_alpha(&[], 2, 0, "nominal") + .unwrap_err() + .contains("subject")); + assert!(kripp_alpha(&ok, 2, 3, "nominal") + .unwrap_err() + .contains("length")); + let inf = [1.0, f64::INFINITY, 2.0, 1.0]; + assert!(kripp_alpha(&inf, 2, 2, "nominal") + .unwrap_err() + .contains("infinit")); + let nan = f64::NAN; + assert!(kripp_alpha(&[nan, nan, nan, nan], 2, 2, "nominal") + .unwrap_err() + .contains("missing")); + // Ratio metric undefined when a level pair sums to zero. + let zsum = [-1.0, 1.0, 1.0, -1.0]; + assert!(kripp_alpha(&zsum, 2, 2, "ratio") + .unwrap_err() + .contains("ratio")); + // Same data is fine for nominal/interval (crate outputs finite). + assert!(kripp_alpha(&zsum, 2, 2, "nominal") + .unwrap() + .value + .is_finite()); + assert!(kripp_alpha(&zsum, 2, 2, "interval") + .unwrap() + .value + .is_finite()); +} + +#[test] +#[ignore = "MC-500: run explicitly with cargo test -- --ignored"] +fn ka_mc_500_permutation_invariance() { + // Alpha is invariant under rater-row and subject-column permutation; + // both sides of every comparison are crate outputs. + let mut state = 0xCA5EEDu64; + let mut next = |s: &mut u64| -> u64 { + *s = s + .wrapping_mul(6364136223846793005) + .wrapping_add(1442695040888963407); + *s >> 11 + }; + let base = ka_k1(); + let methods = ["nominal", "ordinal", "interval", "ratio"]; + for rep in 0..500 { + let m = methods[rep % 4]; + let want = kripp_alpha(&base, 4, 12, m).unwrap(); + // Fisher-Yates over rater rows. + let mut rows: Vec = (0..4).collect(); + for i in (1..4).rev() { + let j = (next(&mut state) % (i as u64 + 1)) as usize; + rows.swap(i, j); + } + let mut cols: Vec = (0..12).collect(); + for i in (1..12).rev() { + let j = (next(&mut state) % (i as u64 + 1)) as usize; + cols.swap(i, j); + } + let mut xp = vec![0.0; 48]; + for (rn, &ro) in rows.iter().enumerate() { + for (cn, &co) in cols.iter().enumerate() { + xp[rn * 12 + cn] = base[ro * 12 + co]; + } + } + let got = kripp_alpha(&xp, 4, 12, m).unwrap(); + assert!( + ka_close(got.value, want.value), + "rep {rep} {m}: {} vs {}", + got.value, + want.value + ); + assert_eq!(got.nmatchval, want.nmatchval, "rep {rep} nmv"); + assert_eq!(got.levels, want.levels, "rep {rep} levels"); + } +}