Add icc intraclass correlation coefficients (Shrout-Fleiss taxonomy, irr 0.85)

CHANGELOG.md

@@ -119,6 +119,9 @@
 
 ### 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.
 - `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.

crates/fast-mlsirm-py/src/lib.rs

@@ -109,8 +109,8 @@ 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,
-    separation_reliability as core_separation_reliability,
+    cronbach_alpha as core_cronbach_alpha, feldt_alpha_ci as core_feldt_alpha_ci, icc as core_icc,
+    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::{
@@ -3578,6 +3578,71 @@ 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<Py<pyo3::types::PyDict>> {
+    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())
+}
+
+/// 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<Py<pyo3::types::PyDict>> {
+    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`.
@@ -5496,13 +5561,17 @@ fn glicko2_rating(
         .iter()
         .map(|&v| usize::try_from(v))
         .collect::<Result<_, _>>()
-        .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<usize> = black
         .as_slice()?
         .iter()
         .map(|&v| usize::try_from(v))
         .collect::<Result<_, _>>()
-        .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 +5725,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 +5879,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
@@ -5848,6 +5915,38 @@ 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<Py<pyo3::types::PyDict>> {
+    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]
@@ -8148,6 +8247,8 @@ 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!(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)?)?;
@@ -8205,6 +8306,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)?)?;