feat(two-tier): orthogonal primary identification with fixed Phi

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<Py<pyo3::types::PyDict>> {
+    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<Vec<bool>> = match &observed {
         Some(o) => Some(o.as_slice()?.to_vec()),
@@ -1913,6 +1922,7 @@ fn fit_two_tier_grm(
         })
         .collect::<PyResult<_>>()?;
     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<Py<pyo3::types::PyDict>> {
+    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<Vec<bool>> = 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,

crates/mlsirm-core/src/two_tier_grm.rs

@@ -195,8 +195,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 +245,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 +280,8 @@ pub struct TwoTierGrmResult {
     pub a_specific: Vec<f64>,
     /// Ordered boundary intercepts `d_ik`, row-major `n_items * (n_cat - 1)`.
     pub threshold: Vec<f64>,
-    /// 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<f64>,
     /// Primary-factor EAPs `E[theta_d | Y_p]`, row-major
     /// `n_persons * n_primary`.
@@ -297,9 +298,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
@@ -541,6 +545,7 @@ fn initial_params(
     observed: Option<&[bool]>,
     seed: u64,
     start: usize,
+    estimate_phi: bool,
 ) -> (Vec<ItemParams>, Vec<f64>) {
     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 +588,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();
         }
@@ -1330,7 +1335,14 @@ fn run_single_start(
     start: usize,
 ) -> Result<SingleStartOutcome, String> {
     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
@@ -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,

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<f64> = self.phi_slots.iter().map(|&s| packed[s]).collect();
+        let z: Vec<f64> = 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<Vec<f64>, 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],

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,

crates/mlsirm-core/tests/two_tier_grm_node_agreement.rs

@@ -87,6 +87,7 @@ fn simulate(n_persons: usize, seed: u64) -> Vec<usize> {
 
 fn fit_config() -> TwoTierGrmConfig {
     TwoTierGrmConfig {
+        estimate_primary_correlation: true,
         q_primary: 15,
         q_specific: 15,
         max_iter: 1000,