Add permutation-column VIP / CV-beta null for the observational relate

CHANGELOG.md

@@ -11,6 +11,24 @@ those changes.
 
 ## [Unreleased]
 
+## [1.58.0] - 2026-07-23
+
+### Added
+
+- `process_improve.sensory.permutation_column_null`: an empirical VIP /
+  cross-validated-beta null for the observational relate's descriptor block. It adds
+  permuted "knockoff" columns - each a row-shuffled copy of a real descriptor, so the
+  null matches the data's own marginals rather than a simulated distribution - refits
+  the PLS relate over many permutations, and reports, per surviving descriptor, whether
+  its VIP and cross-validated beta clear a high quantile of the null band. This
+  calibrates driver magnitude against the data's own permuted columns in the `p >> n`
+  regime, where even an unrelated descriptor earns a non-trivial VIP. The function is
+  decoupled from the influence gate: it takes an `ignore` list of descriptor names (by
+  name, validated - an unknown name raises) that are dropped from the fit entirely, so
+  a caller can remove the gate-demoted spikes before calibrating the survivors. The
+  knockoff count scales with the block (`fraction`, with a `min_knockoffs` floor for
+  narrow data and an optional cap), and `n_iter` sets the permutation count.
+
 ## [1.57.0] - 2026-07-23
 
 ### Added
@@ -2581,7 +2599,8 @@ this entry records them together.
 - Reworked the README with a sharper value proposition and a
   "Why not scikit-learn?" comparison table.
 
-[Unreleased]: https://github.com/kgdunn/process-improve/compare/v1.57.0...HEAD
+[Unreleased]: https://github.com/kgdunn/process-improve/compare/v1.58.0...HEAD
+[1.58.0]: https://github.com/kgdunn/process-improve/compare/v1.57.0...v1.58.0
 [1.57.0]: https://github.com/kgdunn/process-improve/compare/v1.56.0...v1.57.0
 [1.56.0]: https://github.com/kgdunn/process-improve/compare/v1.55.1...v1.56.0
 [1.55.1]: https://github.com/kgdunn/process-improve/compare/v1.55.0...v1.55.1

CITATION.cff

@@ -12,7 +12,7 @@ authors:
 repository-code: "https://github.com/kgdunn/process-improve"
 url: "https://kgdunn.github.io/process-improve/"
 license: MIT
-version: 1.57.0
+version: 1.58.0
 date-released: "2026-07-23"
 keywords:
   - chemometrics

pyproject.toml

@@ -1,6 +1,6 @@
 [project]
 name = "process-improve"
-version = "1.57.0"
+version = "1.58.0"
 description = 'Designed Experiments; Latent Variables (PCA, PLS, multivariate methods with missing data); Process Monitoring; Batch data analysis.'
 readme = "README.md"
 license = "MIT"

src/process_improve/sensory/__init__.py

@@ -20,6 +20,7 @@
     aggregate_to_product,
     analyze_descriptive,
     discriminate_observational,
+    permutation_column_null,
     product_means,
     relate_designed,
     relate_observational,
@@ -48,6 +49,7 @@
     "discriminate_observational",
     "mixed_assessor_model",
     "panel_scorecard",
+    "permutation_column_null",
     "product_means",
     "relate_designed",
     "relate_observational",

src/process_improve/sensory/analysis.py

@@ -38,6 +38,7 @@
 from scipy.stats import pearsonr
 from scipy.stats import t as t_dist
 
+from process_improve.multivariate._common import NotEnoughVarianceError
 from process_improve.multivariate.methods import PCA, PLS, MCUVScaler, selectivity_ratio, vip
 from process_improve.sensory.mam import MAMResult, align_scores, mixed_assessor_model
 from process_improve.sensory.panel import PanelScorecard, apply_correction, panel_scorecard
