diff --git a/CHANGELOG.md b/CHANGELOG.md index 94f66fd..fa7d4f3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,21 @@ those changes. ## [Unreleased] +## [1.57.0] - 2026-07-23 + +### Added + +- `analyze_descriptive` / `relate_observational` gain an `influence_deletions` + parameter (default 1) for the observational relate. It sets how many observations + the marginal-association jackknife removes together: the default is the ordinary + leave-one-out gate, and raising it to 2 also demotes a correlation carried by a + single pair of high-leverage observations, which leave-one-out cannot detect + (deleting either point of the pair leaves the other holding the correlation up). + For `d >= 2` the decision uses a breakdown criterion (the correlation must stay + significant, with the same sign, after removing every subset of `d` observations) + rather than the averaged jackknife variance, which would dilute the single + collapsing subset. `d = 1` behaviour is unchanged. + ## [1.56.0] - 2026-07-22 ### Changed @@ -2566,7 +2581,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.56.0...HEAD +[Unreleased]: https://github.com/kgdunn/process-improve/compare/v1.57.0...HEAD +[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 [1.55.0]: https://github.com/kgdunn/process-improve/compare/v1.54.0...v1.55.0 diff --git a/CITATION.cff b/CITATION.cff index f77374b..87359ff 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.56.0 -date-released: "2026-07-22" +version: 1.57.0 +date-released: "2026-07-23" keywords: - chemometrics - multivariate analysis diff --git a/pyproject.toml b/pyproject.toml index d97343f..29fbbb4 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "process-improve" -version = "1.56.0" +version = "1.57.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 6342d93..91fce76 100644 --- a/src/process_improve/sensory/analysis.py +++ b/src/process_improve/sensory/analysis.py @@ -28,6 +28,7 @@ from __future__ import annotations +import itertools import warnings from dataclasses import dataclass, field from typing import Any @@ -130,46 +131,110 @@ def _attach_fdr(records: list[dict[str, Any]], alpha: float) -> list[dict[str, A 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. +def _survives_all_deletions(x: np.ndarray, y: np.ndarray, alpha: float, d: int, sign: float) -> bool: + """Return whether the correlation stays significant after removing any ``d`` points. + + The breakdown criterion behind the ``max_deletions >= 2`` path of + :func:`_jackknife_correlation`: over every subset of ``d`` observations removed at + once, the remaining correlation must keep the same sign and stay significant at + ``alpha`` (a Pearson-``t`` test on the ``n - d`` retained points). A subset that + removes the whole support of a spike leaves a constant column and fails at once. + """ + n = int(x.size) + m = n - d + t_crit = float(t_dist.ppf(1.0 - alpha / 2.0, df=m - 2)) + idx = np.arange(n) + for drop in itertools.combinations(range(n), d): + keep = np.isin(idx, drop, invert=True) + xi, yi = x[keep], y[keep] + if xi.std() <= 0 or yi.std() <= 0: + return False # deletion removed the whole support of the effect + r_s = float(pearsonr(xi, yi)[0]) + if np.sign(np.arctanh(np.clip(r_s, -_R_CLIP, _R_CLIP))) != sign: + return False # deletion flipped the direction of the effect + t_s = abs(r_s) * np.sqrt((m - 2) / max(1.0 - r_s**2, 1e-12)) + if t_s <= t_crit: + return False # deletion made the remaining correlation non-significant + return True + + +def _jackknife_correlation( + x: np.ndarray, y: np.ndarray, alpha: float, *, max_deletions: int = 1 +) -> tuple[float, bool, int]: + """Return the delete-``d`` 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. + Fisher-z transformed. ``jackknife_se`` is always the ordinary leave-one-out (Tukey) + jackknife standard error. ``influence_robust`` says whether the association survives + removing any ``d = max_deletions`` observations: + + * ``d = 1`` (default): the leave-one-out jackknife confidence interval for the + correlation excludes zero. A predictor non-zero on a single observation has one + deletion that drives the correlation to zero, inflating the standard error until + the interval spans zero, so the pair is demoted. + * ``d >= 2``: leave-one-out is blind to an effect carried by two observations, since + deleting either one leaves the other holding the correlation up. Use the breakdown + criterion instead - the correlation must stay significant, with the same sign, + after removing *every* subset of ``d`` observations (the worst case, not the + averaged jackknife variance, which would dilute the single collapsing subset). + + The rule adds no threshold of its own: it reuses ``alpha``, the number of + observations, and ``d``, so a driver supported by more than ``d`` observations stays + significant while one carried by ``d`` or fewer does not. + + Parameters + ---------- + x, y : numpy.ndarray + Paired, non-constant observation vectors (the caller guarantees non-zero + variance). + alpha : float + Two-sided significance level for the jackknife confidence interval. + max_deletions : int + Number ``d`` of observations removed together in each jackknife subset. Must + be at least 1; the effective cost is ``comb(n, d)`` correlation refits. """ + if max_deletions < 1: + raise ValueError(f"max_deletions must be at least 1, got {max_deletions}.") n = int(x.size) - if n < _MIN_OBS_FOR_JACKKNIFE: + d = max_deletions + # Need at least _MIN_OBS_FOR_JACKKNIFE observations, and enough left after the + # deletion to still estimate a correlation. + if n < _MIN_OBS_FOR_JACKKNIFE or n - d < _MIN_OBS_FOR_JACKKNIFE - 1: return float("nan"), False, n def _z(r_value: float) -> float: return float(np.arctanh(np.clip(r_value, -_R_CLIP, _R_CLIP))) + # Reported ``jackknife_se`` is always the ordinary leave-one-out (Tukey) jackknife + # standard error of the Fisher-z correlation: a stable, interpretable influence + # magnitude independent of ``d``. 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 + + if d == 1: + # Leave-one-out: the shipped jackknife-t decision, unchanged. A single-support + # spike has one deletion that drives the correlation to zero, inflating ``se`` + # until the confidence interval spans zero. + if se <= 0.0: + # Every deletion gives the same correlation: perfectly stable, 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 + + # d >= 2: a correlation carried by a single pair of observations survives every + # *single* deletion, so the averaged jackknife variance stays small and would not + # demote it. Use the breakdown criterion instead (see ``_survives_all_deletions``): + # the correlation must stay significant, with the same sign, after removing *every* + # subset of ``d`` observations. + return se, _survives_all_deletions(x, y, alpha, d, np.sign(z_full)), n def _fit_pls_safe(x: pd.DataFrame, y: pd.DataFrame, n_components: int) -> tuple[PLS | None, int]: @@ -441,14 +506,17 @@ def relate_observational( # noqa: PLR0913 discriminator: bool = True, n_permutations: int = 199, random_state: int = 0, + influence_deletions: int = 1, ) -> dict[str, Any]: """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. + delete-``influence_deletions`` jackknife, ``jackknife_se``, ``influence_robust`` + and ``n_supporting``. ``significant`` requires both FDR rejection and jackknife + robustness, so a correlation created by too few high-leverage observations is not + reported as significant. ``influence_deletions`` (default 1, ordinary leave-one-out) + sets how many observations are removed together: raise it to 2 to also demote a + correlation carried by a single pair of observations. """ x_block = covariates.loc[agg.index].astype(float) y_block = agg.astype(float) @@ -471,7 +539,9 @@ def relate_observational( # noqa: PLR0913 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) + jack_se, robust, n_support = _jackknife_correlation( + xv, yv, alpha, max_deletions=influence_deletions + ) assoc.append( { "attribute": str(attr), @@ -544,6 +614,7 @@ def analyze_descriptive( # noqa: PLR0913 discriminator: bool = True, n_permutations: int = 199, random_state: int = 0, + influence_deletions: int = 1, ) -> AnalysisResult: """Run the descriptive pipeline: panel check, correction, and relate. @@ -579,6 +650,11 @@ def analyze_descriptive( # noqa: PLR0913 Permutations for the discriminator's selectivity-ratio null. random_state : int Seed for the discriminator's permutations and cross-validation folds. + influence_deletions : int + How many observations the marginal-association jackknife removes together + (default 1, ordinary leave-one-out). Raising it to 2 also demotes a + correlation carried by a single pair of high-leverage observations, which + leave-one-out cannot detect. Returns ------- @@ -623,6 +699,7 @@ def analyze_descriptive( # noqa: PLR0913 discriminator=discriminator, n_permutations=n_permutations, random_state=random_state, + influence_deletions=influence_deletions, ) return AnalysisResult( @@ -644,6 +721,7 @@ def analyze_descriptive( # noqa: PLR0913 "discriminator": discriminator, "n_permutations": n_permutations, "random_state": random_state, + "influence_deletions": influence_deletions, "content_hash": validated.content_hash, }, ) diff --git a/tests/test_sensory.py b/tests/test_sensory.py index d543dd9..53a16af 100644 --- a/tests/test_sensory.py +++ b/tests/test_sensory.py @@ -385,6 +385,90 @@ def test_jackknife_correlation_edge_cases(): assert robust +def test_jackknife_delete_two_demotes_two_support_spike(): + """Leave-one-out keeps a two-point spike; delete-two removes it, genuine survives.""" + rng = np.random.default_rng(0) + n = 15 + # A predictor non-zero on exactly two products, both extreme on the response. + xs = np.zeros(n) + xs[13] = xs[14] = 1.0 + ys = rng.normal(0, 1, n) + ys[13] += 6.0 + ys[14] += 6.0 + _, robust_d1, _ = _jackknife_correlation(xs, ys, 0.05, max_deletions=1) + _, robust_d2, _ = _jackknife_correlation(xs, ys, 0.05, max_deletions=2) + assert robust_d1 # leave-one-out is blind to a two-point effect + assert not robust_d2 # removing both supporting points collapses it + + # A genuine broad driver survives both. + xg = np.linspace(0.0, 1.0, n) + yg = 2.0 * xg + rng.normal(0, 0.2, n) + _, g1, _ = _jackknife_correlation(xg, yg, 0.05, max_deletions=1) + _, g2, _ = _jackknife_correlation(xg, yg, 0.05, max_deletions=2) + assert g1 + assert g2 + + with pytest.raises(ValueError, match="max_deletions must be at least 1"): + _jackknife_correlation(xg, yg, 0.05, max_deletions=0) + + +def test_jackknife_breakdown_demotes_via_nonsignificance_and_sign_flip(): + """Delete-two demotes through the non-significance and sign-flip paths, not only support removal.""" + # Moderate continuous correlation (no zeros, so no deletion makes a constant column): + # robust to any single deletion, but removing the two most influential points drops + # the remaining correlation below significance while keeping its sign. + x = np.array([-1.1, -0.81, -0.78, -0.75, -0.73, -0.45, -0.25, 0.13, 0.27, 0.48, 0.84, 0.86]) + y = np.array([-1.69, -0.9, -2.43, -2.68, -0.21, -3.6, -0.55, 0.85, -1.41, -1.42, -0.45, 1.96]) + _, a1, _ = _jackknife_correlation(x, y, 0.05, max_deletions=1) + _, a2, _ = _jackknife_correlation(x, y, 0.05, max_deletions=2) + assert a1 + assert not a2 + + # Two extreme points create a positive correlation; removing them leaves a negative + # trend, so the effect reverses direction under a two-point deletion. + xb = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 20, 21.0]) + yb = np.array([10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 30, 31.0]) + _, b1, _ = _jackknife_correlation(xb, yb, 0.05, max_deletions=1) + _, b2, _ = _jackknife_correlation(xb, yb, 0.05, max_deletions=2) + assert b1 + assert not b2 + + +def test_relate_influence_deletions_two_demotes_two_support_spike(): + """`influence_deletions=2` demotes a two-product spike the default gate keeps.""" + n = len(LEVERAGE_PRODUCTS) + rng = np.random.default_rng(3) + genuine = np.linspace(0.2, 1.0, n) + two_spike = np.zeros(n) + two_spike[_SPIKE_INDEX] = two_spike[_SPIKE_INDEX + 1] = 1.0 + + a_mean = 4.0 * genuine + # Attribute C is flat except at the two products that carry the spike. + c_mean = np.zeros(n) + c_mean[_SPIKE_INDEX] = c_mean[_SPIKE_INDEX + 1] = 5.0 + 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)}) + panel = pd.DataFrame(rows) + cov = pd.DataFrame({"product": LEVERAGE_PRODUCTS, "genuine": genuine, "two_spike": two_spike}) + validated = validate_descriptive(panel, cov, mode="observational") + + def spike_robust(deletions: int) -> bool: + result = analyze_descriptive(validated, discriminator=False, influence_deletions=deletions) + assoc = pd.DataFrame(result.relate["associations"]) + row = assoc[(assoc["attribute"] == "C") & (assoc["descriptor"] == "two_spike")].iloc[0] + return bool(row["influence_robust"]) + + assert spike_robust(1) # default leave-one-out keeps the two-support spike + assert not spike_robust(2) # delete-two demotes it + + def test_relate_marginal_demotes_single_support_spike(): """A single high-leverage observation must not create a significant association.