diff --git a/crates/fast-mlsirm-py/src/lib.rs b/crates/fast-mlsirm-py/src/lib.rs index 607a97c03..91ed1d497 100644 --- a/crates/fast-mlsirm-py/src/lib.rs +++ b/crates/fast-mlsirm-py/src/lib.rs @@ -1828,9 +1828,9 @@ fn fit_bifactor_grm_fipc( } /// Single-group polytomous two-tier graded response model (Cai, 2010, -/// abstract read; Cai, Yang, & Hansen, 2011, eq. 6-7, full text read; +/// pp. 583-584, full text read; Cai, Yang, & Hansen, 2011, eq. 6-7, full text read; /// `mlsirm_core::two_tier_grm::fit_two_tier_grm`). Each item's `n_cat` ORDERED categories load a caller-supplied subset of -/// the `n_primary` correlated primary dimensions (`a_primary`, row-major +/// the `n_primary` primary dimensions (`a_primary`, row-major /// `n_items * n_primary`, unconstrained, `0` at fixed pattern positions) /// and at most one orthogonal specific factor (`a_specific`, /// unconstrained, `0` for specific-free items): @@ -1858,15 +1858,14 @@ fn fit_bifactor_grm_fipc( /// /// Cai, L. (2010). A two-tier full-information item factor analysis model /// with applications. *Psychometrika, 75*(4), 581-612. -/// https://doi.org/10.1007/s11336-010-9178-0 (abstract read; full text not -/// accessible — no equation locator is drawn from it) +/// https://doi.org/10.1007/s11336-010-9178-0 (full text read, pp. 583-584) /// /// Cai, L., Yang, J. S., & Hansen, M. (2011). Generalized full-information /// item bifactor analysis. *Psychological Methods, 16*(3), 221-248. /// https://doi.org/10.1037/a0023350 (full text read) #[pyfunction] #[allow(clippy::too_many_arguments)] -#[pyo3(signature = (y, observed, primary_map, specific_map, n_persons, n_items, n_primary, n_specific, n_cat, q_primary = 15, q_specific = 11, max_iter = 500, tol = 1e-6, n_starts = 1, seed = 0x9E37_79B9_7F4A_7C15))] +#[pyo3(signature = (y, observed, primary_map, specific_map, n_persons, n_items, n_primary, n_specific, n_cat, q_primary = 15, q_specific = 11, max_iter = 500, tol = 1e-6, n_starts = 1, seed = 0x9E37_79B9_7F4A_7C15, primary_correlation = "estimate"))] fn fit_two_tier_grm( py: Python<'_>, y: PyReadonlyArray1<'_, i64>, @@ -1884,7 +1883,17 @@ fn fit_two_tier_grm( tol: f64, n_starts: usize, seed: u64, + primary_correlation: &str, ) -> PyResult> { + let estimate_primary_correlation = match primary_correlation { + "estimate" => true, + "identity" => false, + _ => { + return Err(PyValueError::new_err( + "primary_correlation must be 'estimate' or 'identity'", + )) + } + }; let y_slice = y.as_slice()?; let obs_vec: Option> = match &observed { Some(o) => Some(o.as_slice()?.to_vec()), @@ -1913,6 +1922,7 @@ fn fit_two_tier_grm( }) .collect::>()?; let cfg = TwoTierGrmConfig { + estimate_primary_correlation, q_primary, q_specific, max_iter, @@ -1955,12 +1965,14 @@ fn fit_two_tier_grm( out.set_item("final_loglik_change", res.final_loglik_change)?; out.set_item("best_start", res.best_start)?; out.set_item("n_parameters", res.n_parameters)?; + out.set_item("primary_identification", res.primary_identification)?; Ok(out.into()) } /// Observed-information SEs for the confirmatory two-tier GRM via Oakes /// (1999, eq. 6, p. 480). Free vector: free primary slopes, optional -/// specific slope, thresholds, Fisher-`z` primary correlations. +/// specific slope, thresholds, and Fisher-`z` primary correlations only when +/// `primary_correlation="estimate"` (Cai, 2010, pp. 583-584). /// `q_primary`/`q_specific`/`fd_step` are REQUIRED (ADR-0028 / #1929). /// Non-PD information returns `se=None` (never substituted). /// @@ -1969,12 +1981,16 @@ fn fit_two_tier_grm( /// Statistical Society Series B: Statistical Methodology, 61*(2), 479-482. /// https://doi.org/10.1111/1467-9868.00188; Cai, L., Yang, J. S., & Hansen, /// M. (2011). Generalized full-information item bifactor analysis. -/// *Psychological Methods, 16*(3), 221-248. https://doi.org/10.1037/a0023350 +/// *Psychological Methods, 16*(3), 221-248. https://doi.org/10.1037/a0023350; +/// Cai, L. (2010). A two-tier full-information item factor analysis model +/// with applications. *Psychometrika, 75*(4), 581-612. +/// https://doi.org/10.1007/s11336-010-9178-0 (pp. 583-584) #[pyfunction] #[allow(clippy::too_many_arguments)] #[pyo3(signature = ( a_primary, a_specific, threshold, phi, y, observed, primary_map, specific_map, - n_persons, n_items, n_primary, n_specific, n_cat, q_primary, q_specific, fd_step + n_persons, n_items, n_primary, n_specific, n_cat, q_primary, q_specific, fd_step, + primary_correlation = "estimate" ))] fn two_tier_oakes_se( py: Python<'_>, @@ -1994,7 +2010,17 @@ fn two_tier_oakes_se( q_primary: usize, q_specific: usize, fd_step: f64, + primary_correlation: &str, ) -> PyResult> { + let estimate_primary_correlation = match primary_correlation { + "estimate" => true, + "identity" => false, + _ => { + return Err(PyValueError::new_err( + "primary_correlation must be 'estimate' or 'identity'", + )) + } + }; let y_slice = y.as_slice()?; let obs_vec: Option> = match &observed { Some(o) => Some(o.as_slice()?.to_vec()), @@ -2034,6 +2060,7 @@ fn two_tier_oakes_se( let thr = threshold.as_slice()?.to_vec(); let phi_vec = phi.as_slice()?.to_vec(); let cfg = TwoTierOakesConfig { + estimate_primary_correlation, q_primary, q_specific, fd_step, diff --git a/crates/mlsirm-core/src/two_tier_grm.rs b/crates/mlsirm-core/src/two_tier_grm.rs index 784ce925b..176613d5d 100644 --- a/crates/mlsirm-core/src/two_tier_grm.rs +++ b/crates/mlsirm-core/src/two_tier_grm.rs @@ -45,30 +45,7 @@ //! below means degenerate single-loader inputs valid in stage 1 are rejected //! here — reduction equivalence holds for full-pattern inputs). The two-tier model itself — //! correlated primaries plus orthogonal specifics, subsuming the bifactor -//! and testlet models — is Cai (2010) (abstract; see the source-access note -//! below). -//! -//! # Source-access note (why Cai 2010 has no equation locator here) -//! -//! The Zotero record for Cai (2010) (key `GT3NQ8K8`) holds the abstract, -//! DOI (`10.1007/s11336-010-9178-0`), and bibliographic data — its abstract -//! confirms the model claims used here (the framework "subsumes standard -//! multidimensional IRT models, bifactor IRT models, and testlet response -//! theory models as special cases", "reduction in the dimensionality of the -//! latent variable space", "an EM algorithm for full-information maximum -//! marginal likelihood estimation") — but the attached file is the Springer -//! article landing page, NOT the full text, and the full text could not be -//! obtained in-run (paywalled at the publisher; the institutional proxy -//! serves the article page without entitlement, Springer WAYF rejects the -//! proxy redirect host, and no author manuscript was found). So no Cai -//! (2010) equation or page number is cited: every locator below names a -//! source whose full text was actually read — the local Gibbons et al. -//! (2007) PDF (eq. 9, 11-12, 15), the open-access Cai, Yang, & Hansen -//! (2011) full text (eq. 6-7, 10-11, "Maximum Marginal Likelihood -//! Estimation" section), and the installed mirt 1.46.1 `bfactor` help topic -//! (two-tier covariance, `ncol(G) + 1` integration). The `mirt::bfactor` -//! two-tier oracle, whose own implementation follows Cai (2010), validates -//! the same MLE empirically. +//! and testlet models — is Cai (2010, pp. 583-584). //! //! # Estimation: Bock-Aitkin EM with reduction over the specific tier //! @@ -159,6 +136,10 @@ //! block — negating that dimension's slopes AND the reported primary EAP //! column AND the `Phi` row/column signs for primary flips, but NOT the //! thresholds. +//! When `Phi = I`, identical free-loading item sets for two primary columns +//! leave their orthogonal rotation unidentified. Distinct supports are a +//! necessary condition for the fixed-identity specialization (Cai, 2010, +//! pp. 583-584); the per-column reflection rule above handles signs. //! //! # Caller-owned numerics (no hidden clamps, no magic caps) //! @@ -195,8 +176,7 @@ //! //! Cai, L. (2010). A two-tier full-information item factor analysis model //! with applications. *Psychometrika, 75*(4), 581-612. -//! https://doi.org/10.1007/s11336-010-9178-0 (abstract + metadata read via -//! the Zotero record; full text not accessible — see the source-access note) +//! https://doi.org/10.1007/s11336-010-9178-0 (full text read, pp. 583-584) //! //! Cai, L., Yang, J. S., & Hansen, M. (2011). Generalized full-information //! item bifactor analysis. *Psychological Methods, 16*(3), 221-248. @@ -246,6 +226,8 @@ use crate::poly::{grm_logprobs, grm_node_gradient, solve_small}; /// fixed-table cap), so this module imposes no upper cap of its own. #[derive(Clone, Copy, Debug)] pub struct TwoTierGrmConfig { + /// Estimate primary correlations; false fixes Phi to the identity. + pub estimate_primary_correlation: bool, /// Gauss-Hermite nodes per primary dimension (any `n >= 1`). /// The primary product grid has `q_primary^n_primary` nodes. pub q_primary: usize, @@ -279,8 +261,8 @@ pub struct TwoTierGrmResult { pub a_specific: Vec, /// Ordered boundary intercepts `d_ik`, row-major `n_items * (n_cat - 1)`. pub threshold: Vec, - /// Estimated primary correlation matrix, row-major - /// `n_primary * n_primary` (unit diagonal). + /// Primary correlation matrix, row-major `n_primary * n_primary`; + /// exactly I when `estimate_primary_correlation` is false. pub phi: Vec, /// Primary-factor EAPs `E[theta_d | Y_p]`, row-major /// `n_persons * n_primary`. @@ -297,9 +279,12 @@ pub struct TwoTierGrmResult { pub final_loglik_change: f64, /// Winning start index in `0..n_starts` (deterministic from `seed`). pub best_start: usize, - /// `sum_i (k_i + has_specific(i) + (n_cat - 1)) + P*(P-1)/2` free - /// parameters, where `k_i` is item `i`'s free primary-slope count. + /// `sum_i (k_i + has_specific(i) + (n_cat - 1))` free item parameters + /// (`k_i` is item `i`'s free primary-slope count), plus `P*(P-1)/2` + /// when primary correlations are estimated. pub n_parameters: usize, + /// `"correlated"` when Phi was estimated, `"orthogonal"` when fixed to I. + pub primary_identification: &'static str, } /// Validated problem structure shared by the fitter and the public @@ -420,6 +405,19 @@ pub(crate) fn validate( )); } } + if !cfg.estimate_primary_correlation { + for d in 0..n_primary { + for other in d + 1..n_primary { + if (0..n_items) + .all(|i| primary_map[i * n_primary + d] == primary_map[i * n_primary + other]) + { + return Err(format!( + "identity primary correlation requires distinct free-loading item sets for every pair of primary columns; identical free-loading item sets at columns {d} and {other} admit orthogonal rotation" + )); + } + } + } + } let mut blocks: Vec> = vec![Vec::new(); n_specific]; let mut specific_free = Vec::new(); let mut item_block = vec![None; n_items]; @@ -541,6 +539,7 @@ fn initial_params( observed: Option<&[bool]>, seed: u64, start: usize, + estimate_phi: bool, ) -> (Vec, Vec) { let is_obs = |p: usize, i: usize| observed.is_none_or(|o| o[p * v.n_items + i]); let mut rng = SplitMix64(seed ^ (0x9E37_79B9_7F4A_7C15u64.wrapping_mul(start as u64 + 1))); @@ -583,7 +582,7 @@ fn initial_params( // Primary correlations in Fisher-z space (start 0: Phi = I). let m = v.n_primary * (v.n_primary.saturating_sub(1)) / 2; let mut z = vec![0.0f64; m]; - if start > 0 { + if start > 0 && estimate_phi { for slot in z.iter_mut() { *slot = 0.25 * rng.standard_normal(); } @@ -711,7 +710,12 @@ pub(crate) fn z_from_phi(phi: &[f64], p: usize) -> Result, String> { /// rule order). Fixed-grid Gauss-Hermite quadrature is an implementation /// choice (the embedded rules in `crate::quadrature`); the node counts are /// caller arguments with no upper cap in this module. -pub(crate) fn build_primary_grid(tz: &[f64], wz: &[f64], p: usize, n_grid: usize) -> (Vec, Vec) { +pub(crate) fn build_primary_grid( + tz: &[f64], + wz: &[f64], + p: usize, + n_grid: usize, +) -> (Vec, Vec) { let q = tz.len(); let mut coords = vec![0.0f64; n_grid * p]; let mut log_w0 = vec![0.0f64; n_grid]; @@ -1330,7 +1334,14 @@ fn run_single_start( start: usize, ) -> Result { let p = v.n_primary; - let (mut params, mut z_phi) = initial_params(v, y, observed, cfg.seed, start); + let (mut params, mut z_phi) = initial_params( + v, + y, + observed, + cfg.seed, + start, + cfg.estimate_primary_correlation, + ); // No pre-allocation from `max_iter`: it is caller-owned and unbounded // above, so `with_capacity(max_iter + 1)` could overflow; the trace grows @@ -1347,8 +1358,9 @@ fn run_single_start( .ok_or_else(|| format!("primary correlation became non-PD at iteration {n_iter}"))?; let phi_inv = chol_inverse(&l, p); let log_w = reweighted_log_weights(log_w0, coords, &phi_inv, logdet, p); - let (ll, counts, s_bar_sum) = - e_step(v, y, observed, ¶ms, &log_w, log_ws, coords, ts, n_grid, qs); + let (ll, counts, s_bar_sum) = e_step( + v, y, observed, ¶ms, &log_w, log_ws, coords, ts, n_grid, qs, + ); let previous = loglik_trace.last().copied(); let change = checked_em_loglik_change(ll, previous, n_iter)?; loglik_trace.push(ll); @@ -1368,7 +1380,9 @@ fn run_single_start( for slot in s_bar.iter_mut() { *slot /= v.n_persons as f64; } - z_phi = m_step_phi(z_phi, p, &s_bar, v.n_persons, cfg.ridge, cfg.newton_iter); + if cfg.estimate_primary_correlation { + z_phi = m_step_phi(z_phi, p, &s_bar, v.n_persons, cfg.ridge, cfg.newton_iter); + } for i in 0..v.n_items { let free = &v.free_primaries[i]; let has_specific = v.item_block[i].is_some(); @@ -1426,7 +1440,9 @@ fn run_single_start( /// `primary_map` is row-major `n_items * n_primary` confirmatory /// free-slope pattern; `specific_map` is length `n_items` with `-1` for /// specific-free items and `0..n_specific` otherwise. Runs `n_starts` EM -/// runs and keeps the best loglik. Returns `Err` on malformed input, +/// runs and keeps the best loglik. `estimate_primary_correlation=false` +/// fixes Phi to I (Cai, 2010, pp. 583-584), omitting its Fisher-z M-step; +/// true preserves estimated Phi. Returns `Err` on malformed input or /// unobserved categories (unidentified ordered boundary pair under Cai et /// al., 2011, eq. 7), or total numerical failure; per-start /// non-convergence is reported through the winning run's flags, never @@ -1436,8 +1452,7 @@ fn run_single_start( /// /// Cai, L. (2010). A two-tier full-information item factor analysis model /// with applications. *Psychometrika, 75*(4), 581-612. -/// https://doi.org/10.1007/s11336-010-9178-0 (abstract read; full text not -/// accessible — see the module source-access note) +/// https://doi.org/10.1007/s11336-010-9178-0 (full text read, pp. 583-584) /// /// Cai, L., Yang, J. S., & Hansen, M. (2011). Generalized full-information /// item bifactor analysis. *Psychological Methods, 16*(3), 221-248. @@ -1593,7 +1608,11 @@ pub fn fit_two_tier_grm( let mut a_primary = vec![0.0f64; n_items * p]; let mut a_specific = vec![0.0f64; n_items]; let mut threshold = vec![0.0f64; n_items * v.m1]; - let mut n_parameters = p * (p.saturating_sub(1)) / 2; + let mut n_parameters = if cfg.estimate_primary_correlation { + p * (p.saturating_sub(1)) / 2 + } else { + 0 + }; for (i, par) in params.iter().enumerate() { for &dim in &v.free_primaries[i] { a_primary[i * p + dim] = par.a_p[dim]; @@ -1609,7 +1628,8 @@ pub fn fit_two_tier_grm( // Per-dimension reflection canonicalization (module docs): each primary // over its loading items (flipping slopes, the EAP column, and the Phi - // row/column signs jointly), each specific within its block; + // row/column signs jointly only when Phi is estimated), each specific + // within its block; fixed Phi remains bit-exact I; // thresholds untouched. let mut phi_work = phi; for d in 0..p { @@ -1632,10 +1652,12 @@ pub fn fit_two_tier_grm( for pp in 0..n_persons { theta_p_eap[pp * p + d] = -theta_p_eap[pp * p + d]; } - for q in 0..p { - if q != d { - phi_work[d * p + q] = -phi_work[d * p + q]; - phi_work[q * p + d] = -phi_work[q * p + d]; + if cfg.estimate_primary_correlation { + for q in 0..p { + if q != d { + phi_work[d * p + q] = -phi_work[d * p + q]; + phi_work[q * p + d] = -phi_work[q * p + d]; + } } } } @@ -1677,6 +1699,11 @@ pub fn fit_two_tier_grm( final_loglik_change: outcome.final_loglik_change, best_start, n_parameters, + primary_identification: if cfg.estimate_primary_correlation { + "correlated" + } else { + "orthogonal" + }, }) } @@ -1723,6 +1750,7 @@ pub fn two_tier_grm_marginal_loglik( // `validate` sees them only for its own field-level bounds checks. // (No `..Default()` exists: Project rule, issue #1929.) let cfg = TwoTierGrmConfig { + estimate_primary_correlation: true, q_primary, q_specific, max_iter: 1, @@ -1815,6 +1843,7 @@ pub fn two_tier_grm_marginal_loglik_brute( // `validate` sees them only for its own field-level bounds checks. // (No `..Default()` exists: Project rule, issue #1929.) let cfg = TwoTierGrmConfig { + estimate_primary_correlation: true, q_primary, q_specific, max_iter: 1, diff --git a/crates/mlsirm-core/src/two_tier_oakes.rs b/crates/mlsirm-core/src/two_tier_oakes.rs index 588e6170b..dc3797c8f 100644 --- a/crates/mlsirm-core/src/two_tier_oakes.rs +++ b/crates/mlsirm-core/src/two_tier_oakes.rs @@ -28,8 +28,7 @@ //! # Model and free vector //! //! The response model is the confirmatory two-tier GRM of Cai (2010) -//! (abstract; full text not accessible — locators below use Cai et al., -//! 2011, and Gibbons et al., 2007, which were read in full): correlated +//! (full text read, pp. 583-584): correlated //! primaries plus orthogonal specifics with dimension reduction over the //! specific tier (Gibbons et al., 2007, eq. 15; Chalmers, 2026, mirt //! `bfactor` documentation: `ncol(G) + 1` integration). The free vector is @@ -53,8 +52,7 @@ //! //! Cai, L. (2010). A two-tier full-information item factor analysis model //! with applications. *Psychometrika, 75*(4), 581-612. -//! https://doi.org/10.1007/s11336-010-9178-0 (abstract + metadata read; -//! full text not accessible — no equation locator drawn from it) +//! https://doi.org/10.1007/s11336-010-9178-0 (full text read, pp. 583-584) //! //! Cai, L., Yang, J. S., & Hansen, M. (2011). Generalized full-information //! item bifactor analysis. *Psychological Methods, 16*(3), 221-248. @@ -78,6 +76,8 @@ use crate::two_tier_grm::{ /// nothing is clamped or defaulted (ADR-0028 / #1929). #[derive(Clone, Copy, Debug)] pub struct TwoTierOakesConfig { + /// Estimate primary correlations; false excludes Phi from the free vector. + pub estimate_primary_correlation: bool, /// Gauss–Hermite nodes per primary dimension (any `q >= 1`). pub q_primary: usize, /// Gauss–Hermite nodes per specific factor (any `q >= 1`). @@ -156,6 +156,7 @@ impl Provider { return Err(format!("q_specific must be >= 1; got {}", cfg.q_specific)); } let iter_cfg = TwoTierGrmConfig { + estimate_primary_correlation: cfg.estimate_primary_correlation, q_primary: cfg.q_primary, q_specific: cfg.q_specific, max_iter: 1, @@ -212,11 +213,18 @@ impl Provider { slots, }); } - let n_phi = n_primary * (n_primary.saturating_sub(1)) / 2; + let n_phi = if cfg.estimate_primary_correlation { + n_primary * (n_primary.saturating_sub(1)) / 2 + } else { + 0 + }; let mut phi_slots = Vec::with_capacity(n_phi); let mut t = 0usize; for i in 0..n_primary { for j in (i + 1)..n_primary { + if !cfg.estimate_primary_correlation { + continue; + } phi_slots.push(cursor); labels.push(format!("phi_z:{i}:{j}")); cursor += 1; @@ -255,6 +263,9 @@ impl Provider { let p = self.v.n_primary; let params = pack_params(&self.v, a_primary, a_specific, thresholds); let z = z_from_phi(phi, p)?; + if self.n_phi == 0 && phi != phi_from_z(&vec![0.0; p * (p - 1) / 2], p) { + return Err("fixed primary correlation requires phi = I".into()); + } let mut out = vec![0.0f64; self.free_len()]; for (i, spec) in self.specs.iter().enumerate() { let mut s = 0usize; @@ -297,7 +308,11 @@ impl Provider { let _ = i; items.push(ItemParams { a_p, a_s, d }); } - let z: Vec = self.phi_slots.iter().map(|&s| packed[s]).collect(); + let z: Vec = if self.n_phi == 0 { + vec![0.0; p * (p - 1) / 2] + } else { + self.phi_slots.iter().map(|&s| packed[s]).collect() + }; (items, z) } @@ -588,6 +603,8 @@ fn cholesky_inverse(info: &[f64], k: usize) -> Result, String> { /// Observed-information standard errors via Oakes (1999, eq. 6, p. 480) at /// GIVEN two-tier item and primary-correlation parameters. +/// With fixed Phi=I, only item parameters enter the free vector (Cai, 2010, +/// pp. 583-584; Oakes, 1999, eq. 6, p. 480). #[allow(clippy::too_many_arguments)] pub fn two_tier_oakes_se( a_primary: &[f64], diff --git a/crates/mlsirm-core/tests/two_tier_grm_mirt_agreement.rs b/crates/mlsirm-core/tests/two_tier_grm_mirt_agreement.rs index 060cfb42c..df6335295 100644 --- a/crates/mlsirm-core/tests/two_tier_grm_mirt_agreement.rs +++ b/crates/mlsirm-core/tests/two_tier_grm_mirt_agreement.rs @@ -215,6 +215,7 @@ fn two_tier_grm_agrees_with_mirt_bfactor_two_tier_graded() { assert!(n_persons > 100, "fixture dataset must be non-trivial"); let cfg = TwoTierGrmConfig { + estimate_primary_correlation: true, q_primary: QUADPTS, q_specific: QUADPTS, max_iter: 2000, diff --git a/crates/mlsirm-core/tests/two_tier_grm_node_agreement.rs b/crates/mlsirm-core/tests/two_tier_grm_node_agreement.rs index ca674ab1e..ecc8f0e9e 100644 --- a/crates/mlsirm-core/tests/two_tier_grm_node_agreement.rs +++ b/crates/mlsirm-core/tests/two_tier_grm_node_agreement.rs @@ -16,27 +16,24 @@ //! q_specific`), which still exercises the SAME arbitrary-`n` Gauss-Hermite //! path (`quadrature::gh_rule`) the production `n_primary = 2` config in //! `two_tier_grm_recovery.rs` uses; node-count convergence is a property of -//! the quadrature rule, not of `n_primary`. `n_persons = 300` keeps +//! the quadrature rule, not of `n_primary`. `n_persons = 60` keeps //! `--ignored` runtime bounded; node count, not sample size or dimension -//! count, is what this test exercises. +//! count, is what this test exercises. The identity-mode fit and Oakes SE +//! both use 121 and 241 nodes per dimension. //! Run with `cargo test --release -- --ignored --nocapture`. use mlsirm_core::two_tier_grm::{fit_two_tier_grm, TwoTierGrmConfig}; +use mlsirm_core::two_tier_oakes::{two_tier_oakes_se, TwoTierOakesConfig}; const N_ITEMS: usize = 4; const N_PRIMARY: usize = 1; const N_SPECIFIC: usize = 1; const N_CAT: usize = 3; const PRIMARY_MAP: [bool; N_ITEMS * N_PRIMARY] = [true, true, true, true]; -const SPECIFIC_MAP: [i32; N_ITEMS] = [0, 0, 0, 0]; +const SPECIFIC_MAP: [i32; N_ITEMS] = [-1, -1, 0, 0]; const TRUE_A_PRIMARY: [f64; N_ITEMS] = [1.20, 1.00, 0.90, 1.10]; -const TRUE_A_SPECIFIC: [f64; N_ITEMS] = [0.80, 0.70, 0.90, 0.60]; -const TRUE_D: [[f64; 2]; N_ITEMS] = [ - [1.00, -1.00], - [0.80, -1.20], - [1.10, -0.90], - [0.90, -1.10], -]; +const TRUE_A_SPECIFIC: [f64; N_ITEMS] = [0.0, 0.0, 0.90, 0.60]; +const TRUE_D: [[f64; 2]; N_ITEMS] = [[1.00, -1.00], [0.80, -1.20], [1.10, -0.90], [0.90, -1.10]]; struct Lcg(u64); @@ -87,10 +84,11 @@ fn simulate(n_persons: usize, seed: u64) -> Vec { fn fit_config() -> TwoTierGrmConfig { TwoTierGrmConfig { + estimate_primary_correlation: false, q_primary: 15, q_specific: 15, - max_iter: 1000, - tol: 1e-5, + max_iter: 500, + tol: 1e-4, n_starts: 1, seed: 0x9E37_79B9_7F4A_7C15, newton_iter: 10, @@ -101,13 +99,14 @@ fn fit_config() -> TwoTierGrmConfig { #[test] #[ignore = "slow (121/241 node grids); run with: cargo test --release -- --ignored --nocapture"] fn two_tier_grm_121_vs_241_nodes_agree() { - let n_persons = 300; + let n_persons = 60; let y = simulate(n_persons, 20_260_917); let mut loglik = [0.0f64; 2]; let mut a_primary = [Vec::new(), Vec::new()]; let mut a_specific = [Vec::new(), Vec::new()]; let mut threshold = [Vec::new(), Vec::new()]; + let mut standard_errors = [Vec::new(), Vec::new()]; let mut elapsed_secs = [0.0f64; 2]; for (idx, &q) in [121usize, 241usize].iter().enumerate() { @@ -136,10 +135,38 @@ fn two_tier_grm_121_vs_241_nodes_agree() { "q={q} fit must converge (termination: {})", fit.termination_reason ); + assert_eq!(fit.phi, vec![1.0]); + assert_eq!(fit.primary_identification, "orthogonal"); loglik[idx] = *fit.loglik_trace.last().expect("non-empty trace"); a_primary[idx] = fit.a_primary.clone(); a_specific[idx] = fit.a_specific.clone(); threshold[idx] = fit.threshold.clone(); + let se = two_tier_oakes_se( + &fit.a_primary, + &fit.a_specific, + &fit.threshold, + &fit.phi, + &y, + None, + &PRIMARY_MAP, + &SPECIFIC_MAP, + n_persons, + N_ITEMS, + N_PRIMARY, + N_SPECIFIC, + N_CAT, + &TwoTierOakesConfig { + estimate_primary_correlation: false, + q_primary: q, + q_specific: q, + fd_step: 1e-5, + }, + ) + .expect("identity Oakes information must compute"); + assert!(se.labels.iter().all(|label| !label.starts_with("phi_z:"))); + standard_errors[idx] = se + .se + .expect("identity Oakes information must be positive definite"); eprintln!( "q={q}: n_iter={}, final_loglik={:.6}, elapsed={:.2}s", fit.n_iter, loglik[idx], elapsed_secs[idx] @@ -153,7 +180,7 @@ fn two_tier_grm_121_vs_241_nodes_agree() { elapsed_secs[0], elapsed_secs[1] ); // Same tolerance basis as the bifactor 121-vs-241 regression: the EM - // stopping tolerance is 1e-5, so two well-converged fits at different + // stopping tolerance is 1e-4, so two well-converged fits at different // (already-stabilized) node counts should agree to a small multiple of // that, not to float epsilon (independent EM runs land at slightly // different points on a flat likelihood ridge). @@ -183,4 +210,10 @@ fn two_tier_grm_121_vs_241_nodes_agree() { threshold[1][i] ); } + for i in 0..standard_errors[0].len() { + assert!( + (standard_errors[0][i] - standard_errors[1][i]).abs() < 5e-3, + "SE[{i}] disagrees between 121 and 241 nodes" + ); + } } diff --git a/crates/mlsirm-core/tests/two_tier_grm_recovery.rs b/crates/mlsirm-core/tests/two_tier_grm_recovery.rs index 1e1a1c735..f694ea9d5 100644 --- a/crates/mlsirm-core/tests/two_tier_grm_recovery.rs +++ b/crates/mlsirm-core/tests/two_tier_grm_recovery.rs @@ -177,6 +177,7 @@ fn simulate(n_persons: usize, seed: u64) -> (Vec, Vec) { fn fit_config() -> TwoTierGrmConfig { TwoTierGrmConfig { + estimate_primary_correlation: true, q_primary: 15, q_specific: 11, max_iter: 1000, diff --git a/crates/mlsirm-core/tests/two_tier_oakes_mirt.rs b/crates/mlsirm-core/tests/two_tier_oakes_mirt.rs index 99626e09a..082197215 100644 --- a/crates/mlsirm-core/tests/two_tier_oakes_mirt.rs +++ b/crates/mlsirm-core/tests/two_tier_oakes_mirt.rs @@ -167,6 +167,7 @@ fn rust_two_tier_oakes_se_matches_mirt_fixture() { let pmap = primary_map(); let smap = SPECIFIC_MAP.to_vec(); let cfg = TwoTierOakesConfig { + estimate_primary_correlation: true, q_primary: QUADPTS, q_specific: QUADPTS, fd_step: 1e-5, diff --git a/crates/mlsirm-core/tests/two_tier_reduces_to_bifactor.rs b/crates/mlsirm-core/tests/two_tier_reduces_to_bifactor.rs index 3b9b763da..7e03134ab 100644 --- a/crates/mlsirm-core/tests/two_tier_reduces_to_bifactor.rs +++ b/crates/mlsirm-core/tests/two_tier_reduces_to_bifactor.rs @@ -145,6 +145,7 @@ fn two_tier_with_single_primary_matches_bifactor_fit() { 2, N_CAT, &TwoTierGrmConfig { + estimate_primary_correlation: false, q_primary: 15, q_specific: 11, max_iter: 1000, @@ -161,6 +162,7 @@ fn two_tier_with_single_primary_matches_bifactor_fit() { "two-tier P=1 fit must converge ({}); the P=1 model IS the bifactor model", two_tier.termination_reason ); + assert_eq!(two_tier.primary_identification, "orthogonal"); let mut worst_slope = 0.0f64; for i in 0..N_ITEMS { diff --git a/python/fast_mlsirm/two_tier_grm.py b/python/fast_mlsirm/two_tier_grm.py index 003c007a9..dd45dd4c3 100644 --- a/python/fast_mlsirm/two_tier_grm.py +++ b/python/fast_mlsirm/two_tier_grm.py @@ -2,7 +2,7 @@ Cai, Yang, & Hansen, 2011; Gibbons et al., 2007). Each item's ordered categories load a caller-supplied subset of the -correlated primary dimensions plus at most one orthogonal specific factor. +primary dimensions plus at most one orthogonal specific factor. Estimation is Bock-Aitkin marginal maximum likelihood (Cai et al., 2011, "Maximum Marginal Likelihood Estimation" section) with dimension reduction over the specific tier; the numerical work runs in Rust @@ -19,11 +19,12 @@ comparison. - Two-tier latent covariance ``Sigma = [[G, 0], [0, diag(S)]]``: primaries ``theta_P ~ MVN(0, Phi)`` with ``Phi`` a correlation matrix (unit - diagonal, free off-diagonals — the single-group identification), specifics + diagonal, optionally fixed to identity), specifics orthogonal ``N(0, 1)`` (Chalmers, 2026, mirt ``bfactor`` documentation, "Details" section, which cites Cai, 2010). The bifactor model is the special case of one primary dimension (same source). The two-tier model - itself is Cai (2010) (abstract read; full text not accessible). + itself is Cai (2010, pp. 583-584, full text read). Fixing ``Phi = I`` is + an orthogonal-primary restriction of its covariance structure. - Caller-supplied confirmatory primary pattern (fixed zeros are never estimated); rotation with correlated primaries is the caller's identification responsibility (implementation scope choice; Cai, 2010, is a @@ -71,8 +72,7 @@ Cai, L. (2010). A two-tier full-information item factor analysis model with applications. *Psychometrika, 75*(4), 581-612. - https://doi.org/10.1007/s11336-010-9178-0 (abstract read; full text - not accessible — no equation locator is drawn from it) + https://doi.org/10.1007/s11336-010-9178-0 (full text read) Cai, L., Yang, J. S., & Hansen, M. (2011). Generalized full-information item bifactor analysis. *Psychological Methods, 16*(3), 221-248. @@ -151,12 +151,14 @@ class TwoTierGrmFit: ``a_specific`` the ``n_items`` specific slopes (``0`` for specific-free items, canonicalized within each block); ``threshold`` the ``n_items x (n_cat-1)`` strictly decreasing boundary intercepts; ``phi`` - the ``n_primary x n_primary`` estimated primary correlation matrix (unit + the ``n_primary x n_primary`` primary correlation matrix (unit diagonal); ``theta_p_eap`` / ``theta_p_sd`` the primary-factor EAPs and marginal posterior SDs (``n_persons x n_primary``); ``category_counts`` the observed ``n_items x n_cat`` counts. ``termination_reason`` is ``"tolerance_met"`` or ``"max_iter_reached"``; ``best_start`` the winning - start in ``0..n_starts``. + start in ``0..n_starts``. ``primary_identification`` is + ``"orthogonal"`` when Phi was fixed to I or ``"correlated"`` when + its off-diagonal entries were estimated (Cai, 2010, pp. 583-584). """ a_primary: np.ndarray @@ -176,6 +178,7 @@ class TwoTierGrmFit: final_loglik_change: float best_start: int n_parameters: int + primary_identification: str def fit_two_tier_grm( @@ -191,6 +194,7 @@ def fit_two_tier_grm( tol: float, n_starts: int, seed: int, + primary_correlation: str = "estimate", ) -> TwoTierGrmFit: """Fit the single-group polytomous two-tier GRM (compute in Rust). @@ -216,7 +220,22 @@ def fit_two_tier_grm( See the module docstring for the model, the paper basis of every non-obvious decision, and the APA 7th references. + ``primary_correlation='estimate'`` preserves the existing correlated-primary + fit; ``'identity'`` fixes Phi exactly to I (Cai, 2010, pp. 583-584). + In identity mode, distinct free-loading item sets for each primary pair + are necessary to rule out continuous orthogonal rotations; per-column + reflection canonicalization handles the remaining sign ambiguity. + + References (APA 7th ed.): Cai, L. (2010). A two-tier full-information item + factor analysis model with applications. *Psychometrika, 75*(4), 581-612. + https://doi.org/10.1007/s11336-010-9178-0; Cai, L., Yang, J. S., & + Hansen, M. (2011). Generalized full-information item bifactor analysis. + *Psychological Methods, 16*(3), 221-248. https://doi.org/10.1037/a0023350. """ + if not isinstance(primary_correlation, str) or primary_correlation not in ( + "estimate", "identity" + ): + raise ValueError("primary_correlation must be 'estimate' or 'identity'") n_cat_int = _finite_integer_control(n_cat, "n_cat") if n_cat_int < 2: raise ValueError("n_cat must be >= 2") @@ -259,6 +278,15 @@ def fit_two_tier_grm( if pmap.ndim != 2 or pmap.shape != (n_items, n_primary_int): raise ValueError("primary_map must be an n_items x n_primary boolean array") pmap_bool = np.asarray(pmap, dtype=bool) + if primary_correlation == "identity": + for d in range(n_primary_int): + for other in range(d + 1, n_primary_int): + if np.array_equal(pmap_bool[:, d], pmap_bool[:, other]): + raise ValueError( + "identity primary correlation requires distinct free-loading item sets " + "for every pair of primary columns; identical free-loading item sets " + f"at columns {d} and {other} admit orthogonal rotation" + ) smap = np.asarray(specific_map) if smap.ndim != 1 or smap.shape[0] != n_items: @@ -309,6 +337,7 @@ def fit_two_tier_grm( float(tol_float), int(n_starts_int), int(seed_int), + primary_correlation, ) return TwoTierGrmFit( a_primary=np.asarray(res["a_primary"], dtype=np.float64).reshape( @@ -340,6 +369,7 @@ def fit_two_tier_grm( final_loglik_change=float(res["final_loglik_change"]), best_start=int(res["best_start"]), n_parameters=int(res["n_parameters"]), + primary_identification=str(res["primary_identification"]), ) @@ -382,6 +412,7 @@ def two_tier_oakes_se( q_primary: int, q_specific: int, fd_step: float, + primary_correlation: str = "estimate", ) -> TwoTierOakesSe: """Observed-information SEs via Oakes (1999, eq. 6, p. 480) at given two-tier parameters (valid at every point, not only the MLE). @@ -391,7 +422,20 @@ def two_tier_oakes_se( Implementation basis: Oakes (1999, eq. 6, p. 480); Cai et al. (2011, eq. 6, p. 227); Gibbons et al. (2007, eq. 15). + ``primary_correlation='identity'`` excludes fixed Phi from the information + matrix (Cai, 2010, pp. 583-584; Oakes, 1999, eq. 6, p. 480). + + References (APA 7th ed.): Cai, L. (2010). A two-tier full-information item + factor analysis model with applications. *Psychometrika, 75*(4), 581-612. + https://doi.org/10.1007/s11336-010-9178-0; Oakes, D. (1999). Direct + calculation of the information matrix via the EM algorithm. *Journal of + the Royal Statistical Society: Series B, 61*(2), 479-482. + https://doi.org/10.1111/1467-9868.00188. """ + if not isinstance(primary_correlation, str) or primary_correlation not in ( + "estimate", "identity" + ): + raise ValueError("primary_correlation must be 'estimate' or 'identity'") n_cat_int = _finite_integer_control(n_cat, "n_cat") if n_cat_int < 2: @@ -464,6 +508,7 @@ def two_tier_oakes_se( int(q_primary_int), int(q_specific_int), float(fd_float), + primary_correlation, ) labels = [str(v) for v in res["labels"]] k = len(labels) diff --git a/tests/test_two_tier_grm.py b/tests/test_two_tier_grm.py index 61cec1371..c7613d982 100644 --- a/tests/test_two_tier_grm.py +++ b/tests/test_two_tier_grm.py @@ -11,8 +11,7 @@ --------- Cai, L. (2010). A two-tier full-information item factor analysis model with applications. *Psychometrika, 75*(4), 581-612. -https://doi.org/10.1007/s11336-010-9178-0 (abstract read; full text not -accessible — no equation locator is drawn from it) +https://doi.org/10.1007/s11336-010-9178-0 (full text read, pp. 583-584) """ from __future__ import annotations @@ -20,7 +19,7 @@ import numpy as np import pytest -from fast_mlsirm.two_tier_grm import fit_two_tier_grm +from fast_mlsirm.two_tier_grm import fit_two_tier_grm, two_tier_oakes_se N_PERSONS = 300 N_ITEMS = 6 @@ -63,11 +62,12 @@ ) -def _simulate(seed: int) -> np.ndarray: +def _simulate(seed: int, rho: float = RHO) -> np.ndarray: rng = np.random.default_rng(seed) z0 = rng.normal(0.0, 1.0, N_PERSONS) z1 = rng.normal(0.0, 1.0, N_PERSONS) - theta = np.stack([z0, RHO * z0 + np.sqrt(1.0 - RHO**2) * z1], axis=1) + theta = np.stack([z0, rho * z0 + np.sqrt(1.0 - rho**2) * z1], axis=1) + theta_s = rng.normal(0.0, 1.0, (N_PERSONS, N_SPECIFIC)) y = np.zeros((N_PERSONS, N_ITEMS), dtype=np.int64) for i in range(N_ITEMS): @@ -84,6 +84,139 @@ def _simulate(seed: int) -> np.ndarray: y[:, i] = (draws[:, None] > np.cumsum(probs, axis=1)).sum(axis=1) return y +def test_primary_correlation_validation() -> None: + with pytest.raises(ValueError, match="primary_correlation"): + _fit(_simulate(SEED), primary_correlation="unknown") + + +def test_identity_rejects_identical_primary_support() -> None: + y = _simulate(SEED) + shared = np.ones_like(PRIMARY_MAP) + with pytest.raises(ValueError, match="identical free-loading item sets"): + fit_two_tier_grm( + y, shared, SPECIFIC_MAP, N_CAT, N_PRIMARY, N_SPECIFIC, + q_primary=7, q_specific=7, max_iter=500, tol=1e-5, + n_starts=1, seed=SEED, primary_correlation="identity", + ) + + +# Cross-build tolerance for the stored origin/main literals below. The fixture +# stops on `delta_loglik <= 1e-2 * (1 + |loglik|)` (slack ~3.8e-1 at +# loglik = -37.15), so the pinned iterate is a trajectory point whose low-order +# bits differ between compilers; the observed cross-build spread on this fixture +# was at the 1e-9 scale, while injected behavioural regressions moved these +# values by 3.98e-4 to 1.05e-1. 1e-6 sits between those measured scales. +GOLDEN_CROSS_BUILD_ATOL = 1e-6 +GOLDEN_N_ITER = 6 +GOLDEN_TERMINATION = "tolerance_met" + + +def test_estimate_path_matches_origin_main_golden() -> None: + """Golden from origin/main 99c228a8f50a in a separate clean s1 checkout. + + Built there with Python 3.12 and ``pip install -e .``; ran this test's + deterministic 12 x 4 data through ``fit_two_tier_grm`` with the explicit + arguments below, then ``two_tier_oakes_se`` at the fitted parameters. + Cai (2010, pp. 583-584) defines the two-tier covariance; this checks that + the existing estimated-Phi path stayed numerically unchanged. + + Tolerance contract. Two things are checked with different strictness: + + * Within one build, the default path and ``primary_correlation="estimate"`` + must agree bit-for-bit, and the iteration count, termination reason and + trace length must match the golden exactly. Those are integer/string or + same-arithmetic comparisons, so they carry no cross-platform slack. + * The stored floating-point literals are compared with + ``GOLDEN_CROSS_BUILD_ATOL``. They were produced by one build + (linux-x86_64, pre-AVX host) and reproduce there exactly; a second build + (macos-arm64) was reported to differ at the 1e-9 scale on this fixture, + although that run's assertion text was not preserved. The fixture stops + on the relative rule ``delta_loglik <= tol * (1 + |loglik|)``, which with + ``tol=1e-2`` is a slack of about 3.8e-1 around ``loglik = -37.15``: the + iterate this golden pins is a point on the EM trajectory, not a converged + optimum, so ordinary floating-point reassociation between builds moves it + far more than 1e-12. Any real change to the estimate path moves these + values by far more than the tolerance below - deliberately injected + regressions on this fixture (identity rewiring, one fewer Newton step, a + 10x larger ridge) moved them by 1.05e-1, 2.09e-2 and 3.98e-4 - so the + guard keeps its purpose without asserting bitwise equality across + compilers. Note the smallest of those is 3.98e-4: a 1e-2 tolerance would + have missed it, which is why the bound below sits at 1e-6. + + Reference: Cai, L. (2010). A two-tier full-information item factor analysis + model with applications. *Psychometrika, 75*(4), 581-612. + https://doi.org/10.1007/s11336-010-9178-0 + """ + y = np.array([[(p + i) % 3 for i in range(4)] for p in range(12)], dtype=np.int64) + pmap = np.array([[1, 0], [1, 0], [0, 1], [0, 1]], dtype=bool) + smap = np.zeros(4, dtype=np.int64) + kwargs = dict(q_primary=7, q_specific=7, max_iter=500, tol=1e-2, n_starts=1, seed=20260922) + golden_phi = np.array([[1.0, -0.10480732118999127], [-0.10480732118999127, 1.0]]) + golden_trace = np.array([ + -57.878882056754186, -51.76952034453358, -47.89637529712579, + -43.0974994038811, -38.606846492548215, -37.22539643782572, + -37.15081748073924, + ]) + golden_primary = np.array([ + [-0.14916623496724393, 0.0], [0.2869920793276147, 0.0], + [0.0, 0.2869920842842997], [0.0, -0.14916624583907015], + ]) + golden_specific = np.array([ + 22.307488834567913, -0.8566330375435262, + -0.8566330381941958, 22.307488835442, + ]) + golden_threshold = np.array([ + [11.64599181577608, -11.591569652560967], + [0.6385437446044752, -1.0542766103752086], + [1.05427661081288, -0.6385437452077196], + [11.5915696386629, -11.645991831810742], + ]) + golden_labels = [ + "a_primary:0:0", "a_specific:0", "d:0:0", "d:0:1", + "a_primary:1:0", "a_specific:1", "d:1:0", "d:1:1", + "a_primary:2:1", "a_specific:2", "d:2:0", "d:2:1", + "a_primary:3:1", "a_specific:3", "d:3:0", "d:3:1", "phi_z:0:1", + ] + fits = {} + for correlation in (None, "estimate"): + options = kwargs if correlation is None else {**kwargs, "primary_correlation": correlation} + fit = fit_two_tier_grm(y, pmap, smap, 3, 2, 1, **options) + fits[correlation] = fit + # Exact, platform-independent parts of the contract. + assert fit.n_iter == GOLDEN_N_ITER + assert fit.termination_reason == GOLDEN_TERMINATION + assert fit.loglik_trace.shape == golden_trace.shape + np.testing.assert_array_equal(np.diag(fit.phi), np.ones(2)) + np.testing.assert_array_equal(fit.phi, fit.phi.T) + # Stored literals: cross-build tolerance (see docstring). + for got, expected in ( + (fit.phi, golden_phi), (fit.loglik_trace, golden_trace), + (fit.a_primary, golden_primary), (fit.a_specific, golden_specific), + (fit.threshold, golden_threshold), + ): + np.testing.assert_allclose( + got, expected, rtol=0, atol=GOLDEN_CROSS_BUILD_ATOL + ) + se = two_tier_oakes_se( + fit.a_primary, fit.a_specific, fit.threshold, fit.phi, + y, pmap, smap, 3, 2, 1, q_primary=7, q_specific=7, fd_step=1e-5, + **({} if correlation is None else {"primary_correlation": correlation}), + ) + assert se.labels == golden_labels + assert se.information.shape == (17, 17) + + # Same build, same arithmetic: the default path and the explicit + # "estimate" path must be bit-identical, with no tolerance at all. + default_fit, explicit_fit = fits[None], fits["estimate"] + for left, right in ( + (default_fit.phi, explicit_fit.phi), + (default_fit.loglik_trace, explicit_fit.loglik_trace), + (default_fit.a_primary, explicit_fit.a_primary), + (default_fit.a_specific, explicit_fit.a_specific), + (default_fit.threshold, explicit_fit.threshold), + ): + np.testing.assert_array_equal(left, right) + def _fit(y: np.ndarray, **overrides): kwargs = { diff --git a/tests/unit/two_tier_grm_tests.rs b/tests/unit/two_tier_grm_tests.rs index 66434d5f9..8be1ab0b7 100644 --- a/tests/unit/two_tier_grm_tests.rs +++ b/tests/unit/two_tier_grm_tests.rs @@ -242,6 +242,7 @@ fn single_primary_oracle_matches_stage1_bifactor_oracle() { fn valid_config() -> TwoTierGrmConfig { TwoTierGrmConfig { + estimate_primary_correlation: true, q_primary: 7, q_specific: 7, max_iter: 5, @@ -426,6 +427,49 @@ fn rejects_zero_quadrature_counts() { } } +#[test] +fn identity_rejects_identical_primary_support_but_accepts_nested_support() { + let (y, n_persons) = tiny_data(); + let cfg = TwoTierGrmConfig { + estimate_primary_correlation: false, + max_iter: 500, + ..valid_config() + }; + let shared = [true; TINY_N_ITEMS * TINY_N_PRIMARY]; + let err = fit_two_tier_grm( + &y, + None, + &shared, + &TINY_SPECIFIC_MAP, + n_persons, + TINY_N_ITEMS, + TINY_N_PRIMARY, + TINY_N_SPECIFIC, + TINY_N_CAT, + &cfg, + ) + .expect_err("identical supports admit a continuous rotation"); + assert!(err.contains("identical free-loading item sets"), "{err}"); + + let nested: Vec = (0..TINY_N_ITEMS).flat_map(|i| [true, i < 5]).collect(); + let fit = fit_two_tier_grm( + &y, + None, + &nested, + &TINY_SPECIFIC_MAP, + n_persons, + TINY_N_ITEMS, + TINY_N_PRIMARY, + TINY_N_SPECIFIC, + TINY_N_CAT, + &cfg, + ) + .expect("distinct nested supports must fit"); + assert!(fit.converged, "nested-support fit did not converge"); + assert_ne!(fit.termination_reason, "max_iter_reached"); + assert_eq!(fit.phi, vec![1.0, 0.0, 0.0, 1.0]); +} + #[test] fn rejects_bad_iteration_and_tolerance_budgets() { let (y, n_persons) = tiny_data(); diff --git a/tests/unit/two_tier_oakes_tests.rs b/tests/unit/two_tier_oakes_tests.rs index 108b3e293..0d547f1fb 100644 --- a/tests/unit/two_tier_oakes_tests.rs +++ b/tests/unit/two_tier_oakes_tests.rs @@ -28,6 +28,7 @@ fn tiny_design() -> (Vec, Vec, Vec, usize, usize) { fn fd_step_must_be_positive() { let (y, pmap, smap, n_persons, n_items) = tiny_design(); let cfg = TwoTierOakesConfig { + estimate_primary_correlation: true, q_primary: 5, q_specific: 5, fd_step: 0.0, @@ -49,6 +50,7 @@ fn returns_finite_information_on_tiny_case() { // Need every category observed per item — expand categories carefully. // With n_cat=2 the design above observes both cats on each item. let cfg = TwoTierOakesConfig { + estimate_primary_correlation: true, q_primary: 7, q_specific: 7, fd_step: 1e-5, @@ -113,6 +115,7 @@ fn two_tier_q241_rss_probe() { } let specific_map = vec![0i32, 0, 0, 0, 1, 1, 1, 1]; let cfg = TwoTierGrmConfig { + estimate_primary_correlation: true, q_primary: q, q_specific: q, max_iter: 1,