@@ -573,6 +574,165 @@ def relate_observational(  # noqa: PLR0913
     return result
 
 
+def _knockoff_block(x_block: pd.DataFrame, k: int, rng: np.random.Generator) -> pd.DataFrame:
+    """Return ``k`` null columns, each a row-permutation of a randomly chosen real column.
+
+    Sampling the source columns with replacement lets ``k`` exceed the number of real
+    columns (needed for the ``min_knockoffs`` floor on a narrow block); a permuted copy
+    keeps the source column's own marginal - its spread, sparsity, and support - so the
+    null is matched to the data rather than to a simulated distribution.
+    """
+    n_rows, p = x_block.shape
+    source = rng.integers(0, p, size=k)
+    columns = {
+        f"__null_{j}": x_block.iloc[:, src].to_numpy()[rng.permutation(n_rows)]
+        for j, src in enumerate(source)
+    }
+    return pd.DataFrame(columns, index=x_block.index)
+
+
+def permutation_column_null(  # noqa: PLR0913
+    agg: pd.DataFrame,
+    covariates: pd.DataFrame,
+    *,
+    ignore: list[str] | None = None,
+    n_components: int = 2,
+    fraction: float = 0.15,
+    min_knockoffs: int = 7,
+    max_knockoffs: int | None = None,
+    n_iter: int = 200,
+    quantile: float = 0.95,
+    random_state: int = 0,
+) -> dict[str, Any]:
+    """Empirical VIP / cross-validated-beta null for the descriptor block.
+
+    Adds ``k`` permuted "knockoff" columns - each a row-shuffled copy of a real
+    descriptor (:func:`_knockoff_block`) - to the descriptor block, fits the PLS relate,
+    and reads the VIP and cross-validated beta of every column. Repeated over ``n_iter``
+    permutations, the knockoff columns form an empirical null band: a real descriptor is
+    only credible if its VIP / beta clears a high quantile of the null the permuted
+    columns achieve. Because even a descriptor with no real relationship earns a
+    non-trivial VIP in a ``p >> n`` fit, this calibrates the magnitude against the data's
+    own permuted columns rather than a parametric cutoff.
+
+    This is decoupled from the influence gate: it does not itself remove any descriptors.
+    Pass the descriptors the gate demoted (single or twin-support spikes) as ``ignore``;
+    they are dropped from the fit entirely - not merely skipped when building knockoffs -
+    so they no longer distort the scores or VIP of the survivors.
+
+    Parameters
+    ----------
+    agg : pandas.DataFrame
+        Product-by-attribute mean table (index ``product``).
+    covariates : pandas.DataFrame
+        One row per product with the measured descriptors (index ``product``; a
+        ``product`` column, if present, is ignored).
+    ignore : list of str, optional
+        Descriptor names to drop from the fit before building the null (default: none).
+        A name absent from the descriptor block raises ``ValueError`` so a typo fails
+        loudly instead of silently doing nothing.
+    n_components : int
+        Latent components for each PLS fit.
+    fraction : float
+        Fraction of surviving descriptors used as the knockoff count.
+    min_knockoffs : int
+        Floor on the knockoff count, so a narrow block still gets a usable null.
+    max_knockoffs : int, optional
+        Optional cap on the knockoff count (default: uncapped).
+    n_iter : int
+        Number of permutations (refits); more gives a smoother null threshold.
+    quantile : float
+        Null quantile used as the significance threshold (e.g. 0.95).
+    random_state : int
+        Seed for the permutations.
+
+    Returns
+    -------
+    dict
+        ``descriptors`` (per surviving descriptor: ``vip`` / ``cv_beta`` and the
+        ``*_null_threshold`` and ``*_exceeds_null`` fields), plus the settings used and
+        the counts (``n_descriptors``, ``n_knockoffs``, ``n_iter``, ``ignored``).
+    """
+    if fraction <= 0.0:
+        raise ValueError(f"fraction must be positive, got {fraction}.")
+    if not 0.0 < quantile < 1.0:
+        raise ValueError(f"quantile must be in (0, 1), got {quantile}.")
+    if min_knockoffs < 1 or n_iter < 1:
+        raise ValueError("min_knockoffs and n_iter must both be at least 1.")
+
+    ignore = list(ignore or [])
+    x_all = covariates.loc[agg.index]
+    descriptors = [c for c in x_all.columns if c != "product" and pd.api.types.is_numeric_dtype(x_all[c])]
+    unknown = sorted(set(ignore) - set(descriptors))
+    if unknown:
+        raise ValueError(f"ignore contains descriptor name(s) not in the covariate block: {unknown}.")
+
+    kept = [c for c in descriptors if c not in set(ignore)]
+    p = len(kept)
+    config = {
+        "n_descriptors": p,
+        "n_ignored": len(ignore),
+        "ignored": sorted(ignore),
+        "fraction": fraction,
+        "quantile": quantile,
+        "n_iter_requested": n_iter,
+    }
+    if p < 2:
+        return {"ok": False, "reason": f"need at least 2 descriptors after ignoring; got {p}.", **config}
+
+    y_block = agg.astype(float)
+    x_block = x_all[kept].astype(float)
+    n_rows = x_block.shape[0]
+    k = max(int(min_knockoffs), round(fraction * p))
+    if max_knockoffs is not None:
+        k = min(k, int(max_knockoffs))
+    max_comp = max(1, min(n_components, p, n_rows - 1))
+    null_names = [f"__null_{j}" for j in range(k)]
+
+    rng = np.random.default_rng(random_state)
+    real_vip_runs: list[pd.Series] = []
+    real_beta_runs: list[pd.Series] = []
+    null_vip: list[float] = []
+    null_beta: list[float] = []
+    for _ in range(n_iter):
+        x_aug = pd.concat([x_block, _knockoff_block(x_block, k, rng)], axis=1)
+        try:
+            with warnings.catch_warnings():
+                warnings.simplefilter("ignore", RuntimeWarning)
+                pls = PLS(n_components=max_comp).fit(x_aug, y_block)
+                beta = pls.cross_validate(x_aug, y_block, cv="loo", show_progress=False).beta_mean.abs().max(axis=1)
+        except (np.linalg.LinAlgError, NotEnoughVarianceError):
+            continue  # a degenerate (singular / no-variance) block contributes nothing to the null
+        vips = vip(pls)
+        real_vip_runs.append(vips.reindex(kept))
+        real_beta_runs.append(beta.reindex(kept))
+        null_vip.extend(float(vips[name]) for name in null_names)
+        null_beta.extend(float(beta[name]) for name in null_names)
+
+    if not real_vip_runs:
+        return {"ok": False, "reason": "every permuted fit was singular; no null could be formed.", **config}
+
+    real_vip = pd.concat(real_vip_runs, axis=1).mean(axis=1)
+    real_beta = pd.concat(real_beta_runs, axis=1).mean(axis=1)
+    vip_threshold = float(np.nanquantile(null_vip, quantile))
+    beta_threshold = float(np.nanquantile(null_beta, quantile))
+
+    records: list[dict[str, Any]] = [
+        {
+            "descriptor": desc,
+            "vip": float(real_vip[desc]),
+            "vip_null_threshold": vip_threshold,
+            "vip_exceeds_null": bool(real_vip[desc] > vip_threshold),
+            "cv_beta": float(real_beta[desc]),
+            "cv_beta_null_threshold": beta_threshold,
+            "cv_beta_exceeds_null": bool(real_beta[desc] > beta_threshold),
+        }
+        for desc in kept
+    ]
+    records.sort(key=lambda r: r["vip"], reverse=True)
+    return {"ok": True, "descriptors": records, "n_knockoffs": k, "n_iter": len(real_vip_runs), **config}
+
+
 def product_means(panel: pd.DataFrame, conf_level: float = 0.95) -> pd.DataFrame:
     """Return per product-by-attribute mean with a confidence interval."""
     rows: list[dict[str, Any]] = []

tests/test_sensory.py

@@ -22,6 +22,7 @@
     _collinear_clusters,
     _jackknife_correlation,
     discriminate_observational,
+    permutation_column_null,
     relate_designed,
 )
 from process_improve.sensory.ingest import reshape_to_long
