Add bratt_mm: Bradley-Terry model with additive ties parameter (VGAM bratt)

CHANGELOG.md

@@ -119,6 +119,9 @@
 
 ### 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.
 - `metrics_rating` prediction-quality metrics for binary-outcome forecasts (binomial deviance on capped predictions, RMSE/MAE on raw predictions, optional 0.5-baseline scaling), a Rust reimplementation of CRAN PlayerRatings 1.1-0 `metrics()` with its cap quirk and elementwise NaN semantics preserved; exact-oracle anchor tests, mutation-kill map, and MC-500 invariants.
 - **Multiplayer Elo rating (CRAN PlayerRatings 1.1-0 `elom()`

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

@@ -5741,6 +5741,142 @@ fn fide_rating(
     Ok(d.into())
 }
 
+/// Predicted game outcomes from fitted ratings, PlayerRatings
+/// `predict.rating` two-player branches (see
+/// `mlsirm_core::scaling::predict_rating_two`). `white`/`black` are player
+/// indices with -1 = unmatched. Pass `deviations` for the
+/// Glicko/Glicko-2/Stephenson deviation-shrunk branch.
+#[pyfunction]
+#[pyo3(signature = (ratings, deviations, games, white, black, gamma, tng, trat_rating, trat_deviation, thresh))]
+#[allow(clippy::too_many_arguments)]
+fn predict_rating_two(
+    py: Python<'_>,
+    ratings: PyReadonlyArray1<'_, f64>,
+    deviations: Option<PyReadonlyArray1<'_, f64>>,
+    games: PyReadonlyArray1<'_, u64>,
+    white: PyReadonlyArray1<'_, i64>,
+    black: PyReadonlyArray1<'_, i64>,
+    gamma: PyReadonlyArray1<'_, f64>,
+    tng: u64,
+    trat_rating: Option<f64>,
+    trat_deviation: Option<f64>,
+    thresh: Option<f64>,
+) -> PyResult<Py<PyArray1<f64>>> {
+    let dev_slice = deviations.as_ref().map(|d| d.as_slice()).transpose()?;
+    let trat = trat_rating.map(|t1| (t1, trat_deviation.unwrap_or(f64::NAN)));
+    let out = mlsirm_core::scaling::predict_rating_two(
+        ratings.as_slice()?,
+        dev_slice,
+        games.as_slice()?,
+        white.as_slice()?,
+        black.as_slice()?,
+        gamma.as_slice()?,
+        tng,
+        trat,
+        thresh,
+    )
+    .map_err(PyValueError::new_err)?;
+    Ok(PyArray1::from_slice(py, &out).into())
+}
+
+/// Predicted expected scores for multi-player (EloM) events, PlayerRatings
+/// `predict.rating` EloM branch (see
+/// `mlsirm_core::scaling::predict_rating_multi`). `players` is flattened
+/// row-major nr x np with -1 = empty seat; `placing` returns min-tie ranks.
+#[pyfunction]
+#[pyo3(signature = (ratings, games, players, nr, np, tng, trat, placing))]
+#[allow(clippy::too_many_arguments)]
+fn predict_rating_multi(
+    py: Python<'_>,
+    ratings: PyReadonlyArray1<'_, f64>,
+    games: PyReadonlyArray1<'_, u64>,
+    players: PyReadonlyArray1<'_, i64>,
+    nr: usize,
+    np: usize,
+    tng: u64,
+    trat: Option<f64>,
+    placing: bool,
+) -> PyResult<Py<PyArray1<f64>>> {
+    let out = mlsirm_core::scaling::predict_rating_multi(
+        ratings.as_slice()?,
+        games.as_slice()?,
+        players.as_slice()?,
+        nr,
+        np,
+        tng,
+        trat,
+        placing,
+    )
+    .map_err(PyValueError::new_err)?;
+    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
+/// (worths, alpha[ref_index] == ref_value), alpha0 (tie parameter),
+/// iterations, log_likelihood.
+#[pyfunction]
+#[pyo3(signature = (wins, ties, n, ref_index=0, ref_value=1.0, max_iter=10000, tol=1e-10))]
+#[allow(clippy::too_many_arguments)]
+fn bratt_mm(
+    py: Python<'_>,
+    wins: PyReadonlyArray1<'_, f64>,
+    ties: PyReadonlyArray1<'_, f64>,
+    n: usize,
+    ref_index: usize,
+    ref_value: f64,
+    max_iter: usize,
+    tol: f64,
+) -> PyResult<Py<pyo3::types::PyDict>> {
+    let res = mlsirm_core::scaling::bratt_mm(
+        wins.as_slice()?,
+        ties.as_slice()?,
+        n,
+        ref_index,
+        ref_value,
+        max_iter,
+        tol,
+    )
+    .map_err(PyValueError::new_err)?;
+    let d = pyo3::types::PyDict::new(py);
+    d.set_item("alpha", PyArray1::from_slice(py, &res.alpha))?;
+    d.set_item("alpha0", res.alpha0)?;
+    d.set_item("iterations", res.iterations as u64)?;
+    d.set_item("log_likelihood", res.log_likelihood)?;
+    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]
@@ -8095,6 +8231,10 @@ fn fast_mlsirm_core(m: &Bound<'_, PyModule>) -> PyResult<()> {
     m.add_function(wrap_pyfunction!(elom_rating, m)?)?;
     m.add_function(wrap_pyfunction!(metrics_rating, m)?)?;
     m.add_function(wrap_pyfunction!(fide_rating, m)?)?;
+    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)?)?;

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<f64>,
+    pub category_z: Vec<f64>,
+    pub category_p: Vec<f64>,
+}
+
+/// 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<FleissKappaResult, String> {
+    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<f64>> = 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::<f64>() - nrf) / (nrf * (nrf - 1.0)))
+        .sum::<f64>()
+        / mf;
+
+    let col: Vec<f64> = (0..k)
+        .map(|j| ttab.iter().map(|t| t[j]).sum::<f64>())
+        .collect();
+    let p: Vec<f64> = 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<f64> = (0..nr).map(|r| rater_counts[r][j] / mf).collect();
+            let mean = props.iter().sum::<f64>() / nrf;
+            let s2 = props.iter().map(|&v| (v - mean) * (v - mean)).sum::<f64>() / (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::<f64>())
+        / (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;