rpx_benchmark.phi

   1"""State-Transition Robustness Φ via repeated-measures MANOVA.
   2
   3Implements paper §3.2 from raw H (between-phase) and E (within-phase) SSCP
   4matrices. We do not depend on ``statsmodels``; only ``numpy`` (core) and
   5``scipy.stats.f`` (in the ``[analysis]`` extras) are required.
   6
   7What this module provides
   8-------------------------
   9
  10- :func:`standardise` — z-score each metric column independently across all
  11  ``3N`` (scene, phase) cells (paper §3.2).
  12- :func:`compute_phi_oneway` — one-way repeated-measures MANOVA on the phase
  13  factor. Returns Wilks-based and Pillai-based Φ plus an F-approximation
  14  p-value.
  15- :func:`compute_phi_per_transition` — pairwise Hotelling T² on within-scene
  16  paired differences for each of (Clu→Int, Int→Cln, Clu→Cln). Returns
  17  ``η² = T² / (T² + N − 1)`` and ``Φ = 1 − η²`` per transition.
  18- :func:`phi_interpretation` — Cohen anchor lookup (negligible / small /
  19  medium / large) for a given Φ value.
  20
  21Conventions
  22-----------
  23
  24All array shapes use the convention ``(N_scenes, n_phases, K_metrics)``.
  25Lower-is-better metrics should be **negated before standardisation** by
  26the caller (matching paper §3.2 — sign convention for the MANOVA, separate
  27from the direction-aware bounds policy used by JEDI).
  28"""
  29
  30from __future__ import annotations
  31
  32from dataclasses import dataclass, field
  33from typing import Any, Dict, List, Optional, Sequence, Tuple
  34
  35import numpy as np
  36
  37from .exceptions import MetricError
  38from .logging_utils import get_logger
  39
  40log = get_logger(__name__)
  41
  42
  43# --------------------------------------------------------------------------- #
  44# Standardisation
  45# --------------------------------------------------------------------------- #
  46
  47
  48def standardise(x: np.ndarray) -> np.ndarray:
  49    """Z-score each metric column across all (scene, phase) cells.
  50
  51    Parameters
  52    ----------
  53    x : ndarray, shape ``(N, P, K)``
  54        Raw metric values: ``N`` scenes, ``P`` phases, ``K`` metrics.
  55
  56    Returns
  57    -------
  58    ndarray, shape ``(N, P, K)``
  59        Each metric column normalised to zero mean / unit variance across
  60        the ``N*P`` cells. Constant columns are mapped to zeros (so they
  61        contribute nothing to the MANOVA test).
  62
  63    Notes
  64    -----
  65    Paper §3.2: standardisation is performed *across the 3N (scene, phase)
  66    pairs within a task*, so the K metrics enter the MANOVA on a common
  67    scale and the test is invariant to the units of the underlying raw
  68    scores.
  69    """
  70    if x.ndim != 3:
  71        raise MetricError(
  72            f"standardise expects shape (N, P, K), got {x.shape}",
  73        )
  74    N, P, K = x.shape
  75    flat = x.reshape(N * P, K).astype(np.float64)
  76    mu = flat.mean(axis=0)
  77    sigma = flat.std(axis=0, ddof=0)
  78    # Constant columns: variance 0; zero out instead of div-by-zero.
  79    safe = np.where(sigma > 0, sigma, 1.0)
  80    z = (flat - mu) / safe
  81    z[:, sigma == 0] = 0.0
  82    return z.reshape(N, P, K)
  83
  84
  85# --------------------------------------------------------------------------- #
  86# Collinearity diagnostics
  87# --------------------------------------------------------------------------- #
  88#
  89# Λ = |E| / |H + E| is numerically unstable when E is ill-conditioned, and
  90# ill-conditioned E almost always means two K-vector metrics are near-linear
  91# combinations of each other (paper §3.2 warns against this but does not
  92# gate the computation). We surface a warning at moderate condition numbers
  93# and refuse the run at severe ones — with an actionable hint pointing at
  94# the worst-correlated metric pair if names are available.
  95
  96#: Warn when cond(E) exceeds this — Λ may be unstable but the result is usable.
  97COLLINEARITY_WARN_COND: float = 30.0
  98
  99#: Refuse when cond(E) exceeds this — Λ is numerically unreliable.
 100COLLINEARITY_FAIL_COND: float = 1000.0
 101
 102
 103def _worst_correlated_pair(
 104    z: np.ndarray,
 105    metric_names: Optional[Sequence[str]],
 106) -> str:
 107    """Return a short human-readable hint identifying the most collinear
 108    pair of metrics in the standardised K-vector.
 109
 110    Falls back to an anonymous "column i vs column j" message when
 111    ``metric_names`` is None.
 112    """
 113    N, P, K = z.shape
 114    flat = z.reshape(N * P, K)
 115    if K < 2 or flat.shape[0] < 2:
 116        return ""
 117    with np.errstate(invalid="ignore"):
 118        C = np.corrcoef(flat, rowvar=False)
 119    # np.corrcoef returns NaN when a column is constant; treat those as 0.
 120    C = np.where(np.isfinite(C), C, 0.0)
 121    np.fill_diagonal(C, 0.0)
 122    i, j = np.unravel_index(int(np.argmax(np.abs(C))), C.shape)
 123    r = float(C[i, j])
 124    if metric_names is not None and len(metric_names) == K:
 125        return f"worst pair: ({metric_names[int(i)]!r}, {metric_names[int(j)]!r}) with |r|={abs(r):.3f}"
 126    return f"worst pair: (column {i}, column {j}) with |r|={abs(r):.3f}"
 127
 128
 129# --------------------------------------------------------------------------- #
 130# One-way repeated-measures MANOVA → Φ
 131# --------------------------------------------------------------------------- #
 132
 133
 134@dataclass
 135class PhiOneway:
 136    """Result of the one-way phase MANOVA.
 137
 138    Attributes
 139    ----------
 140    phi_wilks : float
 141        :math:`\\Phi = 1 - (1 - \\Lambda^{1/s})` from Wilks' Λ.
 142    phi_pillai : float
 143        :math:`\\Phi = 1 - V/s` from Pillai's trace V.
 144    p_value : float
 145        F-approximation p-value (Wilks → Rao's F). Treats per-test p as
 146        descriptive; multiple-comparison correction is the caller's
 147        responsibility (paper §3.2, Holm-Bonferroni within (model, task)).
 148    n_eff : int
 149        Number of scenes that contributed (complete-case; scenes with
 150        missing predictions in any phase are dropped upstream).
 151    K : int
 152        Number of metric dimensions.
 153    s : int
 154        Rank of the hypothesis matrix, ``min(K, df_effect)``.
 155    wilks_lambda : float
 156        The raw Λ statistic, in case callers want to recompute.
 157    pillai_V : float
 158        The raw Pillai trace, same reason.
 159    """
 160
 161    phi_wilks: float
 162    phi_pillai: float
 163    p_value: float
 164    n_eff: int
 165    K: int
 166    s: int
 167    wilks_lambda: float
 168    pillai_V: float
 169
 170    def to_dict(self) -> Dict[str, object]:
 171        return {
 172            "phi_wilks": self.phi_wilks,
 173            "phi_pillai": self.phi_pillai,
 174            "phi_conservative": self.phi_conservative,
 175            "p_value": self.p_value,
 176            "n_eff": self.n_eff,
 177            "K": self.K,
 178            "s": self.s,
 179            "wilks_lambda": self.wilks_lambda,
 180            "pillai_V": self.pillai_V,
 181        }
 182
 183    @property
 184    def phi_conservative(self) -> float:
 185        """The smaller (more conservative) of Wilks/Pillai Φ.
 186
 187        Per paper §3.2: *"main results report the more conservative
 188        estimate, with the full comparison in Appendix C."*
 189        """
 190        return min(self.phi_wilks, self.phi_pillai)
 191
 192
 193def _h_e_matrices_oneway(z: np.ndarray) -> Tuple[np.ndarray, np.ndarray, int]:
 194    """Build H (between-phase) and E (within-phase) SSCP matrices.
 195
 196    Repeated-measures structure: scenes are subjects; phase is the
 197    within-subjects factor with P levels.
 198
 199    Parameters
 200    ----------
 201    z : ndarray, shape ``(N, P, K)``
 202        Standardised metric values per (scene, phase).
 203
 204    Returns
 205    -------
 206    H : ndarray, shape ``(K, K)``
 207    E : ndarray, shape ``(K, K)``
 208    df_effect : int
 209        ``P - 1`` for the phase main effect.
 210    """
 211    N, P, K = z.shape
 212
 213    # Per-phase mean across scenes (shape (P, K)) and overall mean (shape (K,)).
 214    phase_means = z.mean(axis=0)
 215    grand_mean = z.reshape(N * P, K).mean(axis=0)
 216
 217    # H: between-phase SSCP, weighted by per-cell count N.
 218    diffs_h = phase_means - grand_mean  # (P, K)
 219    H = N * (diffs_h.T @ diffs_h)
 220
 221    # E: within-phase SSCP — sum over scenes & phases of (cell - phase_mean) outer.
 222    centered = z - phase_means[np.newaxis, :, :]  # (N, P, K)
 223    flat = centered.reshape(N * P, K)
 224    E = flat.T @ flat
 225
 226    return H, E, P - 1
 227
 228
 229def _rao_f_from_wilks(Lambda: float, K: int, q: int, N: int) -> Tuple[float, int, float]:
 230    """Rao's F-approximation for Wilks' Λ in a one-way design.
 231
 232    Returns (F, df1, df2). With P=3 phases (q=2) and any K≥2 this gives an
 233    exact F distribution; for q=1 or K=1 it reduces to a Hotelling/F test
 234    exactly.
 235    """
 236    # p = K, q = df_effect, N_obs = N*P observations but for repeated-measures
 237    # the relevant sample is N subjects; we follow the standard transformation
 238    # used in textbooks (Tabachnick & Fidell ch. 7) with v = N - 1 - (K + q + 1)/2.
 239    p = K
 240    v = N - 1 - (p + q + 1) / 2.0
 241
 242    pq2 = p * p + q * q
 243    if pq2 == 5:
 244        t = 1.0
 245    else:
 246        denom = pq2 - 5.0
 247        if denom <= 0:
 248            t = 1.0
 249        else:
 250            t = np.sqrt((p * p * q * q - 4.0) / denom)
 251
 252    df1 = int(p * q)
 253    df2 = int(round(v * t - (p * q - 2) / 2.0))
 254    if df2 <= 0:
 255        df2 = max(df1, 1)
 256
 257    if Lambda <= 0.0:
 258        F = float("inf")
 259    else:
 260        Lambda_t = Lambda ** (1.0 / t)
 261        F = ((1.0 - Lambda_t) / Lambda_t) * (df2 / df1)
 262
 263    return F, df1, df2
 264
 265
 266def compute_phi_oneway(
 267    z: np.ndarray,
 268    *,
 269    metric_names: Optional[Sequence[str]] = None,
 270) -> PhiOneway:
 271    """One-way repeated-measures MANOVA on the phase factor.
 272
 273    Parameters
 274    ----------
 275    z : ndarray, shape ``(N, P, K)``
 276        Standardised metric values. Use :func:`standardise` first; do
 277        not pass raw scores.
 278    metric_names : sequence of str, optional
 279        Names of the K metrics in column order. When provided, the
 280        collinearity guard's error messages / warnings identify the
 281        offending metric pair by name; otherwise it falls back to
 282        column indices.
 283
 284    Returns
 285    -------
 286    PhiOneway
 287        Wilks-based Φ, Pillai-based Φ, F-approximation p-value, and the
 288        raw test statistics.
 289
 290    Notes
 291    -----
 292    Paper §3.2 Eq. 2: :math:`\\Phi = 1 - \\eta^2_{phase}`, with
 293    :math:`\\eta^2_{phase} = 1 - \\Lambda^{1/s}` (Wilks) or
 294    :math:`\\eta^2_{phase} = V / s` (Pillai), where
 295    :math:`s = \\min(K, df_{eff})` and :math:`df_{eff} = P - 1`.
 296
 297    The collinearity guard proactively rejects K-vectors whose within-
 298    phase SSCP is numerically ill-conditioned (``cond(E) > 1000``) — a
 299    signature of near-duplicate metrics that would otherwise silently
 300    destabilise Λ. A softer warning fires at ``cond(E) > 30``.
 301    """
 302    if z.ndim != 3:
 303        raise MetricError(
 304            f"compute_phi_oneway expects shape (N, P, K), got {z.shape}",
 305        )
 306    N, P, K = z.shape
 307    if N < 2:
 308        raise MetricError(
 309            f"compute_phi_oneway needs at least 2 scenes, got N={N}",
 310            hint="Drop scenes with missing predictions earlier (complete-case).",
 311        )
 312    if P < 2:
 313        raise MetricError(
 314            f"compute_phi_oneway needs at least 2 phases, got P={P}",
 315        )
 316    if metric_names is not None and len(metric_names) != K:
 317        raise MetricError(
 318            f"metric_names has {len(metric_names)} entries, expected K={K}",
 319        )
 320
 321    H, E, df_effect = _h_e_matrices_oneway(z)
 322    s = min(K, df_effect)
 323
 324    # Collinearity guard — cond(E) diagnoses near-duplicate K-vector metrics
 325    # before the log-determinant path silently produces a garbage Λ.
 326    if K >= 2:
 327        cond_E = float(np.linalg.cond(E))
 328        if not np.isfinite(cond_E) or cond_E > COLLINEARITY_FAIL_COND:
 329            hint = _worst_correlated_pair(z, metric_names)
 330            raise MetricError(
 331                f"Within-phase SSCP E is near-singular "
 332                f"(cond(E) = {cond_E:.2e} > {COLLINEARITY_FAIL_COND:.0e})",
 333                hint=(
 334                    "Two or more K-vector metrics are near-linear combinations "
 335                    "of each other; drop one before running MANOVA. "
 336                    + hint
 337                ).strip(),
 338            )
 339        if cond_E > COLLINEARITY_WARN_COND:
 340            pair = _worst_correlated_pair(z, metric_names)
 341            log.warning(
 342                "compute_phi_oneway: cond(E) = %.1f (threshold %.0f) — "
 343                "moderate collinearity may destabilise Λ. %s",
 344                cond_E, COLLINEARITY_WARN_COND, pair,
 345            )
 346
 347    # Wilks' Λ = |E| / |H + E|.
 348    sign_e, logdet_e = np.linalg.slogdet(E)
 349    sign_he, logdet_he = np.linalg.slogdet(H + E)
 350    if sign_e <= 0 or sign_he <= 0:
 351        raise MetricError(
 352            "Wilks' Λ requires positive-definite E and (H+E); got singular matrices",
 353            hint=(
 354                "This usually means K > N (more metrics than scenes) or perfectly "
 355                "collinear metrics. Drop redundant metrics or increase scene count."
 356            ),
 357        )
 358    Lambda = float(np.exp(logdet_e - logdet_he))
 359    eta2_wilks = 1.0 - Lambda ** (1.0 / s)
 360    phi_wilks = 1.0 - eta2_wilks
 361
 362    # Pillai's trace V = tr( H (H + E)^{-1} ).
 363    HE_inv = np.linalg.solve(H + E, H)
 364    V = float(np.trace(HE_inv))
 365    eta2_pillai = V / s
 366    phi_pillai = 1.0 - eta2_pillai
 367
 368    # F-approximation p-value via Wilks → Rao's F.
 369    try:
 370        from scipy.stats import f as f_dist  # lazy: scipy is an [analysis] extra
 371    except ImportError as e:
 372        raise MetricError(
 373            "scipy is required for Φ p-values",
 374            hint="pip install 'rpx-benchmark[analysis]'",
 375        ) from e
 376    F, df1, df2 = _rao_f_from_wilks(Lambda, K, df_effect, N)
 377    p_value = float(1.0 - f_dist.cdf(F, df1, df2)) if np.isfinite(F) else 0.0
 378
 379    return PhiOneway(
 380        phi_wilks=float(np.clip(phi_wilks, 0.0, 1.0)),
 381        phi_pillai=float(np.clip(phi_pillai, 0.0, 1.0)),
 382        p_value=p_value,
 383        n_eff=N,
 384        K=K,
 385        s=s,
 386        wilks_lambda=Lambda,
 387        pillai_V=V,
 388    )
 389
 390
 391# --------------------------------------------------------------------------- #
 392# Bootstrap confidence interval for Φ
 393# --------------------------------------------------------------------------- #
 394#
 395# The one-shot Rao's-F p-value from ``compute_phi_oneway`` says whether the
 396# phase effect is significant; it does not say how tight the Φ estimate is.
 397# A scene-level nonparametric bootstrap answers that: resample scenes (rows
 398# of the cube) with replacement, recompute Φ per bootstrap, then report the
 399# percentile CI. This is the standard reviewer-defensible way to communicate
 400# uncertainty around a small-sample MANOVA effect size.
 401
 402
 403@dataclass
 404class PhiBootstrapCI:
 405    """Bootstrap distribution of Φ over resampled scene sets.
 406
 407    Attributes
 408    ----------
 409    phi_wilks_mean, phi_pillai_mean : float
 410        Mean of the bootstrap Φ estimates (point estimate, robust to the
 411        original single-fit result).
 412    phi_wilks_ci : tuple[float, float]
 413    phi_pillai_ci : tuple[float, float]
 414        Percentile confidence intervals ``(low, high)`` at the requested
 415        confidence level.
 416    phi_wilks_se, phi_pillai_se : float
 417        Bootstrap standard error (standard deviation of the resampled Φ).
 418    n_boot : int
 419        Number of bootstrap resamples requested.
 420    n_success : int
 421        Number of resamples that yielded a valid Φ (non-singular E, no
 422        cond(E) failure). Resamples that fail are skipped and counted
 423        towards ``n_failed``.
 424    n_failed : int
 425    ci_level : float
 426        Confidence level used (e.g. 0.95).
 427    """
 428
 429    phi_wilks_mean: float
 430    phi_pillai_mean: float
 431    phi_wilks_ci: Tuple[float, float]
 432    phi_pillai_ci: Tuple[float, float]
 433    phi_wilks_se: float
 434    phi_pillai_se: float
 435    n_boot: int
 436    n_success: int
 437    n_failed: int
 438    ci_level: float
 439
 440    def to_dict(self) -> Dict[str, object]:
 441        return {
 442            "phi_wilks_mean": self.phi_wilks_mean,
 443            "phi_pillai_mean": self.phi_pillai_mean,
 444            "phi_wilks_ci_low": self.phi_wilks_ci[0],
 445            "phi_wilks_ci_high": self.phi_wilks_ci[1],
 446            "phi_pillai_ci_low": self.phi_pillai_ci[0],
 447            "phi_pillai_ci_high": self.phi_pillai_ci[1],
 448            "phi_wilks_se": self.phi_wilks_se,
 449            "phi_pillai_se": self.phi_pillai_se,
 450            "n_boot": self.n_boot,
 451            "n_success": self.n_success,
 452            "n_failed": self.n_failed,
 453            "ci_level": self.ci_level,
 454        }
 455
 456
 457def bootstrap_phi_oneway(
 458    z: np.ndarray,
 459    *,
 460    n_boot: int = 1000,
 461    ci_level: float = 0.95,
 462    seed: Optional[int] = None,
 463    metric_names: Optional[Sequence[str]] = None,
 464    min_success_frac: float = 0.5,
 465) -> PhiBootstrapCI:
 466    """Scene-level nonparametric bootstrap of one-way Φ.
 467
 468    Parameters
 469    ----------
 470    z : ndarray, shape ``(N, P, K)``
 471        Standardised metric values (use :func:`standardise` first).
 472        Bootstrap resamples along the scene axis with replacement; the
 473        (phase, metric) coupling within each scene is preserved.
 474    n_boot : int
 475        Number of bootstrap resamples. Defaults to 1000, which is
 476        adequate for percentile CIs at the 95% level.
 477    ci_level : float
 478        Confidence level in ``(0, 1)``. Defaults to 0.95.
 479    seed : int, optional
 480        Reproducible seed for the resampler. ``None`` uses fresh entropy.
 481    metric_names : sequence of str, optional
 482        Forwarded to :func:`compute_phi_oneway` on every resample; used
 483        only for collinearity-guard hints if a resample fails.
 484    min_success_frac : float
 485        Raise :class:`MetricError` if fewer than this fraction of
 486        resamples yield a valid Φ. Default 0.5 — protects against pipes
 487        where standardisation upstream produced degenerate columns.
 488
 489    Returns
 490    -------
 491    PhiBootstrapCI
 492    """
 493    if z.ndim != 3:
 494        raise MetricError(f"bootstrap_phi_oneway expects (N, P, K), got {z.shape}")
 495    N, P, K = z.shape
 496    if N < 2:
 497        raise MetricError(
 498            f"bootstrap_phi_oneway needs at least 2 scenes, got N={N}",
 499        )
 500    if not (0.0 < ci_level < 1.0):
 501        raise MetricError(f"ci_level must be in (0, 1), got {ci_level}")
 502    if n_boot < 50:
 503        raise MetricError(
 504            f"n_boot={n_boot} is too small for a stable percentile CI; use ≥ 50",
 505        )
 506
 507    rng = np.random.default_rng(seed)
 508    wilks_samples: List[float] = []
 509    pillai_samples: List[float] = []
 510    n_failed = 0
 511
 512    for _ in range(n_boot):
 513        idx = rng.integers(0, N, size=N)
 514        z_boot = z[idx]
 515        try:
 516            r = compute_phi_oneway(z_boot, metric_names=metric_names)
 517        except MetricError:
 518            n_failed += 1
 519            continue
 520        wilks_samples.append(r.phi_wilks)
 521        pillai_samples.append(r.phi_pillai)
 522
 523    n_success = len(wilks_samples)
 524    if n_success < int(min_success_frac * n_boot):
 525        raise MetricError(
 526            f"bootstrap failed on {n_failed}/{n_boot} resamples "
 527            f"(only {n_success} succeeded)",
 528            hint=(
 529                "The K-vector is likely marginally-conditioned — many "
 530                "scene resamples yield singular E. Reduce K or increase "
 531                "the scene pool before bootstrapping."
 532            ),
 533        )
 534
 535    wilks_arr = np.array(wilks_samples, dtype=np.float64)
 536    pillai_arr = np.array(pillai_samples, dtype=np.float64)
 537    alpha = 1.0 - ci_level
 538    lo_q, hi_q = 100.0 * (alpha / 2.0), 100.0 * (1.0 - alpha / 2.0)
 539
 540    return PhiBootstrapCI(
 541        phi_wilks_mean=float(wilks_arr.mean()),
 542        phi_pillai_mean=float(pillai_arr.mean()),
 543        phi_wilks_ci=(
 544            float(np.percentile(wilks_arr, lo_q)),
 545            float(np.percentile(wilks_arr, hi_q)),
 546        ),
 547        phi_pillai_ci=(
 548            float(np.percentile(pillai_arr, lo_q)),
 549            float(np.percentile(pillai_arr, hi_q)),
 550        ),
 551        phi_wilks_se=float(wilks_arr.std(ddof=1)) if n_success > 1 else 0.0,
 552        phi_pillai_se=float(pillai_arr.std(ddof=1)) if n_success > 1 else 0.0,
 553        n_boot=n_boot,
 554        n_success=n_success,
 555        n_failed=n_failed,
 556        ci_level=ci_level,
 557    )
 558
 559
 560# --------------------------------------------------------------------------- #
 561# Pairwise Hotelling T² per transition
 562# --------------------------------------------------------------------------- #
 563
 564
 565@dataclass
 566class PhiTransition:
 567    """Result of one pairwise Hotelling T² test.
 568
 569    Attributes
 570    ----------
 571    label : str
 572        Human-readable transition label (e.g. ``"0->1"``).
 573    phi : float
 574        :math:`\\Phi_{i \\to j} = 1 - T^2 / (T^2 + N - 1)`.
 575    eta2 : float
 576        :math:`\\eta^2 = T^2 / (T^2 + N - 1)`.
 577    T2 : float
 578        The Hotelling :math:`T^2` statistic on within-scene paired
 579        differences.
 580    p_value : float
 581        F-converted p-value.
 582    p_holm : Optional[float]
 583        Holm-Bonferroni-adjusted p (set by
 584        :func:`compute_phi_per_transition` when multiple transitions
 585        are tested jointly).
 586    n_eff : int
 587        Number of scenes contributing to the paired test.
 588    """
 589
 590    label: str
 591    phi: float
 592    eta2: float
 593    T2: float
 594    p_value: float
 595    p_holm: Optional[float]
 596    n_eff: int
 597
 598    def to_dict(self) -> Dict[str, object]:
 599        return {
 600            "label": self.label,
 601            "phi": self.phi,
 602            "eta2": self.eta2,
 603            "T2": self.T2,
 604            "p_value": self.p_value,
 605            "p_holm": self.p_holm,
 606            "n_eff": self.n_eff,
 607        }
 608
 609
 610def _hotelling_one_sample(d: np.ndarray) -> Tuple[float, float, float, int, int]:
 611    """One-sample Hotelling T² for paired differences.
 612
 613    Parameters
 614    ----------
 615    d : ndarray, shape ``(N, K)``
 616        Paired-difference vectors (one row per scene).
 617
 618    Returns
 619    -------
 620    T2 : float
 621        ``T^2 = N · d̄' S_d^{-1} d̄``.
 622    p_value : float
 623        F-converted, ``F = ((N-K) / ((N-1) K)) · T^2 ~ F(K, N-K)``.
 624    eta2 : float
 625        ``T^2 / (T^2 + N - 1)``.
 626    df1, df2 : int
 627    """
 628    if d.ndim != 2:
 629        raise MetricError(f"paired differences must be 2D, got {d.shape}")
 630    N, K = d.shape
 631    if N <= K:
 632        raise MetricError(
 633            f"Hotelling T² requires N > K, got N={N}, K={K}",
 634            hint="Drop metrics or increase scene count.",
 635        )
 636    d_bar = d.mean(axis=0)
 637    S_d = np.cov(d, rowvar=False, ddof=1)
 638    try:
 639        S_inv = np.linalg.inv(S_d)
 640    except np.linalg.LinAlgError as e:
 641        raise MetricError(
 642            "Within-difference covariance matrix is singular",
 643            hint="Drop collinear metrics before computing Φ_per_transition.",
 644        ) from e
 645    T2 = float(N * d_bar @ S_inv @ d_bar)
 646
 647    # F-distribution conversion.
 648    df1 = K
 649    df2 = N - K
 650    F = ((N - K) / ((N - 1) * K)) * T2
 651
 652    try:
 653        from scipy.stats import f as f_dist
 654    except ImportError as e:
 655        raise MetricError(
 656            "scipy is required for Hotelling T² p-values",
 657            hint="pip install 'rpx-benchmark[analysis]'",
 658        ) from e
 659    p_value = float(1.0 - f_dist.cdf(F, df1, df2))
 660
 661    eta2 = T2 / (T2 + N - 1)
 662    return T2, p_value, eta2, df1, df2
 663
 664
 665def _holm_bonferroni(p_values: list[float]) -> list[float]:
 666    """Holm-Bonferroni step-down adjustment.
 667
 668    Parameters
 669    ----------
 670    p_values : list of float
 671        Raw p-values (any order).
 672
 673    Returns
 674    -------
 675    list of float
 676        Adjusted p-values in the same order as input. Monotonically
 677        non-decreasing under the Holm rule and clipped to ``[0, 1]``.
 678    """
 679    m = len(p_values)
 680    order = np.argsort(p_values)
 681    adjusted = np.empty(m, dtype=np.float64)
 682    running_max = 0.0
 683    for rank_i, idx in enumerate(order):
 684        adj = (m - rank_i) * p_values[idx]
 685        running_max = max(running_max, adj)
 686        adjusted[idx] = min(running_max, 1.0)
 687    return adjusted.tolist()
 688
 689
 690def compute_phi_per_transition(
 691    z: np.ndarray,
 692    *,
 693    transitions: Tuple[Tuple[int, int], ...] = ((0, 1), (1, 2), (0, 2)),
 694    apply_holm: bool = True,
 695) -> Dict[str, PhiTransition]:
 696    """Pairwise Hotelling T² → Φ for each phase transition.
 697
 698    Parameters
 699    ----------
 700    z : ndarray, shape ``(N, P, K)``
 701        Standardised metric values.
 702    transitions : tuple of (p1, p2)
 703        Phase-index pairs to test. Defaults to all three: Clu→Int,
 704        Int→Cln, Clu→Cln.
 705    apply_holm : bool
 706        Apply Holm-Bonferroni correction across the requested
 707        transitions (paper §3.2: family-wise control for the 3 pairwise
 708        transitions per (model, task)).
 709
 710    Returns
 711    -------
 712    dict
 713        ``{"0->1": PhiTransition, "1->2": ..., "0->2": ...}``.
 714    """
 715    if z.ndim != 3:
 716        raise MetricError(f"expects (N, P, K), got {z.shape}")
 717    N, P, _ = z.shape
 718
 719    raw_results: Dict[str, PhiTransition] = {}
 720    raw_ps: list[float] = []
 721    keys: list[str] = []
 722    for p1, p2 in transitions:
 723        if not (0 <= p1 < P) or not (0 <= p2 < P):
 724            raise MetricError(
 725                f"transition ({p1}, {p2}) out of range for P={P}",
 726            )
 727        if p1 == p2:
 728            raise MetricError(f"degenerate transition ({p1}, {p1})")
 729        d = z[:, p2, :] - z[:, p1, :]
 730        T2, p_raw, eta2, _, _ = _hotelling_one_sample(d)
 731        phi = 1.0 - eta2
 732        key = f"{p1}->{p2}"
 733        raw_results[key] = PhiTransition(
 734            label=key,
 735            phi=float(np.clip(phi, 0.0, 1.0)),
 736            eta2=float(np.clip(eta2, 0.0, 1.0)),
 737            T2=T2,
 738            p_value=p_raw,
 739            p_holm=None,
 740            n_eff=N,
 741        )
 742        raw_ps.append(p_raw)
 743        keys.append(key)
 744
 745    if apply_holm and raw_ps:
 746        adj = _holm_bonferroni(raw_ps)
 747        for k, p_h in zip(keys, adj, strict=True):
 748            r = raw_results[k]
 749            raw_results[k] = PhiTransition(
 750                label=r.label,
 751                phi=r.phi,
 752                eta2=r.eta2,
 753                T2=r.T2,
 754                p_value=r.p_value,
 755                p_holm=p_h,
 756                n_eff=r.n_eff,
 757            )
 758
 759    return raw_results
 760
 761
 762# --------------------------------------------------------------------------- #
 763# Interpretation rubric (Cohen anchors)
 764# --------------------------------------------------------------------------- #
 765
 766
 767def phi_interpretation(phi: float) -> str:
 768    """Map a Φ value to a Cohen partial-η² effect-size label.
 769
 770    Thresholds (paper §3.2):
 771
 772    - ``η² < 0.01``  →  Φ > 0.99   →  "negligible"
 773    - ``0.01 ≤ η² < 0.06``  →  0.94 < Φ ≤ 0.99  →  "small"
 774    - ``0.06 ≤ η² < 0.14``  →  0.86 < Φ ≤ 0.94  →  "medium"
 775    - ``η² ≥ 0.14``  →  Φ ≤ 0.86  →  "large"
 776    """
 777    if not (0.0 <= phi <= 1.0):
 778        raise MetricError(f"Φ must be in [0, 1], got {phi}")
 779    if phi > 0.99:
 780        return "negligible"
 781    if phi > 0.94:
 782        return "small"
 783    if phi > 0.86:
 784        return "medium"
 785    return "large"
 786
 787
 788# --------------------------------------------------------------------------- #
 789# Two-way mixed-design MANOVA (paper §3.2)
 790# --------------------------------------------------------------------------- #
 791#
 792# Phase is within-subjects (each scene appears in all P phases), difficulty
 793# tier is between-subjects (each scene belongs to exactly one tier). The
 794# three effects we report:
 795#
 796#   - Phase main effect:           does the per-scene phase profile differ
 797#                                  from a flat one, pooled across tiers?
 798#   - Difficulty main effect:      do tier-averaged scene means differ?
 799#   - Phase × difficulty interaction: does the phase profile differ across
 800#                                  tiers? (The headline number paper §3.2
 801#                                  calls η²_{phase×diff}.)
 802#
 803# Standard recipe: build orthonormal phase contrasts C (P × (P−1)); for each
 804# scene, form the contrast vector c_s = C^T x_s ∈ R^{(P−1) × K}, flatten to
 805# R^{(P−1)K}. Then:
 806#
 807#   - Phase main effect  ≈ one-sample Hotelling on {c_s} vs. zero.
 808#   - Phase × diff       ≈ one-way between-groups MANOVA on {c_s} by tier.
 809#   - Diff main effect   ≈ one-way between-groups MANOVA on
 810#                          {m_s = mean_p x_{s,p}} by tier.
 811#
 812# This is the standard reduction used in repeated-measures multivariate
 813# analysis (Tabachnick & Fidell, ch. 8). It's not the only formulation
 814# possible but it's clean and the effect-size η² values match what the
 815# paper §3.2 asks for.
 816
 817
 818@dataclass
 819class PhiMixedDesign:
 820    """Result of the two-way mixed-design MANOVA.
 821
 822    Attributes
 823    ----------
 824    phi_per_tier : Dict[str, PhiOneway]
 825        Per-tier one-way Φ (computed by restricting the cell tensor to
 826        each tier and running :func:`compute_phi_oneway`).
 827        Paper notation: $\\Phi_\\text{E}, \\Phi_\\text{M}, \\Phi_\\text{H}$.
 828    eta2_phase : float
 829        Main-effect partial $\\eta^2$ for the phase factor, pooled
 830        across tiers.
 831    eta2_diff : float
 832        Main-effect partial $\\eta^2$ for the difficulty factor.
 833    eta2_interaction : float
 834        Phase × difficulty interaction partial $\\eta^2$. The headline
 835        "do harder scenes amplify the phase effect?" number.
 836    phi_phase, phi_diff, phi_interaction : float
 837        Complements of the η² values, ∈ [0, 1].
 838    p_phase, p_diff, p_interaction : float
 839        F-approximation p-values.
 840    K : int
 841        Metric dimension carried into the contrast space.
 842    n_scenes : int
 843        Number of scenes that contributed.
 844    tier_counts : Dict[str, int]
 845        Per-tier scene counts.
 846    notes : List[str]
 847    """
 848
 849    phi_per_tier: Dict[str, PhiOneway]
 850    eta2_phase: float
 851    eta2_diff: float
 852    eta2_interaction: float
 853    phi_phase: float
 854    phi_diff: float
 855    phi_interaction: float
 856    p_phase: float
 857    p_diff: float
 858    p_interaction: float
 859    K: int
 860    n_scenes: int
 861    tier_counts: Dict[str, int]
 862    notes: List[str] = field(default_factory=list)
 863
 864    def to_dict(self) -> Dict[str, object]:
 865        return {
 866            "phi_per_tier": {t: r.to_dict() for t, r in self.phi_per_tier.items()},
 867            "eta2_phase": self.eta2_phase,
 868            "eta2_diff": self.eta2_diff,
 869            "eta2_interaction": self.eta2_interaction,
 870            "phi_phase": self.phi_phase,
 871            "phi_diff": self.phi_diff,
 872            "phi_interaction": self.phi_interaction,
 873            "p_phase": self.p_phase,
 874            "p_diff": self.p_diff,
 875            "p_interaction": self.p_interaction,
 876            "K": self.K,
 877            "n_scenes": self.n_scenes,
 878            "tier_counts": dict(self.tier_counts),
 879            "notes": list(self.notes),
 880        }
 881
 882
 883def _orthonormal_phase_contrasts(P: int) -> np.ndarray:
 884    """Return a $P \\times (P-1)$ matrix whose columns are an
 885    orthonormal basis for the subspace orthogonal to $(1,\\dots,1)$.
 886
 887    Built via QR of the centering matrix; numerically stable for any
 888    reasonable P. Caller does not need to know the specific basis ---
 889    Λ, V, and η² are invariant under orthonormal change of basis on
 890    the contrast subspace.
 891    """
 892    if P < 2:
 893        raise MetricError(f"need P>=2 phases for contrasts, got {P}")
 894    centering = np.eye(P) - np.ones((P, P)) / P
 895    Q, _ = np.linalg.qr(centering[:, : P - 1])
 896    return Q  # shape (P, P-1)
 897
 898
 899def _one_way_between_groups(
 900    X: np.ndarray,
 901    group_labels: Sequence[Any],
 902) -> Tuple[np.ndarray, np.ndarray, int, int]:
 903    """Build H, E matrices for a one-way between-groups MANOVA.
 904
 905    Parameters
 906    ----------
 907    X : (N, D) ndarray
 908        Per-subject observation vectors.
 909    group_labels : length-N sequence
 910        Group assignment per row.
 911
 912    Returns
 913    -------
 914    H : (D, D) between-groups SSCP.
 915    E : (D, D) within-groups (pooled) SSCP.
 916    df_h : int
 917        Hypothesis degrees of freedom: G - 1.
 918    G : int
 919        Number of distinct groups.
 920    """
 921    N, D = X.shape
 922    if N != len(group_labels):
 923        raise MetricError(
 924            f"length mismatch: X has {N} rows, group_labels has {len(group_labels)}"
 925        )
 926    unique = sorted(set(group_labels), key=str)
 927    G = len(unique)
 928    if G < 2:
 929        raise MetricError(
 930            f"between-groups MANOVA needs >=2 groups, got {G}",
 931            hint="The interaction term is undefined with a single tier.",
 932        )
 933
 934    grand = X.mean(axis=0)
 935    H = np.zeros((D, D), dtype=np.float64)
 936    E = np.zeros((D, D), dtype=np.float64)
 937    for g in unique:
 938        mask = np.array([gl == g for gl in group_labels], dtype=bool)
 939        n_g = int(mask.sum())
 940        if n_g < 1:
 941            continue
 942        X_g = X[mask]
 943        mean_g = X_g.mean(axis=0)
 944        diff_g = (mean_g - grand)[:, None]
 945        H += n_g * (diff_g @ diff_g.T)
 946        centered = X_g - mean_g
 947        E += centered.T @ centered
 948
 949    return H, E, G - 1, G
 950
 951
 952def _wilks_eta_squared_from_he(
 953    H: np.ndarray, E: np.ndarray, K: int, df_h: int, n_obs: int
 954) -> Tuple[float, float, float]:
 955    """Convert (H, E, df_h, n_obs) to (eta²_wilks, F, p_value)."""
 956    sign_e, logdet_e = np.linalg.slogdet(E)
 957    sign_he, logdet_he = np.linalg.slogdet(H + E)
 958    if sign_e <= 0 or sign_he <= 0:
 959        raise MetricError(
 960            "Wilks Λ requires positive-definite E and (H+E); got singular matrices",
 961            hint=(
 962                "Drop collinear metrics, increase scene count, or merge "
 963                "tiers — the input is too rank-deficient for the test."
 964            ),
 965        )
 966    Lambda = float(np.exp(logdet_e - logdet_he))
 967    s = max(1, min(K, df_h))
 968    eta2 = 1.0 - Lambda ** (1.0 / s)
 969    F, df1, df2 = _rao_f_from_wilks(Lambda, K, df_h, n_obs)
 970
 971    try:
 972        from scipy.stats import f as f_dist
 973    except ImportError as e:
 974        raise MetricError(
 975            "scipy required for F-approximation p-values",
 976            hint="pip install 'rpx-benchmark[analysis]'",
 977        ) from e
 978
 979    p_value = float(1.0 - f_dist.cdf(F, df1, df2)) if np.isfinite(F) else 0.0
 980    return float(np.clip(eta2, 0.0, 1.0)), F, p_value
 981
 982
 983def compute_phi_mixed_design(
 984    z: np.ndarray,
 985    tier_by_scene: Sequence[str],
 986) -> PhiMixedDesign:
 987    """Two-way mixed-design MANOVA: phase within, tier between.
 988
 989    Parameters
 990    ----------
 991    z : ndarray, shape ``(N, P, K)``
 992        Standardised metric values (use :func:`standardise` first).
 993    tier_by_scene : length-N sequence
 994        Tier label per scene (e.g. ``"easy"``, ``"medium"``, ``"hard"``).
 995        Must align with the scene axis of ``z``.
 996
 997    Returns
 998    -------
 999    PhiMixedDesign
1000    """
1001    if z.ndim != 3:
1002        raise MetricError(f"compute_phi_mixed_design expects (N, P, K), got {z.shape}")
1003    N, P, K = z.shape
1004    if N != len(tier_by_scene):
1005        raise MetricError(
1006            f"tier_by_scene length {len(tier_by_scene)} != N={N}",
1007        )
1008    if P < 2:
1009        raise MetricError(f"need P>=2 phases, got {P}")
1010    if N < 4:
1011        raise MetricError(
1012            f"need at least 4 scenes for a meaningful mixed design, got {N}",
1013        )
1014
1015    notes: List[str] = []
1016
1017    # ---- Per-tier Φ (just restrict compute_phi_oneway to each subset) ----
1018    tier_counts: Dict[str, int] = {}
1019    phi_per_tier: Dict[str, PhiOneway] = {}
1020    tier_indices: Dict[str, List[int]] = {}
1021    for i, t in enumerate(tier_by_scene):
1022        tier_indices.setdefault(str(t), []).append(i)
1023    for tier, indices in tier_indices.items():
1024        tier_counts[tier] = len(indices)
1025        if len(indices) < 2:
1026            notes.append(f"tier {tier!r}: only {len(indices)} scene(s); skipping per-tier Φ")
1027            continue
1028        try:
1029            phi_per_tier[tier] = compute_phi_oneway(z[indices])
1030        except MetricError as e:
1031            notes.append(f"tier {tier!r}: compute_phi_oneway failed ({e})")
1032
1033    # ---- Build contrast vectors c_s ∈ R^{(P-1)*K} per scene ----
1034    C = _orthonormal_phase_contrasts(P)  # (P, P-1)
1035    # c_s = C^T x_s gives (P-1, K); flatten to ((P-1)*K,)
1036    c = np.einsum("ip,npk->nik", C.T, z)  # (N, P-1, K)
1037    D_c = (P - 1) * K
1038    c_flat = c.reshape(N, D_c)
1039
1040    # ---- Phase main effect: one-sample Hotelling on c_s vs zero, pooled
1041    # over all scenes. The pool may include tier-mean shifts; we therefore
1042    # center within-tier before testing the grand mean, which is the
1043    # standard within-subjects-only main-effect formulation.
1044    # H_phase = N * c̄_centered c̄_centered^T,  E_phase = within-tier SSCP.
1045    eta2_phase = phi_phase = p_phase = 0.0
1046    try:
1047        unique_tiers = sorted(set(str(t) for t in tier_by_scene), key=str)
1048        H_phase = np.zeros((D_c, D_c))
1049        E_phase = np.zeros((D_c, D_c))
1050        tier_centered_means_weighted_sum = np.zeros(D_c)
1051        n_total = 0
1052        for tier in unique_tiers:
1053            mask = np.array([str(t) == tier for t in tier_by_scene], dtype=bool)
1054            c_g = c_flat[mask]
1055            if c_g.shape[0] == 0:
1056                continue
1057            mean_g = c_g.mean(axis=0)
1058            tier_centered_means_weighted_sum += c_g.shape[0] * mean_g
1059            n_total += c_g.shape[0]
1060            centered_g = c_g - mean_g
1061            E_phase += centered_g.T @ centered_g
1062        grand_c = tier_centered_means_weighted_sum / max(1, n_total)
1063        # H_phase tests grand_c == 0 (subject-pooled phase profile is flat).
1064        # With df_h = 1, we get a within-subjects T² style test.
1065        H_phase = n_total * np.outer(grand_c, grand_c)
1066        eta2_phase, _F_phase, p_phase = _wilks_eta_squared_from_he(
1067            H_phase, E_phase, K=D_c, df_h=1, n_obs=n_total
1068        )
1069        phi_phase = float(np.clip(1.0 - eta2_phase, 0.0, 1.0))
1070    except MetricError as e:
1071        notes.append(f"phase main effect: {e}")
1072
1073    # ---- Phase × difficulty interaction: one-way between-groups MANOVA
1074    # on the contrast vectors {c_s} with tier as the grouping factor.
1075    eta2_int = phi_int = p_int = 0.0
1076    try:
1077        H_int, E_int, df_h_int, _G = _one_way_between_groups(c_flat, list(map(str, tier_by_scene)))
1078        eta2_int, _F_int, p_int = _wilks_eta_squared_from_he(
1079            H_int, E_int, K=D_c, df_h=df_h_int, n_obs=N
1080        )
1081        phi_int = float(np.clip(1.0 - eta2_int, 0.0, 1.0))
1082    except MetricError as e:
1083        notes.append(f"phase × difficulty interaction: {e}")
1084
1085    # ---- Difficulty main effect: one-way between-groups MANOVA on
1086    # phase-averaged scene means.
1087    eta2_diff = phi_diff = p_diff = 0.0
1088    try:
1089        m_s = z.mean(axis=1)  # (N, K)
1090        H_d, E_d, df_h_d, _G2 = _one_way_between_groups(m_s, list(map(str, tier_by_scene)))
1091        eta2_diff, _F_d, p_diff = _wilks_eta_squared_from_he(
1092            H_d, E_d, K=K, df_h=df_h_d, n_obs=N
1093        )
1094        phi_diff = float(np.clip(1.0 - eta2_diff, 0.0, 1.0))
1095    except MetricError as e:
1096        notes.append(f"difficulty main effect: {e}")
1097
1098    return PhiMixedDesign(
1099        phi_per_tier=phi_per_tier,
1100        eta2_phase=eta2_phase,
1101        eta2_diff=eta2_diff,
1102        eta2_interaction=eta2_int,
1103        phi_phase=phi_phase,
1104        phi_diff=phi_diff,
1105        phi_interaction=phi_int,
1106        p_phase=p_phase,
1107        p_diff=p_diff,
1108        p_interaction=p_int,
1109        K=K,
1110        n_scenes=N,
1111        tier_counts=tier_counts,
1112        notes=notes,
1113    )
1114
1115
1116__all__ = [
1117    "COLLINEARITY_FAIL_COND",
1118    "COLLINEARITY_WARN_COND",
1119    "PhiBootstrapCI",
1120    "PhiMixedDesign",
1121    "PhiOneway",
1122    "PhiTransition",
1123    "bootstrap_phi_oneway",
1124    "compute_phi_mixed_design",
1125    "compute_phi_oneway",
1126    "compute_phi_per_transition",
1127    "phi_interpretation",
1128    "standardise",
1129]