Add per-panelist consistency test via own-score permutation

CHANGELOG.md

@@ -11,6 +11,23 @@ those changes.
 
 ## [Unreleased]
 
+## [1.59.0] - 2026-07-23
+
+### Added
+
+- `process_improve.sensory.panel_consistency`: a per-panelist reliability test that
+  compares each panelist's product signal against a reshuffle of their own scores. For
+  every panelist the product labels are permuted within each attribute (keeping the
+  panelist's own values, only reassigning which product got which) and the mean
+  eta-squared product effect and the agreement with the leave-one-out panel consensus are
+  recomputed; the null is built entirely from the panelist's own data, no simulation. A
+  panelist who carries a real, reproducible signal beats almost every reshuffle (small
+  p-value), while a noise rater is indistinguishable from its own permutations. Returns a
+  `PanelConsistency` with per-panelist `discrimination` / `p_discrimination`, `agreement`
+  / `p_agreement`, and a `consistent` flag (gated on the discrimination p-value; agreement
+  is reported but not gated, since a reproducible rater may still rank the products
+  differently from the panel).
+
 ## [1.58.0] - 2026-07-23
 
 ### Added
@@ -2599,7 +2616,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.58.0...HEAD
+[Unreleased]: https://github.com/kgdunn/process-improve/compare/v1.59.0...HEAD
+[1.59.0]: https://github.com/kgdunn/process-improve/compare/v1.58.0...v1.59.0
 [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

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.58.0
+version: 1.59.0
 date-released: "2026-07-23"
 keywords:
   - chemometrics

pyproject.toml

@@ -1,6 +1,6 @@
 [project]
 name = "process-improve"
-version = "1.58.0"
+version = "1.59.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

@@ -27,7 +27,13 @@
 )
 from process_improve.sensory.ingest import reshape_to_long
 from process_improve.sensory.mam import MAMResult, align_scores, mixed_assessor_model
-from process_improve.sensory.panel import PanelScorecard, apply_correction, panel_scorecard
+from process_improve.sensory.panel import (
+    PanelConsistency,
+    PanelScorecard,
+    apply_correction,
+    panel_consistency,
+    panel_scorecard,
+)
 from process_improve.sensory.recipes import SENSORY_RECIPES
 from process_improve.sensory.validation import (
     DESCRIPTIVE_LONG_COLUMNS,
@@ -40,6 +46,7 @@
     "SENSORY_RECIPES",
     "AnalysisResult",
     "MAMResult",
+    "PanelConsistency",
     "PanelScorecard",
     "ValidationResult",
     "aggregate_to_product",
@@ -48,6 +55,7 @@
     "apply_correction",
     "discriminate_observational",
     "mixed_assessor_model",
+    "panel_consistency",
     "panel_scorecard",
     "permutation_column_null",
     "product_means",

src/process_improve/sensory/panel.py

@@ -197,6 +197,213 @@ def panel_scorecard(panel: pd.DataFrame) -> PanelScorecard:
     return PanelScorecard(table=table, flagged=sorted(reasons), reasons=reasons)
 
 
+#: Minimum products before a panelist's own-permutation null is meaningful.
+_MIN_PRODUCTS_FOR_PERMUTATION = 3
+
+
+@dataclass
+class PanelConsistency:
+    """Outcome of :func:`panel_consistency`.
+
+    Attributes
+    ----------
+    table : pandas.DataFrame
+        One row per panelist, indexed by ``panelist_id``, with the columns
+        ``discrimination``, ``p_discrimination``, ``agreement``, ``p_agreement``,
+        and ``consistent``.
+    inconsistent : list of str
+        Panelist ids whose product signal is not distinguishable from a reshuffle of
+        their own scores (``consistent`` is ``False``).
+    n_permutations : int
+        Number of within-panelist permutations used for each null.
+    """
+
+    table: pd.DataFrame
+    inconsistent: list[str]
+    n_permutations: int
+
+
+def _eta_squared_codes(scores: np.ndarray, codes: np.ndarray, n_groups: int) -> float:
+    """Return the eta-squared of an integer-``codes`` grouping of ``scores``."""
+    grand = float(scores.mean())
+    ss_total = float(((scores - grand) ** 2).sum())
+    if ss_total <= 0.0:
+        return float("nan")
+    sums = np.bincount(codes, weights=scores, minlength=n_groups)
+    counts = np.bincount(codes, minlength=n_groups).astype(float)
+    means = np.divide(sums, counts, out=np.full(n_groups, grand), where=counts > 0)
+    ss_between = float((counts * (means - grand) ** 2).sum())
+    return ss_between / ss_total
+
+
+def _group_means(scores: np.ndarray, codes: np.ndarray, n_groups: int) -> np.ndarray:
+    """Return the per-group mean of ``scores`` (NaN for an empty group)."""
+    sums = np.bincount(codes, weights=scores, minlength=n_groups)
+    counts = np.bincount(codes, minlength=n_groups).astype(float)
+    return np.divide(sums, counts, out=np.full(n_groups, np.nan), where=counts > 0)
+
+
+def _correlation(a: np.ndarray, b: np.ndarray) -> float:
+    """Return the Pearson correlation of two aligned vectors, or NaN if undefined."""
+    if a.size < 2 or a.std() <= 0.0 or b.std() <= 0.0:
+        return float("nan")
+    return float(np.corrcoef(a, b)[0, 1])
+
+
+def _consistency_pvalues(  # noqa: PLR0913
+    blocks: list[tuple[str, np.ndarray, np.ndarray, int]],
+    locators: list[tuple[int, int]],
+    consensus: np.ndarray,
+    real: tuple[float, float],
+    rng: np.random.Generator,
+    n_permutations: int,
+) -> tuple[float, float]:
+    """Return one-sided permutation p-values for discrimination and agreement.
+
+    Each permutation reshuffles the product labels within every attribute (breaking the
+    panelist's product-to-score link while keeping their own score multiset), then
+    recomputes the mean eta-squared and the agreement with the fixed leave-one-out
+    consensus. The p-value is the add-one fraction of permutations reaching at least the
+    real statistic (higher is better for both).
+    """
+    disc_real, agree_real = real
+    ge_disc = 0
+    ge_agree = 0
+    for _ in range(n_permutations):
+        # Permute the product labels once per attribute block and reuse that single
+        # reshuffle for both the eta-squared and the agreement statistic.
+        etas = []
+        means_by_block = []
+        for _attr, scores, codes, ng in blocks:
+            pcodes = rng.permutation(codes)
+            etas.append(_eta_squared_codes(scores, pcodes, ng))
+            means_by_block.append(_group_means(scores, pcodes, ng))
+        null_disc = float(np.nanmean(etas)) if etas else float("nan")
+        if not np.isnan(null_disc) and null_disc >= disc_real:
+            ge_disc += 1
+        if not np.isnan(agree_real):
+            null_own = np.array([means_by_block[bi][gi] for bi, gi in locators])
+            null_agree = _correlation(null_own, consensus)
+            if not np.isnan(null_agree) and null_agree >= agree_real:
+                ge_agree += 1
+    p_disc = (ge_disc + 1.0) / (n_permutations + 1.0)
+    p_agree = (ge_agree + 1.0) / (n_permutations + 1.0) if not np.isnan(agree_real) else float("nan")
+    return p_disc, p_agree
+
+
+def panel_consistency(
+    panel: pd.DataFrame,
+    *,
+    n_permutations: int = 299,
+    random_state: int = 0,
+    alpha: float = 0.05,
+) -> PanelConsistency:
+    """Test each panelist's product signal against a reshuffle of their own scores.
+
+    For every panelist the product labels are permuted within each attribute (keeping the
+    panelist's own score values, only reassigning which product got which), and the mean
+    eta-squared product effect (:func:`_eta_squared_codes`) and the agreement with the
+    leave-one-out panel consensus are recomputed. A panelist who carries a real,
+    reproducible product signal beats almost every reshuffle of their own scores (small
+    p-value); one whose scores are noise is indistinguishable from their own permutations
+    (p near 1). The null is built entirely from the panelist's own data - no simulated or
+    parametric reference.
+
+    ``consistent`` gates on the discrimination p-value (does the panelist separate the
+    products beyond their own noise). The agreement p-value is reported alongside but not
+    gated: a panelist can discriminate reproducibly yet rank the products differently from
+    the panel, which is disagreement rather than inconsistency.
+
+    Parameters
+    ----------
+    panel : pandas.DataFrame
+        Validated ``descriptive_long`` panel data.
+    n_permutations : int
+        Within-panelist permutations per null.
+    random_state : int
+        Seed for the permutations.
+    alpha : float
+        Significance level for the ``consistent`` flag.
+
+    Returns
+    -------
+    PanelConsistency
+        See the class docstring.
+
+    Examples
+    --------
+    >>> result = panel_consistency(validated.normalized_df, n_permutations=299)
+    >>> result.inconsistent
+    ['P8']
+    """
+    if n_permutations < 1:
+        raise ValueError(f"n_permutations must be at least 1, got {n_permutations}.")
+
+    cell = panel.groupby(["panelist_id", "product", "attribute"], observed=True)["score"].mean().reset_index()
+    rng = np.random.default_rng(random_state)
+    records: dict[str, dict[str, float | bool]] = {}
+    for pid in panel["panelist_id"].unique():
+        pan = panel[panel["panelist_id"] == pid]
+        blocks: list[tuple[str, np.ndarray, np.ndarray, int]] = []
+        block_labels: list[np.ndarray] = []
+        etas_real: list[float] = []
+        for attr, grp in pan.groupby("attribute", observed=True):
+            scores = grp["score"].to_numpy(dtype=float)
+            valid = ~np.isnan(scores)
+            labels, codes = np.unique(grp["product"].to_numpy()[valid], return_inverse=True)
+            if valid.sum() < 2 or labels.size < 2:
+                continue
+            blocks.append((str(attr), scores[valid], codes, labels.size))
+            block_labels.append(labels)
+            etas_real.append(_eta_squared_codes(scores[valid], codes, labels.size))
+
+        # Agreement against the leave-one-out consensus, over the (product, attribute)
+        # cells this panelist and the rest of the panel share.
+        loo = cell[cell["panelist_id"] != pid].groupby(["product", "attribute"], observed=True)["score"].mean()
+        loo_map = {(str(p), str(a)): float(v) for (p, a), v in loo.items()}
+        locators: list[tuple[int, int]] = []
+        consensus: list[float] = []
+        real_own: list[float] = []
+        for bi, (attr, scores, codes, ng) in enumerate(blocks):
+            means = _group_means(scores, codes, ng)
+            labels = block_labels[bi]
+            for gi in range(ng):
+                key = (str(labels[gi]), attr)
+                if key in loo_map and not np.isnan(means[gi]):
+                    locators.append((bi, gi))
+                    consensus.append(loo_map[key])
+                    real_own.append(float(means[gi]))
+        finite_etas = [e for e in etas_real if not np.isnan(e)]
+        discrimination = float(np.mean(finite_etas)) if finite_etas else float("nan")
+        agreement = _correlation(np.array(real_own), np.array(consensus)) if len(real_own) >= 2 else float("nan")
+
+        if not blocks or pan["product"].nunique() < _MIN_PRODUCTS_FOR_PERMUTATION or np.isnan(discrimination):
+            records[str(pid)] = {
+                "discrimination": discrimination,
+                "p_discrimination": float("nan"),
+                "agreement": agreement,
+                "p_agreement": float("nan"),
+                "consistent": False,
+            }
+            continue
+
+        p_disc, p_agree = _consistency_pvalues(
+            blocks, locators, np.array(consensus), (discrimination, agreement), rng, n_permutations
+        )
+        records[str(pid)] = {
+            "discrimination": discrimination,
+            "p_discrimination": p_disc,
+            "agreement": agreement,
+            "p_agreement": p_agree,
+            "consistent": bool(p_disc <= alpha),
+        }
+
+    table = pd.DataFrame.from_dict(records, orient="index")
+    table.index.name = "panelist_id"
+    inconsistent = sorted(str(pid) for pid in table.index if not bool(table.loc[pid, "consistent"]))
+    return PanelConsistency(table=table, inconsistent=inconsistent, n_permutations=n_permutations)
+
+
 def apply_correction(panel: pd.DataFrame, drop: list[str]) -> pd.DataFrame:
     """Return ``panel`` with the listed panelists removed.
 

tests/test_sensory.py

@@ -15,6 +15,7 @@
     align_scores,
     analyze_descriptive,
     mixed_assessor_model,
+    panel_consistency,
     panel_scorecard,
     validate_descriptive,
 )
@@ -271,6 +272,77 @@ def test_scorecard_clean_panel_has_no_flags():
     assert card.flagged == []
 
 
+def test_consistency_flags_the_noise_rater():
+    """A panelist scoring at random is indistinguishable from a reshuffle of its own scores."""
+    result = panel_consistency(_panel(anomalous="P8"), n_permutations=199, random_state=0)
+    table = result.table
+    assert {"discrimination", "p_discrimination", "agreement", "p_agreement", "consistent"} == set(table.columns)
+
+    # The random rater cannot beat its own permutations: high discrimination p-value,
+    # not consistent, and the least consistent of the whole panel.
+    assert "P8" in result.inconsistent
+    assert table.loc["P8", "p_discrimination"] > 0.5
+    assert table["p_discrimination"].idxmax() == "P8"
+    assert table.loc["P8", "agreement"] < 0.5
+
+    # The strongest genuine raters clear their own-permutation null.
+    assert table.loc["P5", "consistent"]
+    assert table.loc["P4", "consistent"]
+
+
+def _row(panelist: str, product: str, attribute: str, replicate: int, score: float) -> dict:
+    """Build a single descriptive_long row."""
+    return {
+        "panelist_id": panelist,
+        "session": 1,
+        "product": product,
+        "attribute": attribute,
+        "replicate": replicate,
+        "score": score,
+    }
+
+
+def test_consistency_handles_degenerate_panelists():
+    """A flat rater and a single-product rater cannot be assessed and are not consistent."""
+    products = list("WXYZ")
+    truth = {"W": -1, "X": 0, "Y": 1, "Z": 2}
+    rows = [_row("good", prod, "A", rep, 5.0 + 2.0 * truth[prod]) for prod in products for rep in (1, 2)]
+    rows += [_row("flat", prod, "A", rep, 5.0) for prod in products for rep in (1, 2)]  # no variance
+    rows += [_row("one", "W", "A", rep, 6.0) for rep in (1, 2)]  # only one product
+    result = panel_consistency(pd.DataFrame(rows), n_permutations=50, random_state=0)
+    table = result.table
+
+    assert np.isnan(table.loc["flat", "p_discrimination"])
+    assert np.isnan(table.loc["flat", "agreement"])
+    assert not bool(table.loc["flat", "consistent"])
+    assert not bool(table.loc["one", "consistent"])
+    assert {"flat", "one"}.issubset(result.inconsistent)
+
+
+def test_consistency_rater_without_panel_overlap_has_no_agreement():
+    """A rater discriminating over products the rest of the panel never rated gets no agreement."""
+    truth = {"W": -1, "X": 0, "Y": 1, "Z": 2}
+    rows = [_row("solo", p, "A", rep, 5.0 + 2.0 * truth[p]) for p in "WXYZ" for rep in (1, 2)]
+    # Two other raters cover only product W, so ``solo`` shares a single cell with the panel.
+    rows += [_row("b1", "W", "A", rep, 5.0) for rep in (1, 2)]
+    rows += [_row("b2", "W", "A", rep, 5.2) for rep in (1, 2)]
+    row = panel_consistency(pd.DataFrame(rows), n_permutations=30, random_state=0).table.loc["solo"]
+    assert not np.isnan(row["p_discrimination"])  # solo still gets a discrimination p-value
+    assert np.isnan(row["agreement"])  # but no agreement can be formed against the panel
+    assert np.isnan(row["p_agreement"])
+
+
+def test_consistency_is_deterministic_and_validates():
+    """Same seed gives identical output; a non-positive permutation count is rejected."""
+    panel = _panel(anomalous="P8")
+    first = panel_consistency(panel, n_permutations=99, random_state=1)
+    second = panel_consistency(panel, n_permutations=99, random_state=1)
+    assert first.table.equals(second.table)
+    assert first.inconsistent == second.inconsistent
+    with pytest.raises(ValueError, match="n_permutations must be at least 1"):
+        panel_consistency(panel, n_permutations=0)
+
+
 def test_dropping_panelist_changes_means():
     validated = validate_descriptive(_panel(), _obs(), mode="observational")
     kept = analyze_descriptive(validated, drop_panelists=None, discriminator=False)