Add delete-d option to the observational relate influence gate

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

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

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"

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,
         },
     )

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.