@@ -552,6 +553,112 @@ def test_discriminator_small_sample_skips_jackknife_gate():
     assert not desc["discriminator_significant"].any()
 
 
+def _null_case(seed: int = 0) -> tuple[pd.DataFrame, pd.DataFrame]:
+    """Product-mean table + covariates with two genuine drivers and six noise columns."""
+    rng = np.random.default_rng(seed)
+    n = 20
+    products = [f"P{i}" for i in range(n)]
+    driver1 = np.linspace(0.0, 1.0, n) + rng.normal(0, 0.03, n)
+    driver2 = rng.normal(0, 1, n)
+    cov = pd.DataFrame(
+        {"driver1": driver1, "driver2": driver2, **{f"noise{j}": rng.normal(0, 1, n) for j in range(6)}},
+        index=products,
+    )
+    cov.index.name = "product"
+    agg = pd.DataFrame(
+        {"A": 3.0 * driver1 + rng.normal(0, 0.2, n), "B": 2.0 * driver2 + rng.normal(0, 0.2, n)},
+        index=products,
+    )
+    return agg, cov
+
+
+def test_permutation_null_flags_drivers_over_noise():
+    """Genuine drivers clear the permuted-column null; noise columns fall below them."""
+    agg, cov = _null_case()
+    result = permutation_column_null(agg, cov, n_iter=12, min_knockoffs=7, random_state=0)
+    assert result["ok"]
+    assert result["n_descriptors"] == 8
+    assert result["n_knockoffs"] == 7
+    by_name = {r["descriptor"]: r for r in result["descriptors"]}
+
+    for driver in ("driver1", "driver2"):
+        assert by_name[driver]["vip_exceeds_null"]
+        assert by_name[driver]["cv_beta_exceeds_null"]
+
+    # Every noise column scores below both drivers on VIP (robust to platform jitter near
+    # the threshold), and no noise column clears the beta null.
+    min_driver_vip = min(by_name["driver1"]["vip"], by_name["driver2"]["vip"])
+    for j in range(6):
+        noise = by_name[f"noise{j}"]
+        assert noise["vip"] < min_driver_vip
+        assert not noise["cv_beta_exceeds_null"]
+
+
+def test_permutation_null_ignore_drops_columns():
+    """`ignore` removes named descriptors from the fit and the output."""
+    agg, cov = _null_case()
+    result = permutation_column_null(agg, cov, ignore=["noise0", "noise1"], n_iter=5, random_state=0)
+    names = [r["descriptor"] for r in result["descriptors"]]
+    assert result["n_descriptors"] == 6
+    assert result["ignored"] == ["noise0", "noise1"]
+    assert "noise0" not in names
+    assert "noise1" not in names
+
+
+def test_permutation_null_unknown_ignore_name_raises():
+    """A descriptor name that is not in the block fails loudly rather than silently."""
+    agg, cov = _null_case()
+    with pytest.raises(ValueError, match="not in the covariate block"):
+        permutation_column_null(agg, cov, ignore=["not_a_real_column"], n_iter=2)
+
+
+def test_permutation_null_is_deterministic():
+    """The same seed gives identical thresholds and flags."""
+    agg, cov = _null_case()
+    first = permutation_column_null(agg, cov, n_iter=5, random_state=1)
+    second = permutation_column_null(agg, cov, n_iter=5, random_state=1)
+    assert first["descriptors"] == second["descriptors"]
+
+
+def test_permutation_null_validates_parameters():
+    """Out-of-range fraction / quantile / counts are rejected."""
+    agg, cov = _null_case()
+    with pytest.raises(ValueError, match="fraction must be positive"):
+        permutation_column_null(agg, cov, fraction=0.0, n_iter=2)
+    with pytest.raises(ValueError, match="quantile must be in"):
+        permutation_column_null(agg, cov, quantile=1.5, n_iter=2)
+    with pytest.raises(ValueError, match="at least 1"):
+        permutation_column_null(agg, cov, min_knockoffs=0, n_iter=2)
+
+
+def test_permutation_null_max_knockoffs_caps_the_count():
+    """`max_knockoffs` bounds the knockoff count even when the fraction is larger."""
+    agg, cov = _null_case()
+    result = permutation_column_null(agg, cov, fraction=0.9, max_knockoffs=3, n_iter=3, random_state=0)
+    assert result["n_knockoffs"] == 3
+
+
+def test_permutation_null_too_few_descriptors_returns_not_ok():
+    """Fewer than two descriptors after ignoring cannot form a null."""
+    agg, cov = _null_case()
+    keep_two = cov[["driver1", "driver2"]]
+    result = permutation_column_null(agg, keep_two, ignore=["driver2"], n_iter=2)
+    assert not result["ok"]
+    assert "at least 2 descriptors" in result["reason"]
+
+
+def test_permutation_null_degenerate_block_returns_not_ok():
+    """A no-variance descriptor block degrades gracefully instead of crashing."""
+    products = [f"P{i}" for i in range(8)]
+    cov = pd.DataFrame(
+        {"c1": np.ones(8), "c2": np.full(8, 2.0), "c3": np.full(8, 3.0)}, index=products
+    )
+    agg = pd.DataFrame({"A": np.arange(8.0)}, index=products)
+    result = permutation_column_null(agg, cov, n_iter=3, min_knockoffs=2, random_state=0)
+    assert not result["ok"]
+    assert "singular" in result["reason"]
+
+
 def test_relate_observational_requires_numeric_descriptors():
     obs = pd.DataFrame({"product": PRODUCTS, "grade": list("AABBCCA")})
     result = validate_descriptive(_panel(), obs, mode="observational")