diff --git a/CHANGELOG.md b/CHANGELOG.md index 81094b27..94f66fde 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,23 @@ those changes. ## [Unreleased] +## [1.56.0] - 2026-07-22 + +### Changed + +- The observational sensory relate (`process_improve.sensory.analyze_descriptive`) + is now influence-robust for sparse, wide predictor blocks. Both the marginal + associations and the cross-validated discriminator are gated on a leave-one-out + jackknife, so an association or predictive coefficient that rests on a single + high-leverage observation (a descriptor that is non-zero on only one product) is + demoted instead of reported as a driver. Marginal `significant` now requires both + Benjamini-Hochberg rejection and jackknife robustness; `discriminator_significant` + additionally requires the per-coefficient jackknife interval to exclude zero. The + jackknife reuses the existing `alpha` and the number of observations, adding no new + threshold, so genuine multi-observation drivers are unaffected. Each association + now carries `jackknife_se`, `influence_robust` and `n_supporting`; each + discriminator descriptor now carries `jackknife_significant`. + ## [1.55.1] - 2026-07-21 ### Documentation @@ -2549,7 +2566,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.55.1...HEAD +[Unreleased]: https://github.com/kgdunn/process-improve/compare/v1.56.0...HEAD +[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 [1.55.0]: https://github.com/kgdunn/process-improve/compare/v1.54.0...v1.55.0 [1.54.0]: https://github.com/kgdunn/process-improve/compare/v1.53.0...v1.54.0 diff --git a/CITATION.cff b/CITATION.cff index f0676276..f77374b4 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -12,8 +12,8 @@ authors: repository-code: "https://github.com/kgdunn/process-improve" url: "https://kgdunn.github.io/process-improve/" license: MIT -version: 1.55.1 -date-released: "2026-07-21" +version: 1.56.0 +date-released: "2026-07-22" keywords: - chemometrics - multivariate analysis diff --git a/pyproject.toml b/pyproject.toml index dd4c3fce..d97343f5 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "process-improve" -version = "1.55.1" +version = "1.56.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/analysis.py b/src/process_improve/sensory/analysis.py index 7c86af55..6342d939 100644 --- a/src/process_improve/sensory/analysis.py +++ b/src/process_improve/sensory/analysis.py @@ -17,7 +17,13 @@ The observational relate corrects across the family of tests with Benjamini-Hochberg FDR and returns supporting product means with confidence -intervals and a PCA sensory map. +intervals and a PCA sensory map. Both the marginal associations and the +cross-validated discriminator are additionally gated on a leave-one-out +jackknife, so an association or predictive coefficient that rests on a single +high-leverage observation (a predictor that is non-zero on only one product, +common in sparse, wide descriptor blocks) is demoted rather than reported. The +jackknife adds no threshold of its own: it reuses the same ``alpha`` and the +number of observations, so a genuine multi-observation driver is unaffected. """ from __future__ import annotations @@ -29,6 +35,7 @@ import numpy as np import pandas as pd from scipy.stats import pearsonr +from scipy.stats import t as t_dist from process_improve.multivariate.methods import PCA, PLS, MCUVScaler, selectivity_ratio, vip from process_improve.sensory.mam import MAMResult, align_scores, mixed_assessor_model @@ -98,6 +105,16 @@ def aggregate_to_product(panel: pd.DataFrame) -> pd.DataFrame: return wide +#: Minimum usable observations before the leave-one-out jackknife of an +#: association is defined; below this an association cannot be certified as +#: influence-robust (mirrors the ``len < 4`` guard in ``panel._mad_bands``). +_MIN_OBS_FOR_JACKKNIFE = 4 + +#: Correlations are clipped off +/-1 before the Fisher-z transform so ``arctanh`` +#: stays finite. +_R_CLIP = 1.0 - 1e-12 + + def _attach_fdr(records: list[dict[str, Any]], alpha: float) -> list[dict[str, Any]]: """Attach Benjamini-Hochberg q-values (and reject flags) to ``records``.""" pvals = [r["p_value"] for r in records] @@ -106,10 +123,55 @@ def _attach_fdr(records: list[dict[str, Any]], alpha: float) -> list[dict[str, A bh = benjamini_hochberg(np.asarray(pvals), alpha=alpha) for rec, q, rej in zip(records, bh.p_adjusted, bh.reject, strict=True): rec["q_value"] = float(q) - rec["significant"] = bool(rej) + # Harden the marginal significance: an association counts only when it also + # survives the leave-one-out jackknife, so a single high-leverage + # observation cannot manufacture a "significant" correlation. + rec["significant"] = bool(rej) and bool(rec.get("influence_robust", True)) return records +def _jackknife_correlation(x: np.ndarray, y: np.ndarray, alpha: float) -> tuple[float, bool, int]: + """Leave-one-out jackknife significance of a Pearson correlation. + + Returns ``(jackknife_se, influence_robust, n_supporting)``. The correlation is + Fisher-z transformed, and ``influence_robust`` is ``True`` when the two-sided + ``alpha`` jackknife confidence interval for the transformed correlation excludes + zero, i.e. the association does not collapse when any single observation is + removed. A predictor that is non-zero on a single high-leverage observation has + one deletion that drives the correlation to zero, which inflates the jackknife + standard error until the interval spans zero and the association is demoted. The + rule introduces no threshold of its own: it reuses ``alpha`` and the number of + observations, so a genuine multi-observation driver stays significant while a + single-support spike does not. + """ + n = int(x.size) + if n < _MIN_OBS_FOR_JACKKNIFE: + return float("nan"), False, n + + def _z(r_value: float) -> float: + return float(np.arctanh(np.clip(r_value, -_R_CLIP, _R_CLIP))) + + z_full = _z(float(pearsonr(x, y)[0])) + pseudo = np.empty(n) + for i in range(n): + keep = np.arange(n) != i + xi, yi = x[keep], y[keep] + # Removing the sole support of a spike leaves a constant column: no variance + # means no correlation, so the deleted-sample estimate is zero. + r_i = float(pearsonr(xi, yi)[0]) if xi.std() > 0 and yi.std() > 0 else 0.0 + pseudo[i] = n * z_full - (n - 1) * _z(r_i) + z_mean = float(pseudo.mean()) + # ``se`` is always finite here: the correlations are clipped before the Fisher-z + # transform and the caller guarantees non-constant inputs, so no NaN/inf arises. + se = float(pseudo.std(ddof=1) / np.sqrt(n)) + if se <= 0.0: + # Every leave-one-out deletion gives the same correlation, so it is + # perfectly stable: influence-robust iff that correlation is non-zero. + return se, bool(abs(z_mean) > 0.0), n + t_crit = float(t_dist.ppf(1.0 - alpha / 2.0, df=n - 1)) + return se, bool(abs(z_mean) > t_crit * se), n + + def _fit_pls_safe(x: pd.DataFrame, y: pd.DataFrame, n_components: int) -> tuple[PLS | None, int]: """Fit PLS, stepping the component count down on a near-collinear (singular) block. @@ -221,8 +283,11 @@ def discriminate_observational( # noqa: PLR0913, PLR0915 dict ``per_attribute`` (the Q-squared gate per attribute), ``descriptors`` (the per attribute-descriptor selectivity ratio, permutation q-value, - ``discriminator_significant`` flag and ``cluster_id``), ``clusters`` - (the descriptor-to-cluster map), and the settings used. + ``jackknife_significant`` flag from the leave-one-out beta confidence + interval, ``discriminator_significant`` flag and ``cluster_id``), + ``clusters`` (the descriptor-to-cluster map), and the settings used. A + descriptor is ``discriminator_significant`` only when it also survives the + jackknife, so a coefficient carried by a single product is demoted. """ 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])] @@ -249,17 +314,25 @@ def discriminate_observational( # noqa: PLR0913, PLR0915 pls, a = _fit_pls_safe(x_scaled, y_scaled, a) # 1. Leave-one-out cross-validated Q-squared gate: is the attribute - # predictable from the descriptor block out of sample? + # predictable from the descriptor block out of sample? The same LOO + # refit also yields a jackknife confidence interval per descriptor + # coefficient (Martens' uncertainty test); it is reused below so a + # descriptor whose predictive weight rests on a single high-leverage + # product is demoted even when it survives the permutation null. predictable = False q2_cv = float("nan") rmsep_cv = float("nan") + jack_significant: dict[str, bool] = {} if pls is not None and n_rows >= 5 and y_attr.std() > 0: with warnings.catch_warnings(): warnings.simplefilter("ignore", RuntimeWarning) - cv = pls.cross_validate(x_scaled, y_scaled, cv="loo", show_progress=False) + cv = pls.cross_validate( + x_scaled, y_scaled, cv="loo", conf_level=1.0 - alpha, show_progress=False + ) q2_cv = float(cv.q_squared.iloc[0]) rmsep_cv = float(cv.rmse_cv.iloc[0]) predictable = q2_cv > 0.0 + jack_significant = {str(name): bool(flag) for name, flag in cv.significant.iloc[:, 0].items()} per_attribute.append( { "attribute": str(attr), @@ -309,6 +382,7 @@ def discriminate_observational( # noqa: PLR0913, PLR0915 p_raw = np.ones(k) p_maxt = np.ones(k) for i, desc in enumerate(descriptors): + desc_robust = jack_significant.get(str(desc), False) records.append( { "attribute": str(attr), @@ -316,7 +390,8 @@ def discriminate_observational( # noqa: PLR0913, PLR0915 "selectivity_ratio": float(sr_values[i]), "p_value": float(p_raw[i]), "q_value": float(p_maxt[i]), - "discriminator_significant": bool(p_maxt[i] <= alpha and predictable), + "jackknife_significant": bool(desc_robust), + "discriminator_significant": bool(p_maxt[i] <= alpha and predictable and desc_robust), "cluster_id": clusters[str(desc)], } ) @@ -367,7 +442,14 @@ def relate_observational( # noqa: PLR0913 n_permutations: int = 199, random_state: int = 0, ) -> dict[str, Any]: - """Relate the attribute block to measured descriptors with PLS plus correlations.""" + """Relate the attribute block to measured descriptors with PLS plus correlations. + + Each marginal association carries a Pearson ``r``, an FDR ``q_value`` and, from a + leave-one-out jackknife, ``jackknife_se``, ``influence_robust`` and + ``n_supporting``. ``significant`` requires both FDR rejection and jackknife + robustness, so a correlation created by a single high-leverage observation is not + reported as significant. + """ x_block = covariates.loc[agg.index].astype(float) y_block = agg.astype(float) max_comp = max(1, min(n_components, x_block.shape[1], x_block.shape[0] - 1)) @@ -386,9 +468,20 @@ def relate_observational( # noqa: PLR0913 for desc in x_block.columns: pair = pd.concat([y_block[attr], x_block[desc]], axis=1).dropna() if pair.shape[0] >= 3 and pair.iloc[:, 0].std() > 0 and pair.iloc[:, 1].std() > 0: - r, p = pearsonr(pair.iloc[:, 0], pair.iloc[:, 1]) + yv = pair.iloc[:, 0].to_numpy(dtype=float) + xv = pair.iloc[:, 1].to_numpy(dtype=float) + r, p = pearsonr(yv, xv) + jack_se, robust, n_support = _jackknife_correlation(xv, yv, alpha) assoc.append( - {"attribute": str(attr), "descriptor": str(desc), "r": float(r), "p_value": float(p)} + { + "attribute": str(attr), + "descriptor": str(desc), + "r": float(r), + "p_value": float(p), + "jackknife_se": float(jack_se), + "influence_robust": bool(robust), + "n_supporting": int(n_support), + } ) assoc = _attach_fdr(assoc, alpha) result: dict[str, Any] = { diff --git a/tests/test_sensory.py b/tests/test_sensory.py index dbd8cdce..d543dd97 100644 --- a/tests/test_sensory.py +++ b/tests/test_sensory.py @@ -18,7 +18,12 @@ panel_scorecard, validate_descriptive, ) -from process_improve.sensory.analysis import _collinear_clusters, discriminate_observational, relate_designed +from process_improve.sensory.analysis import ( + _collinear_clusters, + _jackknife_correlation, + discriminate_observational, + relate_designed, +) from process_improve.sensory.ingest import reshape_to_long from process_improve.univariate.metrics import benjamini_hochberg @@ -74,6 +79,54 @@ def _panel(*, anomalous: str | None = "P8", seed: int = 0) -> pd.DataFrame: return pd.DataFrame(rows) +# Products for the high-leverage case study. ``Q5`` is both a response outlier on +# attribute C and the sole product carrying the ``spike`` descriptor. +LEVERAGE_PRODUCTS = [f"Q{i}" for i in range(12)] +_SPIKE_INDEX = 5 + + +def _leverage_case(seed: int = 1) -> tuple[pd.DataFrame, pd.DataFrame]: + """Panel + covariates isolating a genuine driver from a single-support spike. + + ``genuine`` varies smoothly across all products and drives attribute A. + ``spike`` is zero on every product except one (``Q5``); attribute C is flat + except at that same product, where it is a large outlier. The spike-vs-C + correlation is therefore created entirely by one high-leverage product: strong + in-sample, but it collapses the moment that product is removed. An + influence-robust relate must keep genuine-on-A and demote spike-on-C. + """ + rng = np.random.default_rng(seed) + n = len(LEVERAGE_PRODUCTS) + genuine = np.linspace(0.2, 1.0, n) + spike = np.zeros(n) + spike[_SPIKE_INDEX] = 1.0 + cov = pd.DataFrame( + { + "product": LEVERAGE_PRODUCTS, + "genuine": genuine, + "spike": spike, + "noise": rng.normal(0.0, 1.0, n), + } + ) + a_mean = 4.0 * genuine + c_mean = np.zeros(n) + c_mean[_SPIKE_INDEX] = 5.0 # one product is a large outlier on attribute C + rows = [] + for pid in [f"J{i}" for i in range(6)]: + bias = rng.normal(0.0, 0.2) + for j, prod in enumerate(LEVERAGE_PRODUCTS): + for rep in (1, 2): + rows.append( + {"panelist_id": pid, "session": 1, "product": prod, "attribute": "A", + "replicate": rep, "score": 5.0 + a_mean[j] + bias + rng.normal(0, 0.15)} + ) + rows.append( + {"panelist_id": pid, "session": 1, "product": prod, "attribute": "C", + "replicate": rep, "score": 5.0 + c_mean[j] + bias + rng.normal(0, 0.15)} + ) + return pd.DataFrame(rows), cov + + # --------------------------------------------------------------------------- # Validation # --------------------------------------------------------------------------- @@ -300,6 +353,121 @@ def test_discriminator_gate_and_clusters(): assert not desc[desc["descriptor"] == "d3"]["discriminator_significant"].any() +def test_jackknife_correlation_edge_cases(): + """The jackknife helper's degenerate branches behave sensibly.""" + rng = np.random.default_rng(0) + + # Too few observations: the jackknife is undefined, so not robust. + se, robust, n = _jackknife_correlation(np.arange(3.0), np.arange(3.0) + 1.0, 0.05) + assert n == 3 + assert not robust + assert np.isnan(se) + + # A smooth relationship supported across all points is robust. + x = np.linspace(0.0, 1.0, 12) + y = 2.0 * x + rng.normal(0, 0.05, 12) + se, robust, n = _jackknife_correlation(x, y, 0.05) + assert robust + assert se > 0.0 + + # A single-support spike collapses when its one point is removed: not robust. + xs = np.zeros(12) + xs[5] = 1.0 + ys = rng.normal(0, 1, 12) + ys[5] += 5.0 + _, robust, _ = _jackknife_correlation(xs, ys, 0.05) + assert not robust + + # A perfectly stable correlation (zero jackknife spread) is robust when non-zero. + z = np.linspace(0.0, 1.0, 8) + se, robust, _ = _jackknife_correlation(z, z, 0.05) + assert se == 0.0 + assert robust + + +def test_relate_marginal_demotes_single_support_spike(): + """A single high-leverage observation must not create a significant association. + + The ``spike`` descriptor is non-zero on exactly one product, which is also a + large outlier on attribute C, so its in-sample Pearson correlation is strong and + would pass a bare FDR gate. The leave-one-out jackknife demotes it (removing that + one product collapses the correlation), while the genuine multi-product driver on + attribute A survives. + """ + panel, cov = _leverage_case() + validated = validate_descriptive(panel, cov, mode="observational") + result = analyze_descriptive(validated, discriminator=False) + assoc = pd.DataFrame(result.relate["associations"]) + + # New influence-robustness fields are surfaced per association. + assert {"jackknife_se", "influence_robust", "n_supporting"}.issubset(assoc.columns) + + spike_c = assoc[(assoc["attribute"] == "C") & (assoc["descriptor"] == "spike")].iloc[0] + genuine_a = assoc[(assoc["attribute"] == "A") & (assoc["descriptor"] == "genuine")].iloc[0] + + # The spike correlates strongly in-sample yet rests on one product: not robust, + # so not significant, even though its raw p-value is tiny. + assert abs(spike_c["r"]) > 0.8 + assert spike_c["p_value"] < 0.05 + assert not spike_c["influence_robust"] + assert not spike_c["significant"] + + # The genuine driver is supported across all products: robust and significant. + assert genuine_a["influence_robust"] + assert genuine_a["significant"] + + +def test_discriminator_demotes_single_support_spike(): + """The cross-validated discriminator must also demote a single-support spike.""" + n = len(LEVERAGE_PRODUCTS) + rng = np.random.default_rng(2) + genuine = np.linspace(0.0, 1.0, n) + spike = np.zeros(n) + spike[_SPIKE_INDEX] = 1.0 + agg = pd.DataFrame( + { + "A": 4.0 * genuine + rng.normal(0, 0.1, n), + "C": np.where(np.arange(n) == _SPIKE_INDEX, 5.0, 0.0) + rng.normal(0, 0.1, n), + }, + index=LEVERAGE_PRODUCTS, + ) + cov = pd.DataFrame( + {"genuine": genuine, "spike": spike, "noise": rng.normal(0, 1, n)}, index=LEVERAGE_PRODUCTS + ) + disc = discriminate_observational(agg, cov, n_components=1, n_permutations=99, random_state=0) + desc = pd.DataFrame(disc["descriptors"]) + + assert "jackknife_significant" in desc.columns + # The spike's predictive weight rests on one product, so its jackknife interval + # spans zero and it is never confirmed. + spike_rows = desc[desc["descriptor"] == "spike"] + assert not spike_rows["jackknife_significant"].any() + assert not spike_rows["discriminator_significant"].any() + # The genuine driver stays jackknife-stable on the attribute it drives, so the + # jackknife gate demotes the single-support spike without over-pruning a real + # driver's predictive coefficient. + genuine_a = desc[(desc["descriptor"] == "genuine") & (desc["attribute"] == "A")].iloc[0] + assert genuine_a["jackknife_significant"] + + +def test_discriminator_small_sample_skips_jackknife_gate(): + """With too few products the LOO gate is skipped, so nothing is confirmed. + + Below the cross-validation floor the per-coefficient jackknife cannot be run, so + ``jackknife_significant`` defaults to ``False`` and no descriptor is confirmed. + """ + products = [f"P{i}" for i in range(4)] + rng = np.random.default_rng(0) + u = np.linspace(0.0, 1.0, 4) + agg = pd.DataFrame({"A": 2.0 * u + rng.normal(0, 0.05, 4)}, index=products) + cov = pd.DataFrame({"d1": u, "d2": rng.normal(0, 1, 4)}, index=products) + + disc = discriminate_observational(agg, cov, n_components=1, n_permutations=19, random_state=0) + desc = pd.DataFrame(disc["descriptors"]) + assert not desc["jackknife_significant"].any() + assert not desc["discriminator_significant"].any() + + def test_relate_observational_requires_numeric_descriptors(): obs = pd.DataFrame({"product": PRODUCTS, "grade": list("AABBCCA")}) result = validate_descriptive(_panel(), obs, mode="observational")