diff --git a/.claude/board/EPIPHANIES.md b/.claude/board/EPIPHANIES.md index 4ffa6005..171b28bf 100644 --- a/.claude/board/EPIPHANIES.md +++ b/.claude/board/EPIPHANIES.md @@ -1,10 +1,41 @@ +## 2026-08-04 — E-THE-CROSS-IDENTITY-SUITE-INHERITED-MY-BLIND-SPOT-1 — five independent cross-checks all passed while ω was inflated 3× and R² was unit-dependent + +**Status:** FINDING (both defects reproduced, fixed, regression-tested). **Confidence:** High — measured before and after. Found by external review of merged code, not by the suite that shipped with it. + +**What the suite could not see.** D-KIA-C1b's headline discipline was that every new estimator is cross-checked against an *independently computed* quantity — φ vs `pearson`, R² vs `pearson²`, η² vs R² on a dummy, η² vs `t²/(t²+df)`, ω vs α. All five held. Two P1 defects were live the whole time: + +- **ω erased loading signs.** The triad identity returns `λ_i²`, so the positive root was taken for every item; ω uses `(Σλ)²`, where a negatively-keyed item must *subtract*. On a signed fixture: **0.75 reported against a true 0.25**. +- **R² was unit-dependent.** Normal equations on raw columns plus an *absolute* `1e-12` pivot cutoff — a statement about units, not rank. An exact linear fit at 1e-8 magnitude returned **`None`**; the identical relationship at unit scale returned **1.0**. + +**Why the cross-identities were blind, and this is the transferable part.** A cross-identity checks that two computations *agree*, and agreement is only informative over the region the fixtures span. **Every ω fixture I wrote had all-positive loadings; every R² fixture was O(1).** The identities were not weak — they were evaluated at a single point of a space with unexplored dimensions. `ω = α` under tau-equivalence is true whether or not signs are handled, because tau-equivalent items are all positively keyed. `R² = pearson²` is true at any scale the guard happens to admit. + +**So the rule is not "add more cross-checks."** It is: **name the parameter space of the estimator, then confirm each dimension is actually varied by some fixture.** For these two the dimensions were obvious in hindsight and absent in fact — *sign* of the loadings, *magnitude* of the predictors. A cross-identity suite inherits the blind spots of the fixtures it runs on, and inherits them silently, because every check reports green. + +**The concrete practice, cheap enough to be worth doing:** for any numeric estimator, before claiming coverage, list the transformations under which the result should be **invariant** (affine rescaling, global sign flip, permutation, unit change) and the ones under which it should **change**, and write one test for each. Those tests found both defects immediately once written, and one of them — `omega_is_invariant_to_a_global_sign_flip` — is now the thing that proves the *fix* introduced no new bias. + +**Second lesson, on the earlier self-congratulation.** #887's write-up said the doc examples caught two defects the unit tests missed, and drew the moral that readable examples sample the degenerate case. True — but incomplete, and I stopped there. The fuller statement is that **fixtures are a sample of the input space and every suite is only as good as that sample's coverage**; doc examples happened to cover *one* uncovered corner. Treating that as the lesson rather than as an instance is why the next two corners stayed uncovered until someone else looked. + +## 2026-08-04 — E-A-STATISTIC-THAT-RETURNS-ZERO-FOR-UNDEFINED-CANNOT-FAIL-VISIBLY-1 — consolidating a utility does nothing while 56 copies with a *weaker contract* survive + +**Status:** FINDING (measured workspace-wide at `a9f813c`). **Confidence:** High — a grep census, reproducible. Filed as `TD-STATS-DEGENERACY-CONTRACT-DIVERGENCE`; **not paid** in this pass. + +**The census.** Hand-rolled `pearson` / `spearman` / `cronbach_alpha` definitions outside `jc`: **47 returning bare `f64`, 9 returning bare `f32`, 3 returning `Option`** — and all three `Option`s are `jc::reliability`'s own. At least one copy (`perturbation-sim/src/stats.rs:11`) returns **`0.0`** for ragged input, for `n < 2`, *and* for zero variance. + +**Why `0.0` is the sharp end.** It is byte-identical to a real "perfectly uncorrelated" result. A caller cannot distinguish *"r = 0 because these variables are unrelated"* — a finding — from *"r = 0 because I handed you a constant vector"* — a bug in the caller. The undefined case is not merely unreported; it is **disguised as the most publishable value in the range**. Several copies also omit the non-finite guard, so a `NaN` input yields a finite-looking number rather than a rejection. + +**This is the falsifiability rule one level down.** That rule governs guards that cannot fire; this is a *measurement* that cannot fail. Both defects share the shape: **the artifact has no way to say "no".** A statistic returning `Option` can; a statistic returning `f64` has spent its entire output range on answers and kept none for "I don't know". + +**The second-order finding, which is the transferable one.** `reliability.rs`'s module header states that callers rolled their own "**until now**" — presenting consolidation as accomplished. Measurement says consolidation reached the *four new* callers and **zero of the 56 pre-existing definitions**. The doc was written from the intent of the change rather than from the state of the tree, and nothing forced a re-read. **Generalisation: a "we consolidated X" claim is an inventory claim, and an inventory claim decays silently** — it is true on the day it is written and drifts every day after, because new duplicates cost nothing to add. Such claims need either a census in the same breath or the word *aspiration*. Cf. `CLAUDE.md` § falsifiability, "a doc-comment claim is not a behaviour". + +**Consequence, deliberately narrow:** the debt is recorded with a risk-ordered paydown (silent-`0.0` copies first, missing-non-finite-guard second, plain `f64` third, `f32` hot paths last or never) and explicitly **not** auto-fixable — each migration is a decision about what that call site should do when the estimate is undefined, which is exactly the information the current code discards. + ## 2026-08-04 — E-EXACT-FIT-IS-WHERE-ABSOLUTE-ZERO-GUARDS-BREAK-1 — a misfit guard written as `x < 0.0` rejects the PERFECT case, and the doc example is what caught it **Status:** FINDING (two measured defects, both fixed in the same commit). **Confidence:** High — reproduced, fixed, and each now carries a regression test. Code: `crates/jc/src/stats.rs` (D-KIA-C1b). **The defect.** `omega_total` guards against a **Heywood case** — estimated common variance exceeding an item's total variance, `ψ_i = σ_ii − λ_i² < 0`, which means the single-factor model does not fit. Written as `if psi < 0.0 { return None }`, that guard **rejects the best-fitting input there is**: items that are exact multiples of one factor have `ψ = 0` in real arithmetic and land a few ulps *either side* of zero in f64. The perfect model was reported as misfit. Fix: a tolerance **relative to the item's own variance** (`1e-9 · max(|σ_ii|, 1)`), then clamp — only a violation larger than rounding is a real one. -**The generalisation, which is the reusable part.** A guard testing a quantity that is **exactly zero at the ideal** must never use an absolute `< 0` (or `== 0`) comparison. The failure is invisible on ordinary noisy inputs — where `ψ` is comfortably positive — and appears **only on the cleanest possible data**, which is exactly the case a doc example or a smoke test reaches for. Sibling instances to check for: any `variance < 0`, `residual < 0`, `determinant == 0`, or `SS_within <= 0` test in a fitting path. +**The generalisation, which is the reusable part** *(narrowed 2026-08-04 after external review — the first statement was too broad).* The rule is about **theoretically non-negative FITTED quantities**: a residual variance, an estimated error term, a sum of squared residuals. Such a quantity must distinguish **rounding-scale negativity** (accept, clamp) from a **material model violation** (reject) using a SCALE-AWARE tolerance — never a bare `< 0`. It is **not** a licence to stop testing zero: a zero determinant may genuinely mean singularity, zero within-group variance genuinely leaves `F` undefined, and zero variance genuinely leaves `t` undefined. Those are real degeneracies to reject, not rounding to absorb. The earlier phrasing — *"a quantity exactly zero at the ideal must never use an absolute comparison"*, with determinants and `SS_within` listed as siblings — invited exactly that misreading. The failure is invisible on ordinary noisy inputs — where `ψ` is comfortably positive — and appears **only on the cleanest possible data**, which is exactly the case a doc example or a smoke test reaches for. Sibling instances to check for: any `variance < 0`, `residual < 0`, `determinant == 0`, or `SS_within <= 0` test in a fitting path. **The second defect, same commit, different lesson.** The `multiple_r_squared` doc example used predictors `a` and `a + 1` and asserted `R² = 1`. It returned `None` — **correctly**: the design carries an intercept, so `[1, a, a+1]` is rank-deficient. Here the *example* was wrong and the code was right. A duplicated column is the collinearity everyone remembers; **collinearity with the intercept is the one that gets written by accident**, because the two predictor columns are visibly different. diff --git a/.claude/board/STATUS_BOARD.md b/.claude/board/STATUS_BOARD.md index cc88a546..959cce29 100644 --- a/.claude/board/STATUS_BOARD.md +++ b/.claude/board/STATUS_BOARD.md @@ -6,12 +6,14 @@ Plan: `.claude/plans/kanban-64k-inverted-awareness-v1.md` (operator anchors a/b; |---|---|---|---|---| | D-KIA-0 | jc capability map + dichotomous-statistics decision note (phi/KR-20/kappa naming; Spearman dropped at view 2) | lance-graph | Queued | plan W0 | | D-KIA-A1 | ⊘ RESCOPED 2026-08-04 (E-ACTOR-IS-NOT-THE-PHASE-PATH-1): #879 is the complete phase-progression path; KanbanActor has no assigned architectural responsibility (legacy compatibility code). SHIPPED: held-owner reschedule/wake. OPEN: run_cycle drained-writer retry guard; missing-owner counter in cognitive_pass | lance-graph | Queued | plan W1 | -| D-KIA-C1b | jc additive-only extension: kappa + McDonald's omega + r-family effect size (R/R-squared, eta-squared = explained variance) + t-test (t/df/p) + a named phi wrapper. Cohen's d explicitly OUT — calculated separately if ever wanted. HARD CONSTRAINT: additive only — pearson/spearman/cronbach_alpha/icc keep their arithmetic, signature and semantics; any diff changing an existing jc statistic is an automatic reject. ONE sanctioned edit: widening reliability.rs private helpers (mean/all_finite/average_ranks/pop_var) to pub(crate) for reuse, visibility only, no body change. C1 audit found phi = pearson-on-binaries (already present in substance) and KR-20 = alpha-on-dichotomous (naming only); kappa absent = the real gap. SHIPPED as crates/jc/src/stats.rs: cohen_kappa, omega_total, phi, multiple_r/multiple_r_squared, eta_squared, t_test_one_sample/paired/welch/student, anova_one_way; 31 new tests (107 lib + 11 doctests green), clippy-clean. Existing-file diff is visibility-only (mean/all_finite -> pub(crate); average_ranks/pop_var NOT widened, unused). Unblocks D3's fusion falsifier | lance-graph | In PR | plan W0/C1b | +| D-KIA-C1b | jc additive-only extension: kappa + McDonald's omega + r-family effect size (R/R-squared, eta-squared = explained variance) + t-test (t/df/p) + a named phi wrapper. Cohen's d explicitly OUT — calculated separately if ever wanted. HARD CONSTRAINT: additive only — pearson/spearman/cronbach_alpha/icc keep their arithmetic, signature and semantics; any diff changing an existing jc statistic is an automatic reject. ONE sanctioned edit: widening reliability.rs private helpers (mean/all_finite/average_ranks/pop_var) to pub(crate) for reuse, visibility only, no body change. C1 audit found phi = pearson-on-binaries (already present in substance) and KR-20 = alpha-on-dichotomous (naming only); kappa absent = the real gap. SHIPPED as crates/jc/src/stats.rs: cohen_kappa, omega_total, phi, multiple_r/multiple_r_squared, eta_squared, t_test_one_sample/paired/welch/student, anova_one_way; 31 new tests (107 lib + 11 doctests green), clippy-clean. Existing-file diff is visibility-only (mean/all_finite -> pub(crate); average_ranks/pop_var NOT widened, unused). Unblocks D3a (overlap MEASUREMENT) — NOT a fusion claim: kappa is chance-corrected agreement under the observed marginals and says nothing about incremental value, so fusion still needs D3b's external criterion per the plan's own C3. Corrective slice (external review): omega sign-erasure + R-squared scale-dependence fixed; BinaryAssociation/kr20 added | lance-graph | Shipped (#887) + corrective slice | plan W0/C1b | | D-KIA-A2 | parallelism falsifier (protocol pre-registered: median-of-5, >=2x at >=4k owners, +/-10% stay-silent; kill = regrade claim (a)) | lance-graph | Queued | plan W2 | | D-KIA-B1 | catalog binary-range criterion contract type + generalized catalog-mirror drift guard | lance-graph | Queued | plan W3 | | D-KIA-C5 | cohort-statistic witness type under the ELEVATED carve-out + held-out anti-circularity gate | lance-graph | Queued | plan W4 | | D-KIA-D1 | observer/observed as two Locus categories over one arena (cheapest-first) | lance-graph | Queued | plan W5 | -| D-KIA-D3 | Horizontverschmelzung fusion falsifier, middle band pre-registered; corpus-side Synthesis producer (un-blocks session task #65 gate 1 — session-local task list, not a GitHub number) | lance-graph | Queued | plan W6 | +| D-KIA-C2 | Name the dichotomous statistics correctly (Pearson->phi, alpha->KR-20, kappa NOT a renamed ICC, Spearman dropped on binaries). AUDIT RESULT 2026-08-04: the jc reliability battery has exactly 4 consumers (style_table_agreement, rung_divergence_reliability, partof_isa_vs_palette256, l9_loci_real_text) and NONE is dichotomous — style columns, rung levels 1-10, palette/taxonomy distances, i4 loci offsets are all continuous/ordinal, so Pearson/alpha/ICC are correctly named at every existing call site and there is ZERO rename work today. The discipline binds PROSPECTIVELY at the first binary-criteria witness (D3). Surfaced instead: TD-STATS-DEGENERACY-CONTRACT-DIVERGENCE | lance-graph | Audited (no rename work; binds at D3) | plan W0/C2 | +| D-KIA-D3a | DESCRIPTIVE binary overlap: contingency counts + BOTH marginals + observed/expected agreement + kappa + phi, via jc::stats::binary_association. Claim ceiling is overlap / disagreement / marginal asymmetry / redundancy-or-complementarity CANDIDATE. No fusion or validity claim | lance-graph | Queued (unblocked) | plan W6 | +| D-KIA-D3b | HELD-OUT fusion falsifier — BLOCKED until an external criterion and a criterion-appropriate scoring rule are chosen. Continuous criterion: pre-registered delta-R-squared = R2(A+B) - max(R2(A),R2(B)). Binary criterion: a proper held-out score, NOT R-squared forced onto it. Reliability is not validity (plan C3) | lance-graph | Blocked (needs external criterion) | plan W6 | ## PROBE-BABEL-STANCES — two Rosetta stones + four-channel phase split (IN PR — slice 2, 2026-07-28) diff --git a/.claude/board/TECH_DEBT.md b/.claude/board/TECH_DEBT.md index 29d1b5d7..e5c88287 100644 --- a/.claude/board/TECH_DEBT.md +++ b/.claude/board/TECH_DEBT.md @@ -1,5 +1,50 @@ # Technical Debt Log — Open + Paid (double-entry, append-only) +## TD-STATS-DEGENERACY-CONTRACT-DIVERGENCE (2026-08-04) + +**Measured, not estimated** (at `a9f813c`, workspace-wide grep over `crates/`): + +| return type of hand-rolled `pearson` / `spearman` / `cronbach_alpha` | definitions | +|---|---| +| bare `f64` | 47 | +| bare `f32` | 9 | +| **`Option` (the degeneracy contract)** | **3 — all of them jc's own** | + +**The debt is a CONTRACT divergence, not duplication.** `jc::reliability`'s +whole documented design point is that degenerate input (constant series, +ragged lengths, `n < 2`, non-finite values, overflowed denominators) returns +`None`, so a caller can never mistake an undefined estimate for a measured +one. The 56 hand-rolled copies drop that guarantee, and at least one +**collapses every degeneracy to `0.0`** — `perturbation-sim/src/stats.rs:11` +returns `0.0` for ragged input, for `n < 2`, and for zero variance, which is +byte-identical to a genuine "perfectly uncorrelated" result. A caller cannot +distinguish *"r = 0 because the variables are unrelated"* from *"r = 0 because +I handed you a constant vector"*. Several also omit the non-finite guard +entirely, so a `NaN` in the input propagates to a finite-looking or `NaN` +result rather than being rejected. + +**This is the falsifiability rule's defect class one level down:** a statistic +that reports `0.0` for "undefined" **cannot fail visibly**, exactly like a +guard that never fires. It is worse than a wrong number because it is a wrong +number wearing the costume of a right one. + +**Also a correction to a shipped doc claim.** `reliability.rs`'s module header +says callers rolled their own "until now", implying consolidation happened. +The measurement says otherwise: consolidation happened for *new* callers only +(4 jc examples), while 56 pre-existing definitions remain. The header +overstates; treat it as *aspiration*, not inventory, until this is paid. + +**Not paid here, deliberately.** Migrating 56 definitions across ~26 files is +its own wave, each call site owned by its author per the workspace's +clippy-tier convention, and several are `f32` in hot paths where the `Option` +wrapper is a real (if small) signature change. **Suggested paydown order by +risk:** (1) any duplicate returning `0.0` for degeneracy — the silent-wrong +class; (2) duplicates without a non-finite guard; (3) plain `f64` duplicates +that already reject degeneracy some other way; (4) `f32` hot-path copies last, +if at all. **Do not auto-fix** — each migration is a decision about what the +call site should do when the estimate is undefined, which is precisely the +information the current code throws away. + ## TD-WORKSPACE-FMT-DRIFT (2026-07-30) **Measured, not estimated.** `cargo fmt --all -- --check` at diff --git a/.claude/plans/kanban-64k-inverted-awareness-v1.md b/.claude/plans/kanban-64k-inverted-awareness-v1.md index 05c532f1..636d4ce3 100644 --- a/.claude/plans/kanban-64k-inverted-awareness-v1.md +++ b/.claude/plans/kanban-64k-inverted-awareness-v1.md @@ -271,27 +271,105 @@ Blocks D3. > > **Two defects found by the doc examples, both fixed** — see > `E-EXACT-FIT-IS-WHERE-ABSOLUTE-ZERO-GUARDS-BREAK-1`. +> +> **⊘ CORRECTIVE SLICE (external review, 2026-08-04) — two P1 numerical +> defects in the shipped module, both reproduced before fixing:** +> +> 1. **ω erased loading SIGNS.** The triad identity yields `λ_i²`, so taking +> the positive root for every item dropped the sign of a negatively-keyed +> one — and ω depends on `(Σλ)²`, where a negative loading must *subtract*. +> Measured on an exact signed-congeneric fixture (λ = [+1,−1,+1]): reported +> **0.75** against a true **0.25**. Signs are now recovered from the +> covariance row (`sign(σ_ij) = s_i·s_j`, anchored `s_0 = +1`, which is free +> because ω is invariant to a global flip — asserted as a test). **Why the +> 31-test suite missed it: every ω fixture had all-positive loadings.** +> 2. **R² was scale-DEPENDENT.** The normal equations were built on raw columns +> and a pivot was called singular below an ABSOLUTE `1e-12` — a statement +> about units, not rank. Measured: an exact linear fit at 1e-8 magnitude +> returned `None` while the identical relationship at unit scale returned +> 1.0. The design is now centered and unit-normed (R² is affine-invariant, +> so no correct answer changes), the rank test is relative, and the `[0,1]` +> clamp is bounded to a rounding-scale band so a failed solve surfaces as +> `None` rather than as a plausible 0 or 1. +> +> **Also corrected:** ω's doc claimed rejections mean "the congeneric model +> does not fit" — it checks three *necessary* conditions (non-negative λ², no +> Heywood case, sign consistency) and **cannot certify a one-factor matrix**; +> the full vanishing-tetrad constraints are not tested and `k ≥ 4` misfit can +> still return a number. **Added:** `BinaryAssociation` + `binary_association` +> (counts + BOTH marginals + p_o/p_e alongside κ and φ — because the shipped +> φ doc said marginals are required to interpret it while the function returned +> one scalar), `kr20` (the dichotomous naming surface C2 asked for and C1b did +> not ship), and `betacf` now returns `None` on non-convergence instead of +> presenting the 300th iterate as a p-value. **116 lib + 13 doctests green.** **C2 — name the dichotomous statistics correctly.** Over binary catalog criteria: Pearson→**φ** (report the marginal-capped ceiling), Cronbach's α→**KR-20**; **κ is a SEPARATE estimator, not a renamed ICC** — where a continuous workflow would reach for ICC on binary criteria, compute **κ** instead, and keep **ICC as ICC** for the non-binary jc escalation only; -Spearman **degenerates and is dropped** at view 2 (it returns only in jc's -non-binary escalation). The implementation and every doc name the dichotomous +Spearman is **omitted as redundant** at view 2, not as degenerate: on two +non-constant binary variables the average-rank transform is affine, so ρ +carries exactly the same information as φ. (It returns in jc's non-binary +escalation, where the ranks stop being an affine image of the values.) The implementation and every doc name the dichotomous forms; reporting "Pearson" while computing φ is the defect class this arm exists to prevent. +> **✅ C2 AUDIT RESULT (2026-08-04) — no rename work exists today; the +> discipline binds prospectively at D3.** The `jc::reliability` battery has +> exactly **four** consumers, all in `crates/jc/examples/`, and **none feeds +> dichotomous data**: +> +> | consumer | what it correlates | domain | +> |---|---|---| +> | `style_table_agreement.rs` | style-table resonance / fan_out / exploration columns | continuous `f64` | +> | `rung_divergence_reliability.rs` | `dispatch_rung` vs `loci_rung` over the 34 recipes | ordinal 1–10 | +> | `partof_isa_vs_palette256.rs` | palette256² distances vs `part_of:is_a` path distances | continuous + multi-valued discrete (self-labelled in its own output) | +> | `l9_loci_real_text.rs` | i4 loci offsets (\|offset\|, 0..8) | multi-valued | +> +> So **Pearson / α / ICC are correctly named at every existing call site** — +> C2's defect ("reporting Pearson while computing φ") has **zero instances** +> to fix. Checked, not assumed: `probe_p1.rs` and `sigma_probe.rs` import +> other `jc` modules and never touch the battery, and +> `lance-graph/examples/g0_graph_loadbearing.rs` matched only a `println!` +> string naming `jc::reliability` as future work. Binary 0/1 codings DO exist +> in the tree (`probe_babel_stances.rs:784,792`) but feed a NARS belief arena, +> not a reliability statistic. +> +> **C2 is therefore a forward discipline, not a migration.** It binds at the +> first binary-criteria witness — which is D3's. `stats.rs` already makes half +> of it structural: `phi` takes `&[bool]`, so "compute φ, call it Pearson" +> cannot be written; the remaining half (feeding 0/1 `f64` to `pearson` and +> *reporting* "Pearson") stays a naming obligation on D3's author, with the +> marginal-capped ceiling caveat documented at `phi`. +> +> **What the audit surfaced instead — a real and larger finding:** +> `TD-STATS-DEGENERACY-CONTRACT-DIVERGENCE`. **56** hand-rolled +> `pearson`/`spearman`/`cronbach_alpha` definitions across ~26 files return a +> bare `f64`/`f32` against **3** returning `Option` (jc's own), and at least +> one collapses *every* degeneracy — ragged input, `n < 2`, zero variance — to +> **`0.0`**, indistinguishable from a genuine "uncorrelated" measurement. That +> is the falsifiability rule's defect class one level down: a statistic that +> cannot report its own undefinedness. Not paid here (it is its own wave); the +> risk-ordered paydown is in the TECH_DEBT entry. + **C3 — reliability vs validity split (hard gate).** α/KR-20/ICC/κ = **reliability**, claimable from the cohort alone. **Validity requires an external criterion** (an external gold-standard criterion, defined on the private consumer board) and is NOT claimed until one is wired. The plan's public claim ceiling until then: *"measurable reliability as a first step toward measurable awareness."* -**C4 — Jirak noise floors.** Binary criteria within one catalog panel are -domain-correlated — weak dependence *by construction*, so every -significance statement cites Jirak 2016 rates per `I-NOISE-FLOOR-JIRAK`; -classical IID Berry-Esseen is forbidden here exactly as for fingerprints. +**C4 — dependence-aware significance (scope corrected 2026-08-04).** Binary +criteria within one catalog panel are domain-correlated, so classical IID +significance is wrong here exactly as it is for fingerprints — that much +stands. **What does NOT follow is that `jirak.rs` is the answer.** The +p-values in `jc::stats` are classical independent-sample p-values and are +labelled as such at the module; `jirak.rs` is a fingerprint-specific empirical +probe, not a general uncertainty engine for κ / α / ω / ICC / R² / η² / φ / +t. A dependent-cohort significance claim therefore needs **its own justified +dependence model**, named at the claim site — citing `I-NOISE-FLOOR-JIRAK` is +a pointer to the *problem*, not a licence to reuse that implementation as the +*solution*. **C5 — witness storage under the ELEVATED carve-out.** The cohort statistic is a cross-input derivation of a different KIND than any observation → it may diff --git a/crates/jc/src/stats.rs b/crates/jc/src/stats.rs index 63e8b820..4087bc40 100644 --- a/crates/jc/src/stats.rs +++ b/crates/jc/src/stats.rs @@ -143,7 +143,7 @@ fn ln_gamma(x: f64) -> f64 { /// Continued fraction for the incomplete beta function (Numerical Recipes /// §6.4, modified Lentz). Converges for `x < (a+1)/(a+b+2)`. -fn betacf(a: f64, b: f64, x: f64) -> f64 { +fn betacf(a: f64, b: f64, x: f64) -> Option { const MAXIT: usize = 300; const EPS: f64 = 3.0e-16; const FPMIN: f64 = 1.0e-300; @@ -187,10 +187,12 @@ fn betacf(a: f64, b: f64, x: f64) -> f64 { let del = d * c; h *= del; if (del - 1.0).abs() < EPS { - break; + return h.is_finite().then_some(h); } } - h + // Exhausted MAXIT without meeting EPS. Returning the last iterate would + // present an unconverged value as a p-value; report failure instead. + None } /// Regularised incomplete beta `I_x(a, b)`, the CDF machinery behind every @@ -206,15 +208,23 @@ fn reg_inc_beta(a: f64, b: f64, x: f64) -> Option { if x == 1.0 { return Some(1.0); } - let front = - (ln_gamma(a + b) - ln_gamma(a) - ln_gamma(b) + a * x.ln() + b * (1.0 - x).ln()).exp(); + // `ln_1p(-x)` rather than `(1-x).ln()`: near x = 1 the subtraction loses + // most of its significant digits before the log ever sees it. + let front = (ln_gamma(a + b) - ln_gamma(a) - ln_gamma(b) + a * x.ln() + b * (-x).ln_1p()).exp(); let v = if x < (a + 1.0) / (a + b + 2.0) { - front * betacf(a, b, x) / a + front * betacf(a, b, x)? / a } else { // Symmetry: I_x(a,b) = 1 − I_{1−x}(b,a) - 1.0 - front * betacf(b, a, 1.0 - x) / b + 1.0 - front * betacf(b, a, 1.0 - x)? / b }; - v.is_finite().then(|| v.clamp(0.0, 1.0)) + // Clamp only a rounding-scale excursion; a materially out-of-domain result + // means the continued fraction failed and must surface as None, not as a + // plausible 0 or 1. + const SLACK: f64 = 1e-9; + if !v.is_finite() || !(-SLACK..=1.0 + SLACK).contains(&v) { + return None; + } + Some(v.clamp(0.0, 1.0)) } /// Two-tailed p-value of Student's `t` with `df` degrees of freedom: @@ -310,13 +320,35 @@ pub fn cohen_kappa(a: &[usize], b: &[usize]) -> Option { /// `λ_i² = σ_ij·σ_ik / σ_jk` for any pair `j,k ≠ i`; the estimate averages /// every admissible triad. Residual variances are `ψ_i = σ_ii − λ_i²`. /// +/// **Loading SIGNS matter and are recovered.** The triad identity gives only +/// `λ_i²`, hence `|λ_i|`; since ω depends on `(Σλ)²`, a negatively-keyed item +/// must subtract from that sum. Signs are read off the first row +/// (`sign(σ_ij) = s_i·s_j`, anchored at `s_0 = +1`, which is free because the +/// model is identified only up to a global flip and `(Σλ)²` is invariant to +/// one). Taking the positive root everywhere — as this function did before +/// 2026-08-04 — inflates ω badly; the regression fixture reports 0.25, where +/// the unsigned version returned 0.75. +/// +/// # What is and is NOT verified (read before quoting this as ω_t) +/// +/// This is a **triad-derived single-factor estimate**, not a fitted +/// factor-analysis solution. It checks three necessary conditions — +/// non-negative `λ_i²`, non-negative `ψ_i` (no Heywood case), and a +/// **sign-consistent** covariance pattern — and rejects with `None` when any +/// fails. It does **not** test the full vanishing-tetrad constraints, so it +/// cannot certify that a covariance matrix is one-factor: for `k ≥ 4` a +/// multidimensional or otherwise inconsistent matrix can still produce a +/// finite number. Treat a returned value as *ω under an assumed single factor*, +/// and establish the factor structure separately if the claim depends on it. +/// /// Requires `k ≥ 3` — with two items the single-factor model is **not /// identified** (one covariance, two unknown loadings), so `None` is returned /// rather than a fabricated estimate. Also returns `None` on ragged input, /// `n < 2`, non-finite input, a triad set with no usable denominator, a -/// negative `λ_i²`, or a **Heywood case** (`ψ_i < 0`, i.e. estimated common -/// variance exceeding the item's total variance) — each of which means the -/// congeneric model does not fit, not that reliability is low. +/// negative `λ_i²`, a **Heywood case** (`ψ_i < 0`, estimated common variance +/// exceeding the item's total variance), or an inconsistent sign pattern — +/// each of which means the assumed model is contradicted by the data, not that +/// reliability is low. /// /// ``` /// use jc::stats::omega_total; @@ -390,7 +422,51 @@ pub fn omega_total(items: &[Vec]) -> Option { if lam_sq < -tol || !lam_sq.is_finite() { return None; // single-factor model violated (negative common variance) } - lambda[i] = lam_sq.max(0.0).sqrt(); + lambda[i] = lam_sq.max(0.0).sqrt(); // MAGNITUDE only — signed below + } + + // ── recover the loading SIGNS (P1 fix, 2026-08-04) ────────────────────── + // + // The triad identity yields λ_i² and therefore only |λ_i|. Taking the + // positive root for every item is WRONG whenever the true loadings differ + // in sign: ω depends on `(Σλ)²`, and a negatively-keyed item must SUBTRACT + // from that sum. Erasing its sign inflates ω badly — on the signed fixture + // in the tests, 0.25 was reported as 0.75. + // + // Under a single factor `σ_ij = λ_iλ_j`, so `sign(σ_ij) = s_i·s_j`. The + // model is identified only up to a GLOBAL flip (and ω is invariant to one, + // since it uses `(Σλ)²`), so anchor `s_0 = +1` and read every other sign + // off the first row. + let mut sign = vec![1.0f64; k]; + for j in 1..k { + // A materially-zero covariance means λ_0 or λ_j is ~0; that item + // contributes ~nothing to Σλ, so either sign is harmless. Keep +1. + let scale = (cov[0][0].abs() * cov[j][j].abs()).sqrt().max(1.0); + if cov[0][j] < -1e-9 * scale { + sign[j] = -1.0; + } + } + // Sign CONSISTENCY is a genuine (partial) one-factor structure test: every + // off-diagonal covariance must agree with the product of the two recovered + // signs. A conflict means no single set of loadings reproduces the sign + // pattern, i.e. the data is not one-factor — reject rather than return a + // number the model cannot support. (This does NOT verify the full tetrad + // constraints; see the function's doc for what is and is not checked.) + for i in 0..k { + for j in (i + 1)..k { + let scale = (cov[i][i].abs() * cov[j][j].abs()).sqrt().max(1.0); + let cutoff = 1e-9 * scale; + if cov[i][j].abs() <= cutoff { + continue; // too near zero to carry a sign + } + let expected = sign[i] * sign[j]; + if cov[i][j].signum() != expected { + return None; // sign pattern is not one-factor + } + } + } + for i in 0..k { + lambda[i] *= sign[i]; } let sum_lambda: f64 = lambda.iter().sum(); @@ -444,6 +520,129 @@ pub fn phi(x: &[bool], y: &[bool]) -> Option { pearson(&xf, &yf) } +/// A 2×2 contingency table with **both marginals**, the agreement decomposition, +/// and the two association coefficients that read off it. +/// +/// This exists because **κ and φ are not interpretable as bare scalars.** φ's +/// attainable maximum is capped by how far the two marginals differ, so a +/// "low" φ may be sitting at its own ceiling; κ is a ratio whose denominator +/// `1 − p_e` collapses as the marginals become extreme, which makes κ unstable +/// exactly where the table is most lopsided. Reporting either number without +/// its marginals hides the information needed to read it — so the binary +/// surface returns the table, and the scalar functions remain conveniences. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct BinaryAssociation { + /// Count of `(false, false)`. + pub n00: u64, + /// Count of `(false, true)`. + pub n01: u64, + /// Count of `(true, false)`. + pub n10: u64, + /// Count of `(true, true)`. + pub n11: u64, + /// Rate of `true` in the first rater — the marginal φ's ceiling depends on. + pub positive_rate_a: f64, + /// Rate of `true` in the second rater. + pub positive_rate_b: f64, + /// `p_o` — proportion of cells where the two agree. + pub observed_agreement: f64, + /// `p_e` — agreement expected from the marginals alone. + pub expected_agreement: f64, + /// Cohen's κ, or `None` when `p_e == 1` (undefined, `0/0`). + pub kappa: Option, + /// φ, or `None` when either variable is constant (zero variance). + pub phi: Option, +} + +/// Cross-tabulate two binary vectors into a [`BinaryAssociation`] — the +/// **preferred** entry point for binary agreement, since it carries the +/// marginals κ and φ cannot be read without. +/// +/// Returns `None` only on structurally unusable input (length mismatch or +/// empty); a degenerate *table* still returns the counts, with `kappa` / `phi` +/// individually `None`, because the counts remain informative even where the +/// coefficients are undefined. +/// +/// ``` +/// use jc::stats::binary_association; +/// let a = [true, true, false, false]; +/// let b = [true, false, false, false]; +/// let t = binary_association(&a, &b).unwrap(); +/// assert_eq!((t.n11, t.n10, t.n01, t.n00), (1, 1, 0, 2)); +/// assert!((t.observed_agreement - 0.75).abs() < 1e-12); +/// ``` +pub fn binary_association(a: &[bool], b: &[bool]) -> Option { + if a.len() != b.len() || a.is_empty() { + return None; + } + let (mut n00, mut n01, mut n10, mut n11) = (0u64, 0u64, 0u64, 0u64); + for (&x, &y) in a.iter().zip(b.iter()) { + match (x, y) { + (false, false) => n00 += 1, + (false, true) => n01 += 1, + (true, false) => n10 += 1, + (true, true) => n11 += 1, + } + } + let n = a.len() as f64; + let pa = (n10 + n11) as f64 / n; + let pb = (n01 + n11) as f64 / n; + let p_o = (n00 + n11) as f64 / n; + // Expected agreement from the marginals alone. + let p_e = pa * pb + (1.0 - pa) * (1.0 - pb); + let kappa = { + let denom = 1.0 - p_e; + if denom == 0.0 || !denom.is_finite() { + None + } else { + let k = (p_o - p_e) / denom; + k.is_finite().then_some(k) + } + }; + Some(BinaryAssociation { + n00, + n01, + n10, + n11, + positive_rate_a: pa, + positive_rate_b: pb, + observed_agreement: p_o, + expected_agreement: p_e, + kappa, + phi: phi(a, b), + }) +} + +/// KR-20 — the Kuder-Richardson formula 20, which **is Cronbach's α computed +/// on dichotomous items**. +/// +/// Takes `&[Vec]` so the dichotomous precondition is enforced by the +/// type, converts once, and **delegates the arithmetic** to +/// [`crate::reliability::cronbach_alpha`] — this is a naming surface, not a +/// second implementation. Reporting "α" while the items are binary is the +/// mislabel this function exists to remove. +/// +/// Same shape and degeneracy conditions as `cronbach_alpha`: `k ≥ 2` items, +/// equal-length non-empty rows, and non-zero variance of the per-subject +/// totals (all-identical totals → `None`). +/// +/// ``` +/// use jc::stats::kr20; +/// let items = vec![ +/// vec![true, true, false, false], +/// vec![true, true, false, false], +/// vec![true, false, false, false], +/// ]; +/// assert!(kr20(&items).unwrap() > 0.8); +/// ``` +pub fn kr20(items: &[Vec]) -> Option { + let numeric: Vec> = items + .iter() + .map(|it| it.iter().map(|&v| if v { 1.0 } else { 0.0 }).collect()) + .collect(); + crate::reliability::cronbach_alpha(&numeric) +} + /// Solve `A·b = rhs` by Gaussian elimination with partial pivoting. /// `None` if the system is singular to working precision. fn solve(mut a: Vec>, mut rhs: Vec) -> Option> { @@ -459,7 +658,11 @@ fn solve(mut a: Vec>, mut rhs: Vec) -> Option> { .partial_cmp(&y.abs()) .unwrap_or(std::cmp::Ordering::Equal) })?; - if max.abs() < 1e-12 { + // RELATIVE rank test. Callers standardise their design so the initial + // diagonal is 1, making this a statement about rank rather than about + // the units the data happened to be measured in — an absolute cutoff + // here silently rejected full-rank designs at small magnitudes. + if max.abs() < 1e-10 { return None; // singular → collinear predictors } a.swap(col, piv); @@ -530,51 +733,82 @@ pub fn multiple_r_squared(y: &[f64], predictors: &[Vec]) -> Option { return None; } - // Design matrix with intercept: columns [1, p_0, .., p_{k-1}]. - let m = k + 1; - let col = |c: usize, i: usize| -> f64 { - if c == 0 { - 1.0 - } else { - predictors[c - 1][i] + // CENTER the response and every predictor, then SCALE each predictor to + // unit norm (P1 fix, 2026-08-04). + // + // R² is invariant under any affine transform of the predictors, so this + // changes no correct answer — but it is what makes the rank test + // meaningful. The previous code built `XᵀX` on the RAW columns and called + // a pivot singular below an ABSOLUTE `1e-12`; that threshold is a + // statement about units, not about rank, so merely measuring the same + // quantity in different units flipped a full-rank design to `None` + // (measured: `x` on the order of 1e-8 with an exact linear `y` returned + // `None`, while the identical relationship at unit scale returned 1.0). + // + // Centering also absorbs the intercept, so the design drops to `k` + // columns and the "collinear with the intercept" case (`b = a + 1`) + // correctly reappears as two identical centered columns. + let my = mean(y)?; + let yc: Vec = y.iter().map(|v| v - my).collect(); + let ss_tot: f64 = yc.iter().map(|v| v * v).sum(); + if ss_tot <= 0.0 || !ss_tot.is_finite() { + return None; // constant y → R² undefined + } + + let mut cols: Vec> = Vec::with_capacity(k); + for p in predictors { + let mp = mean(p)?; + let mut c: Vec = p.iter().map(|v| v - mp).collect(); + let norm = c.iter().map(|v| v * v).sum::().sqrt(); + if !norm.is_finite() || norm <= 0.0 { + return None; // constant predictor carries no rank } - }; + for v in &mut c { + *v /= norm; + } + cols.push(c); + } - // Normal equations XᵀX b = Xᵀy. - let mut xtx = vec![vec![0.0f64; m]; m]; - let mut xty = vec![0.0f64; m]; + // Normal equations on the standardised centered design (no intercept + // column). Every diagonal entry is now exactly 1, so a RELATIVE pivot + // threshold is well defined. + let mut xtx = vec![vec![0.0f64; k]; k]; + let mut xty = vec![0.0f64; k]; for (r, xtx_row) in xtx.iter_mut().enumerate() { for (c, cell) in xtx_row.iter_mut().enumerate() { - *cell = (0..n).map(|i| col(r, i) * col(c, i)).sum(); + *cell = (0..n).map(|i| cols[r][i] * cols[c][i]).sum(); } - xty[r] = (0..n).map(|i| col(r, i) * y[i]).sum(); + xty[r] = (0..n).map(|i| cols[r][i] * yc[i]).sum(); } if xtx.iter().any(|row| row.iter().any(|v| !v.is_finite())) || xty.iter().any(|v| !v.is_finite()) { - return None; // overflowed on large finite input + return None; } let beta = solve(xtx, xty)?; - let my = mean(y)?; - let ss_tot: f64 = y.iter().map(|&v| (v - my) * (v - my)).sum(); - if ss_tot == 0.0 || !ss_tot.is_finite() { - return None; // constant y → R² undefined - } let ss_res: f64 = (0..n) .map(|i| { - let pred: f64 = (0..m).map(|c| beta[c] * col(c, i)).sum(); - let e = y[i] - pred; + let pred: f64 = (0..k).map(|c| beta[c] * cols[c][i]).sum(); + let e = yc[i] - pred; e * e }) .sum(); - if !ss_res.is_finite() { + if !ss_res.is_finite() || ss_res < 0.0 { + return None; + } + let r2 = 1.0 - ss_res / ss_tot; + // `ss_res` is a sum of squares, so it cannot be negative; the only + // excursions possible are `ss_res` marginally exceeding `ss_tot` (R² + // slightly below 0) or cancellation pushing it a few ulps past 1. Clamp + // ONLY that rounding-scale band — a materially out-of-range value means + // the solve failed and must surface as `None`, not as a plausible 0 or 1. + const R2_SLACK: f64 = 1e-9; + if !(-R2_SLACK..=1.0 + R2_SLACK).contains(&r2) || !r2.is_finite() { return None; } - // Clamp: with an exact fit `ss_res` can land a few ulps below zero. - let r2 = (1.0 - ss_res / ss_tot).clamp(0.0, 1.0); - r2.is_finite().then_some(r2) + Some(r2.clamp(0.0, 1.0)) } /// Multiple correlation `R = √R²`. See [`multiple_r_squared`]. @@ -1030,6 +1264,127 @@ mod tests { assert_eq!(omega_total(&items), None); } + #[test] + fn omega_respects_loading_signs() { + // P1 REGRESSION (external review, 2026-08-04). The triad identity + // yields λ_i², hence only |λ_i|; taking the positive root for every + // item erased the sign of a negatively-keyed one. ω uses `(Σλ)²`, so a + // negative loading must SUBTRACT — erasing it inflated ω from 0.25 to + // 0.75 on exactly this fixture. + // + // Four mutually orthogonal mean-zero vectors give an exact one-factor + // construction with signed loadings λ = [+1, −1, +1] and equal + // orthogonal residuals: + // common variance = (1 − 1 + 1)²·V = V ; residual = 3V + // ω = V / (V + 3V) = 0.25 + let f = [1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0]; + let e1 = [1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0]; + let e2 = [1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0]; + let e3 = [1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0]; + let items = vec![ + (0..8).map(|i| f[i] + e1[i]).collect::>(), + (0..8).map(|i| -f[i] + e2[i]).collect::>(), + (0..8).map(|i| f[i] + e3[i]).collect::>(), + ]; + let w = omega_total(&items).unwrap(); + assert!( + approx(w, 0.25, 1e-12), + "signed-loading ω was {w}, expected 0.25" + ); + } + + #[test] + fn omega_is_invariant_to_a_global_sign_flip() { + // The one-factor model is identified only up to a GLOBAL flip, and ω + // uses (Σλ)², so flipping every item must leave ω unchanged. This also + // proves the sign anchor (s_0 = +1) introduces no arbitrary bias. + let items = congeneric_items(); + let flipped: Vec> = items + .iter() + .map(|it| it.iter().map(|v| -v).collect()) + .collect(); + let a = omega_total(&items).unwrap(); + let b = omega_total(&flipped).unwrap(); + assert!(approx(a, b, 1e-12), "global flip changed ω: {a} vs {b}"); + } + + #[test] + fn omega_rejects_an_inconsistent_sign_pattern() { + // Can-it-fire for the sign-consistency structure test, and written so + // it can ONLY pass if that specific guard is what rejects. + // + // At k = 3 the check is provably REDUNDANT: with one triad per item, + // sign(λ_0²) = sign(σ01)·sign(σ02)·sign(σ12), which is +1 for every + // consistent pattern and −1 for every inconsistent one — so an + // inconsistent triple always trips the negative-λ² guard first. Only + // at k ≥ 4, where λ² is AVERAGED over several triads, can the sign + // pattern conflict while every λ² stays non-negative. This fixture + // (found by search) is exactly that case. + let items = vec![ + vec![1.0, -2.0, 0.0, 3.0, 2.0, 2.0], + vec![2.0, 3.0, 1.0, 0.0, 0.0, 3.0], + vec![1.0, 2.0, -3.0, 3.0, -3.0, 0.0], + vec![2.0, -2.0, 1.0, -1.0, 3.0, -2.0], + ]; + let k = items.len(); + let cv = |x: &[f64], y: &[f64]| -> f64 { + let (mx, my) = (mean(x).unwrap(), mean(y).unwrap()); + x.iter() + .zip(y) + .map(|(a, b)| (a - mx) * (b - my)) + .sum::() + / (x.len() - 1) as f64 + }; + let c: Vec> = (0..k) + .map(|i| (0..k).map(|j| cv(&items[i], &items[j])).collect()) + .collect(); + + // PRE-REGISTER that neither pre-existing guard can be the rejecter: + // every λ² is non-negative and every ψ is non-negative. Without this, + // the test would pass even if the negative-λ² or Heywood guard fired, + // and would prove nothing about the sign check. + for i in 0..k { + let (mut acc, mut cnt) = (0.0, 0usize); + for j in 0..k { + if j == i { + continue; + } + for l in (j + 1)..k { + if l == i || c[j][l] == 0.0 { + continue; + } + acc += c[i][j] * c[i][l] / c[j][l]; + cnt += 1; + } + } + let lam_sq = acc / cnt as f64; + assert!( + lam_sq >= 0.0, + "item {i}: λ²={lam_sq} would trip the NEGATIVE guard" + ); + assert!( + c[i][i] - lam_sq >= 0.0, + "item {i}: ψ={} would trip the HEYWOOD guard", + c[i][i] - lam_sq + ); + } + // …and that the sign pattern really is inconsistent. + let mut sign = vec![1.0f64; k]; + for j in 1..k { + if c[0][j] < 0.0 { + sign[j] = -1.0; + } + } + let conflicts = (0..k) + .flat_map(|i| ((i + 1)..k).map(move |j| (i, j))) + .filter(|&(i, j)| c[i][j] != 0.0 && c[i][j].signum() != sign[i] * sign[j]) + .count(); + assert!(conflicts > 0, "fixture must carry a sign conflict"); + + // Therefore only the sign-consistency guard can produce the rejection. + assert_eq!(omega_total(&items), None); + } + #[test] fn omega_degenerate_returns_none() { // k < 3 → single-factor model unidentified. @@ -1089,6 +1444,103 @@ mod tests { assert_eq!(phi(&[true, false], &[true]), None); // ragged } + // ─────────────── binary association table + KR-20 ─────────────── + + #[test] + fn binary_association_agrees_with_the_scalar_functions() { + // Cross-identity: the table's κ and φ must equal the standalone + // functions computed independently, or the two surfaces have drifted. + let a = [ + true, true, false, true, false, false, true, false, true, false, + ]; + let b = [ + true, false, false, true, true, false, true, false, false, false, + ]; + let t = binary_association(&a, &b).unwrap(); + let ka: Vec = a.iter().map(|&v| v as usize).collect(); + let kb: Vec = b.iter().map(|&v| v as usize).collect(); + assert!(approx( + t.kappa.unwrap(), + cohen_kappa(&ka, &kb).unwrap(), + 1e-12 + )); + assert!(approx(t.phi.unwrap(), phi(&a, &b).unwrap(), 1e-12)); + // Counts must reconstruct the inputs exactly. + assert_eq!(t.n00 + t.n01 + t.n10 + t.n11, a.len() as u64); + assert!(approx( + t.observed_agreement, + (t.n00 + t.n11) as f64 / a.len() as f64, + 1e-12 + )); + } + + #[test] + fn binary_association_exposes_the_marginal_asymmetry_phi_alone_hides() { + // The reason this type exists. Two tables with a SIMILAR φ but very + // different marginals are not equally interpretable; the scalar cannot + // say so and the table can. Anti-vacuity: the marginals must actually + // differ materially between the two cases. + let bal_a = [true, true, false, false, true, false, true, false]; + let bal_b = [true, true, false, false, false, true, true, false]; + let skew_a = [true, true, true, true, true, true, true, false]; + let skew_b = [true, true, true, true, true, true, false, true]; + let bal = binary_association(&bal_a, &bal_b).unwrap(); + let skew = binary_association(&skew_a, &skew_b).unwrap(); + assert!( + (bal.positive_rate_a - 0.5).abs() < 1e-12 && skew.positive_rate_a > 0.85, + "fixtures must differ in marginals: {} vs {}", + bal.positive_rate_a, + skew.positive_rate_a + ); + // Expected agreement is what the marginals buy for free — far higher in + // the skewed table, which is exactly the caveat a bare φ omits. + assert!( + skew.expected_agreement > bal.expected_agreement + 0.2, + "skewed p_e {} should far exceed balanced p_e {}", + skew.expected_agreement, + bal.expected_agreement + ); + } + + #[test] + fn binary_association_degenerate_table_keeps_counts_but_drops_coefficients() { + // Constant inputs: φ undefined (zero variance) and κ undefined + // (p_e == 1), but the counts are still real information — the table + // must survive where the coefficients cannot. + let a = [true, true, true, true]; + let t = binary_association(&a, &a).unwrap(); + assert_eq!(t.n11, 4); + assert!(approx(t.observed_agreement, 1.0, 1e-12)); + assert_eq!(t.kappa, None); + assert_eq!(t.phi, None); + assert_eq!(binary_association(&[], &[]), None); + assert_eq!(binary_association(&[true], &[]), None); + } + + #[test] + fn kr20_is_cronbach_alpha_on_the_dichotomous_coding() { + // KR-20 IS α on binary items — assert the delegation rather than + // trusting it, computing the reference through the proven function. + let items = vec![ + vec![true, true, false, false, true], + vec![true, false, false, true, true], + vec![true, true, false, false, false], + ]; + let numeric: Vec> = items + .iter() + .map(|it| it.iter().map(|&v| if v { 1.0 } else { 0.0 }).collect()) + .collect(); + let expect = crate::reliability::cronbach_alpha(&numeric).unwrap(); + assert!(approx(kr20(&items).unwrap(), expect, 1e-15)); + } + + #[test] + fn kr20_degenerate_returns_none() { + assert_eq!(kr20(&[vec![true, false]]), None); // k < 2 + // Every subject total identical → no between-subject variance. + assert_eq!(kr20(&[vec![true, false], vec![false, true]]), None); + } + // ────────────────────────── R / R² ────────────────────────── #[test] @@ -1121,6 +1573,36 @@ mod tests { assert!(r2 > 0.05 && r2 < 0.99, "expected partial fit, got R²={r2}"); } + #[test] + fn r_squared_is_scale_and_offset_invariant() { + // P1 REGRESSION (external review, 2026-08-04). R² is invariant under + // any affine transform of the predictors, but the old code built the + // normal equations on RAW columns and rejected a pivot below an + // ABSOLUTE 1e-12 — a statement about UNITS, not rank. Measured: an + // exact linear fit at 1e-8 magnitude returned `None` while the + // identical relationship at unit scale returned 1.0. + let base: Vec = (1..=6).map(|k| k as f64).collect(); + let fit = |x: &Vec| -> Option { + let y: Vec = x.iter().map(|v| 3.0 + 2.0 * v).collect(); + multiple_r_squared(&y, std::slice::from_ref(x)) + }; + let unit = fit(&base).expect("unit scale must fit"); + assert!(approx(unit, 1.0, 1e-9)); + for factor in [1e-8, 1e-4, 1e4, 1e8] { + let scaled: Vec = base.iter().map(|v| v * factor).collect(); + let got = + fit(&scaled).unwrap_or_else(|| panic!("scale {factor} lost a full-rank design")); + assert!( + approx(got, unit, 1e-9), + "scale {factor}: R²={got} vs {unit}" + ); + } + // Large common offset — the case that stresses centering. + let shifted: Vec = base.iter().map(|v| v + 1e12).collect(); + let got = fit(&shifted).expect("offset lost a full-rank design"); + assert!(approx(got, unit, 1e-6), "offset: R²={got} vs {unit}"); + } + #[test] fn r_squared_never_decreases_when_a_predictor_is_added() { let y = vec![1.0, 3.0, 2.0, 6.0, 4.0, 9.0, 7.0, 11.0];