diff --git a/CHANGELOG.md b/CHANGELOG.md index fa7d4f3..a4a9826 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/CITATION.cff b/CITATION.cff index 87359ff..a544588 100644 --- a/CITATION.cff +++ b/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 diff --git a/pyproject.toml b/pyproject.toml index 29fbbb4..3996999 100644 --- a/pyproject.toml +++ b/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" diff --git a/src/process_improve/sensory/__init__.py b/src/process_improve/sensory/__init__.py index 5eb28b3..63443bb 100644 --- a/src/process_improve/sensory/__init__.py +++ b/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", diff --git a/src/process_improve/sensory/analysis.py b/src/process_improve/sensory/analysis.py index 91fce76..4e821d1 100644 --- a/src/process_improve/sensory/analysis.py +++ b/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]] = [] diff --git a/tests/test_sensory.py b/tests/test_sensory.py index 53a16af..35af93e 100644 --- a/tests/test_sensory.py +++ b/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")