From c6b21511393119834a367a91db2477a8da4601e5 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 23 Jul 2026 20:33:59 +0000 Subject: [PATCH 1/3] Add designed-mode panel comparison: factorial ANOVA + Tukey/Dunnett post-hoc New module process_improve.sensory.designed complements the observational relate with the designed-mode question: which product treatments differ, and by how much, on each attribute, for a randomized complete block design (the same panelists score every treatment, panelist as the block). Public API (generic in the factor column names): - factorial_anova: per-attribute Type III factorial ANOVA (score ~ C(factor_1) * C(factor_2) * ... + C(block)), handling unbalanced grids and testing factor-by-factor interactions. - tukey_hsd: all-pairwise Tukey HSD off the blocked-model error mean square and the studentized-range distribution, with a compact-letter display. - dunnett_vs_control: Dunnett two-sided test of each treatment vs one control. - compare_products: orchestrates the ANOVA + post-hoc tests and returns a ComparisonResult; a within= argument runs simple-effects post-hoc within each level of another factor. Built on the existing statsmodels/patsy/scipy dependencies (no new deps). Tests in tests/test_sensory_designed.py drive a synthetic RCBD with planted effects and cross-check the Tukey maths against scipy.stats.tukey_hsd. Version bump 1.58.0 -> 1.59.0 with matching CITATION.cff and CHANGELOG entry. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_018Ra2KAdJeGtnUzi68yDFaR --- CHANGELOG.md | 28 +- CITATION.cff | 2 +- pyproject.toml | 2 +- src/process_improve/sensory/__init__.py | 12 + src/process_improve/sensory/designed.py | 532 ++++++++++++++++++++++++ tests/test_sensory_designed.py | 200 +++++++++ 6 files changed, 773 insertions(+), 3 deletions(-) create mode 100644 src/process_improve/sensory/designed.py create mode 100644 tests/test_sensory_designed.py diff --git a/CHANGELOG.md b/CHANGELOG.md index a4a9826..adc6174 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,31 @@ those changes. ## [Unreleased] +## [1.59.0] - 2026-07-23 + +### Added + +- New module `process_improve.sensory.designed` for designed-mode comparison of + descriptive panel data, the complement to the existing observational relate. It + answers "which product treatments differ, and by how much, on each attribute" for a + randomized complete block design (the same panelists score every treatment, with + panelist as the block). Public entry points, all generic in the factor column names: + - `factorial_anova(panel, *, factors, block, interactions)`: per-attribute Type III + factorial ANOVA (`score ~ C(factor_1) * C(factor_2) * ... + C(block)`), so an + unbalanced grid is handled correctly and the interaction terms test whether one + factor's effect depends on another (e.g. does aging change some formulations more + than others). + - `tukey_hsd(panel, *, factor, block, alpha)`: all-pairwise Tukey HSD using the + blocked-model error mean square and the studentized-range distribution (the block + variance is removed), with a compact-letter display grouping treatments that are + not separable. + - `dunnett_vs_control(panel, *, factor, control, alpha)`: Dunnett's two-sided test of + every treatment against one named control. + - `compare_products(...)` orchestrates the ANOVA plus post-hoc tests and returns a + `ComparisonResult`; a `within` argument runs the post-hoc tests as simple effects + within each level of another factor (the right follow-up once an interaction is + significant). + ## [1.58.0] - 2026-07-23 ### Added @@ -2599,7 +2624,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 diff --git a/CITATION.cff b/CITATION.cff index a544588..fca9fa5 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.58.0 +version: 1.59.0 date-released: "2026-07-23" keywords: - chemometrics diff --git a/pyproject.toml b/pyproject.toml index 3996999..7fc5fc3 100644 --- a/pyproject.toml +++ b/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" diff --git a/src/process_improve/sensory/__init__.py b/src/process_improve/sensory/__init__.py index 63443bb..f2df141 100644 --- a/src/process_improve/sensory/__init__.py +++ b/src/process_improve/sensory/__init__.py @@ -25,6 +25,13 @@ relate_designed, relate_observational, ) +from process_improve.sensory.designed import ( + ComparisonResult, + compare_products, + dunnett_vs_control, + factorial_anova, + tukey_hsd, +) 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 @@ -39,6 +46,7 @@ "DESCRIPTIVE_LONG_COLUMNS", "SENSORY_RECIPES", "AnalysisResult", + "ComparisonResult", "MAMResult", "PanelScorecard", "ValidationResult", @@ -46,7 +54,10 @@ "align_scores", "analyze_descriptive", "apply_correction", + "compare_products", "discriminate_observational", + "dunnett_vs_control", + "factorial_anova", "mixed_assessor_model", "panel_scorecard", "permutation_column_null", @@ -54,5 +65,6 @@ "relate_designed", "relate_observational", "reshape_to_long", + "tukey_hsd", "validate_descriptive", ] diff --git a/src/process_improve/sensory/designed.py b/src/process_improve/sensory/designed.py new file mode 100644 index 0000000..8d950ae --- /dev/null +++ b/src/process_improve/sensory/designed.py @@ -0,0 +1,532 @@ +"""(c) Kevin Dunn, 2010-2026. MIT License. + +Designed-mode comparison of descriptive panel data: factorial ANOVA of the +product (formulation) effect, with all-pairwise and against-a-control post-hoc +multiple comparisons. + +The observational half of this subpackage relates attributes to *measured* +covariates of products whose formulation is unknown. This module covers the +complementary **designed** question: the products are controlled treatments +(formulations, aging conditions, ...), the same panelists score every treatment +(a randomized complete block design, with panelist as the block), and we want to +know which treatments differ, and by how much, on each attribute. + +Per attribute the workhorse is a fixed-effects ANOVA + + score ~ C(factor_1) * C(factor_2) * ... + C(block) + +fitted by ordinary least squares with **Type III** sums of squares, so an +unbalanced grid (missing cells, a panelist who skipped a sample) is handled +correctly and the interaction terms test whether one factor's effect depends on +another (e.g. does aging change some formulations more than others). When the +omnibus factor effect is real, two post-hoc procedures answer the follow-up +question: + +* :func:`tukey_hsd` - all-pairwise Tukey HSD, using the blocked-model error + mean square and the studentized-range distribution, so the block (panelist) + variance is removed from the yardstick. Answers "which treatments differ from + which". A compact-letter display groups treatments that are not separable. +* :func:`dunnett_vs_control` - Dunnett's two-sided test of every treatment + against a single named control, controlling the family-wise error for the + "many treatments, one control" comparison only (tighter than Tukey). + +:func:`compare_products` runs the whole sequence and returns a +:class:`ComparisonResult`. Everything is generic in the factor column names, so +the same code serves a one-factor formulation screen or a +formulation-by-aging-condition stability study. + +References +---------- +Tukey, J. W. (1949). Comparing individual means in the analysis of variance. +*Biometrics*, 5(2), 99-114. + +Dunnett, C. W. (1955). A multiple comparison procedure for comparing several +treatments with a control. *JASA*, 50(272), 1096-1121. + +Piepho, H.-P. (2004). An algorithm for a letter-based representation of +all-pairwise comparisons. *Journal of Computational and Graphical Statistics*, +13(2), 456-466. + +Naes, T., Brockhoff, P. B. & Tomic, O. (2010). *Statistics for Sensory and +Consumer Science*. Wiley. +""" + +from __future__ import annotations + +import re +from dataclasses import dataclass +from itertools import combinations +from typing import Any + +import numpy as np +import pandas as pd +from scipy.stats import dunnett as _scipy_dunnett +from scipy.stats import studentized_range +from statsmodels.formula.api import ols +from statsmodels.stats.anova import anova_lm + +from process_improve.univariate.metrics import confidence_interval + +#: Regular expression that recovers the original column name from a patsy +#: ``C(Q('name'))`` term, so ANOVA tables read in the caller's own vocabulary. +_Q_TERM = re.compile(r"Q\('([^']+)'\)") + + +@dataclass +class ComparisonResult: + """Outcome of :func:`compare_products`. + + Attributes + ---------- + anova : pandas.DataFrame + Type III ANOVA table, one row per (attribute, source) with ``df``, + ``sum_sq``, ``mean_sq``, ``F`` and ``p_value``. The ``Residual`` row + carries the error mean square used by the post-hoc tests. + tukey : pandas.DataFrame + All-pairwise Tukey HSD contrasts (see :func:`tukey_hsd`), prefixed with + the stratum column when the comparison was stratified. + dunnett : pandas.DataFrame + Dunnett-vs-control contrasts (see :func:`dunnett_vs_control`); empty when + no ``control`` was given. + letters : pandas.DataFrame + Compact-letter display: one row per (stratum, attribute, level) with a + ``letters`` string. Treatments that share a letter are not separable at + the chosen ``alpha``. + means : pandas.DataFrame + Per (stratum, attribute, level) mean, confidence interval and ``n``. + config : dict + The resolved call arguments, for provenance. + """ + + anova: pd.DataFrame + tukey: pd.DataFrame + dunnett: pd.DataFrame + letters: pd.DataFrame + means: pd.DataFrame + config: dict[str, Any] + + +def _clean_term(term: str) -> str: + """Turn a patsy term such as ``C(Q('a')):C(Q('b'))`` back into ``a:b``.""" + names = _Q_TERM.findall(term) + return ":".join(names) if names else term + + +def _formula(factors: list[str], block: str | None, *, interactions: bool) -> str: + """Build the OLS formula for one attribute.""" + joiner = " * " if (interactions and len(factors) > 1) else " + " + terms = joiner.join(f"C(Q('{f}'))" for f in factors) + if block is not None: + terms = f"{terms} + C(Q('{block}'))" + return f"score ~ {terms}" + + +def factorial_anova( + panel: pd.DataFrame, + *, + factors: list[str], + block: str | None = "panelist_id", + interactions: bool = True, +) -> pd.DataFrame: + """Type III factorial ANOVA of the panel scores, one model per attribute. + + Parameters + ---------- + panel : pandas.DataFrame + Descriptive-long panel data (columns ``attribute`` and ``score`` plus the + factor and block columns named below). + factors : list of str + Column names of the fixed factors of interest (e.g. + ``["formulation", "condition"]``). With more than one factor and + ``interactions=True`` their full crossed model is fitted. + block : str or None + Column treated as a blocking factor (default ``"panelist_id"``). Pass + ``None`` for no block. + interactions : bool + Include the factor-by-factor interaction terms (default ``True``). + + Returns + ------- + pandas.DataFrame + One row per (attribute, source) with ``df``, ``sum_sq``, ``mean_sq``, + ``F`` and ``p_value``. The ``Residual`` row gives the error term. An + attribute whose model cannot be fitted (too few observations, a singular + design) yields a single ``source="(model failed)"`` row rather than + aborting the sweep. + + Examples + -------- + >>> factorial_anova(panel, factors=["formulation", "condition"]).head() + """ + missing = [c for c in [*factors, *([block] if block else []), "attribute", "score"] if c not in panel.columns] + if missing: + raise KeyError(f"panel is missing required column(s): {missing}") + + formula = _formula(factors, block, interactions=interactions) + rows: list[dict[str, Any]] = [] + for attribute in sorted(panel["attribute"].unique()): + sub = panel[panel["attribute"] == attribute].dropna(subset=["score"]) + try: + model = ols(formula, data=sub).fit() + table = anova_lm(model, typ=3) + except Exception as exc: # noqa: BLE001 - degenerate design should not abort the sweep + rows.append( + { + "attribute": str(attribute), + "source": "(model failed)", + "df": float("nan"), + "sum_sq": float("nan"), + "mean_sq": float("nan"), + "F": float("nan"), + "p_value": float("nan"), + "note": str(exc).splitlines()[0][:120], + } + ) + continue + for term, record in table.iterrows(): + if term == "Intercept": + continue + df_term = float(record["df"]) + ss_term = float(record["sum_sq"]) + rows.append( + { + "attribute": str(attribute), + "source": _clean_term(str(term)), + "df": df_term, + "sum_sq": ss_term, + "mean_sq": ss_term / df_term if df_term > 0 else float("nan"), + "F": float(record["F"]) if not pd.isna(record["F"]) else float("nan"), + "p_value": float(record["PR(>F)"]) if not pd.isna(record["PR(>F)"]) else float("nan"), + "note": "", + } + ) + return pd.DataFrame(rows) + + +def _group_means(sub: pd.DataFrame, factor: str) -> pd.DataFrame: + """Per-level mean and non-missing count for one attribute subset.""" + grouped = sub.dropna(subset=["score"]).groupby(factor, observed=True)["score"] + return pd.DataFrame({"mean": grouped.mean(), "n": grouped.count()}) + + +def _blocked_error(sub: pd.DataFrame, factor: str, block: str | None) -> tuple[float, float]: + """Return (mean-square error, error degrees of freedom) of ``score ~ factor + block``.""" + formula = _formula([factor], block, interactions=False) + model = ols(formula, data=sub.dropna(subset=["score"])).fit() + return float(model.mse_resid), float(model.df_resid) + + +def tukey_hsd( + panel: pd.DataFrame, + *, + factor: str, + block: str | None = "panelist_id", + alpha: float = 0.05, +) -> pd.DataFrame: + """All-pairwise Tukey HSD of ``factor`` levels, one comparison per attribute. + + The critical difference uses the error mean square of the blocked model + ``score ~ C(factor) + C(block)`` and the studentized-range distribution, so + the block (panelist) variance is removed. Unequal group sizes use the + Tukey-Kramer standard error. + + Parameters + ---------- + panel : pandas.DataFrame + Descriptive-long panel data. + factor : str + Column whose levels are compared all-pairwise. + block : str or None + Blocking column (default ``"panelist_id"``); ``None`` for no block. + alpha : float + Family-wise significance level (default ``0.05``). + + Returns + ------- + pandas.DataFrame + One row per (attribute, pair) with ``group1``, ``group2``, ``meandiff`` + (group1 minus group2), ``se``, ``q_stat``, ``p_value``, the + ``(ci_low, ci_high)`` simultaneous interval and ``reject``. + """ + rows: list[dict[str, Any]] = [] + for attribute in sorted(panel["attribute"].unique()): + sub = panel[panel["attribute"] == attribute].dropna(subset=["score"]) + means = _group_means(sub, factor) + levels = list(means.index) + n_levels = len(levels) + if n_levels < 2: + continue + try: + mse, df_err = _blocked_error(sub, factor, block) + except Exception: # noqa: BLE001, S112 - degenerate attribute is skipped + continue + if not (df_err > 0 and mse > 0): + continue + q_crit = float(studentized_range.ppf(1.0 - alpha, n_levels, df_err)) + for a, b in combinations(levels, 2): + mean_a, mean_b = float(means.loc[a, "mean"]), float(means.loc[b, "mean"]) + n_a, n_b = float(means.loc[a, "n"]), float(means.loc[b, "n"]) + se = float(np.sqrt(mse / 2.0 * (1.0 / n_a + 1.0 / n_b))) + diff = mean_a - mean_b + q_stat = abs(diff) / se if se > 0 else float("nan") + p_value = float(studentized_range.sf(q_stat, n_levels, df_err)) if se > 0 else float("nan") + half = q_crit * se + rows.append( + { + "attribute": str(attribute), + "group1": str(a), + "group2": str(b), + "meandiff": diff, + "se": se, + "q_stat": q_stat, + "p_value": p_value, + "ci_low": diff - half, + "ci_high": diff + half, + "reject": bool(p_value < alpha), + } + ) + return pd.DataFrame(rows) + + +def dunnett_vs_control( + panel: pd.DataFrame, + *, + factor: str, + control: str, + alpha: float = 0.05, +) -> pd.DataFrame: + """Dunnett's two-sided test of every ``factor`` level against ``control``. + + Uses :func:`scipy.stats.dunnett`, which pools the within-level variance and + controls the family-wise error for the many-treatments-vs-one-control family + (tighter than all-pairwise Tukey when the control is the only reference of + interest). Unlike :func:`tukey_hsd` it does not remove a block term. + + Parameters + ---------- + panel : pandas.DataFrame + Descriptive-long panel data. + factor : str + Column whose levels are compared to the control. + control : str + The level of ``factor`` used as the reference. + alpha : float + Significance level (default ``0.05``). + + Returns + ------- + pandas.DataFrame + One row per (attribute, level) for every non-control level, with + ``meandiff`` (level minus control), ``statistic``, ``p_value`` and + ``reject``. + """ + rows: list[dict[str, Any]] = [] + for attribute in sorted(panel["attribute"].unique()): + sub = panel[panel["attribute"] == attribute].dropna(subset=["score"]) + if control not in set(sub[factor].unique()): + continue + control_scores = sub.loc[sub[factor] == control, "score"].to_numpy() + others = [lvl for lvl in sorted(sub[factor].unique()) if lvl != control] + samples = [sub.loc[sub[factor] == lvl, "score"].to_numpy() for lvl in others] + usable = [(lvl, arr) for lvl, arr in zip(others, samples, strict=True) if arr.size >= 1] + if len(control_scores) < 2 or not usable or any(arr.size < 2 for _, arr in usable): + continue + control_mean = float(np.mean(control_scores)) + try: + result = _scipy_dunnett(*[arr for _, arr in usable], control=control_scores) + except Exception: # noqa: BLE001, S112 - degenerate attribute is skipped + continue + stats = np.atleast_1d(result.statistic) + pvals = np.atleast_1d(result.pvalue) + for (lvl, arr), stat, pval in zip(usable, stats, pvals, strict=True): + rows.append( + { + "attribute": str(attribute), + "level": str(lvl), + "control": str(control), + "meandiff": float(np.mean(arr)) - control_mean, + "statistic": float(stat), + "p_value": float(pval), + "reject": bool(pval < alpha), + } + ) + return pd.DataFrame(rows) + + +def _compact_letter_display(levels_by_mean: list[str], sig_pairs: set[frozenset[str]]) -> dict[str, str]: + """Map each level to a letter string (Piepho 2004 insert-and-absorb). + + ``levels_by_mean`` is ordered by descending mean so the first letter tends to + mark the highest group. ``sig_pairs`` holds the frozenset ``{a, b}`` for every + pair judged significantly different. Levels sharing a letter are not separable. + """ + columns: list[set[str]] = [set(levels_by_mean)] + for pair in sig_pairs: + a, b = tuple(pair) + rebuilt: list[set[str]] = [] + for col in columns: + if a in col and b in col: + rebuilt.append(col - {a}) + rebuilt.append(col - {b}) + else: + rebuilt.append(col) + # Absorb: drop empties and any column that is a subset of another. + rebuilt = [c for c in rebuilt if c] + keep: list[set[str]] = [] + for col in rebuilt: + if any(col < other for other in rebuilt if col is not other) or any(col == other for other in keep): + continue + keep.append(col) + columns = keep + + # Order columns by the position of their highest-mean member, then letter them. + order = {lvl: i for i, lvl in enumerate(levels_by_mean)} + columns.sort(key=lambda col: min(order[lvl] for lvl in col)) + letters: dict[str, list[str]] = {lvl: [] for lvl in levels_by_mean} + for idx, col in enumerate(columns): + tag = chr(ord("a") + idx) if idx < 26 else f"({idx + 1})" + for lvl in col: + letters[lvl].append(tag) + return {lvl: "".join(tags) for lvl, tags in letters.items()} + + +def _letters_table(tukey: pd.DataFrame, means: pd.DataFrame, factor: str, alpha: float) -> pd.DataFrame: + """Build the compact-letter-display table from a Tukey frame and level means.""" + rows: list[dict[str, Any]] = [] + for attribute, grp in means.groupby("attribute", observed=True): + ordered = list(grp.sort_values("mean", ascending=False)[factor].astype(str)) + pairs = tukey[tukey["attribute"] == attribute] + sig = {frozenset((str(r.group1), str(r.group2))) for r in pairs.itertuples() if bool(r.reject)} + mapping = _compact_letter_display(ordered, sig) + rows.extend({"attribute": str(attribute), factor: lvl, "letters": mapping.get(lvl, "")} for lvl in ordered) + return pd.DataFrame(rows) + + +def _means_table(panel: pd.DataFrame, factor: str, conf_level: float) -> pd.DataFrame: + """Per (attribute, level) mean with a confidence interval and count.""" + rows: list[dict[str, Any]] = [] + for (attr, lvl), grp in panel.groupby(["attribute", factor], observed=True): + scores = grp[["score"]].dropna() + center = float(scores["score"].mean()) if not scores.empty else float("nan") + if scores.shape[0] >= 2: + lo, hi = confidence_interval(scores, "score", conflevel=conf_level, style="regular") + else: + lo = hi = float("nan") + rows.append( + { + "attribute": str(attr), + factor: str(lvl), + "mean": center, + "ci_low": float(lo), + "ci_high": float(hi), + "n": int(scores.shape[0]), + } + ) + return pd.DataFrame(rows) + + +def compare_products( # noqa: PLR0913 + panel: pd.DataFrame, + *, + factors: list[str], + block: str | None = "panelist_id", + primary: str | None = None, + within: str | None = None, + control: str | None = None, + interactions: bool = True, + alpha: float = 0.05, + conf_level: float = 0.95, +) -> ComparisonResult: + """Compare product treatments per attribute: ANOVA plus post-hoc contrasts. + + Fits the factorial ANOVA over ``factors`` (with a blocking factor), then runs + all-pairwise Tukey HSD and, when ``control`` is given, Dunnett vs the control, + on the ``primary`` factor. When ``within`` is set the post-hoc tests are run + as *simple effects* separately within each level of that factor (e.g. compare + formulations within each aging condition), which is the right follow-up once + the ``primary x within`` interaction is significant. + + Parameters + ---------- + panel : pandas.DataFrame + Descriptive-long panel data. + factors : list of str + Fixed factors for the ANOVA (e.g. ``["formulation", "condition"]``). + block : str or None + Blocking column (default ``"panelist_id"``). + primary : str or None + Factor whose levels the post-hoc tests compare. Defaults to + ``factors[0]``. + within : str or None + If given, run the post-hoc tests separately within each level of this + factor (simple effects). If ``None``, they pool over the other factors. + control : str or None + Level of ``primary`` used as the Dunnett reference. If ``None``, the + Dunnett table is empty. + interactions : bool + Include interaction terms in the ANOVA (default ``True``). + alpha : float + Post-hoc significance level (default ``0.05``). + conf_level : float + Confidence level for the reported means (default ``0.95``). + + Returns + ------- + ComparisonResult + The ANOVA table, Tukey and Dunnett contrasts, compact-letter display and + per-level means; see the class docstring. + + Examples + -------- + >>> res = compare_products( + ... panel, factors=["formulation", "condition"], within="condition", control="Control" + ... ) + >>> res.letters.query("condition == 'REF'").head() + """ + primary = primary or factors[0] + stratum_col = within or "stratum" + anova = factorial_anova(panel, factors=factors, block=block, interactions=interactions) + + strata = [None] if within is None else sorted(panel[within].dropna().unique()) + tukey_frames: list[pd.DataFrame] = [] + dunnett_frames: list[pd.DataFrame] = [] + letter_frames: list[pd.DataFrame] = [] + means_frames: list[pd.DataFrame] = [] + for stratum in strata: + sub = panel if stratum is None else panel[panel[within] == stratum] + label = "all" if stratum is None else str(stratum) + + means = _means_table(sub, primary, conf_level) + tukey = tukey_hsd(sub, factor=primary, block=block, alpha=alpha) + letters = _letters_table(tukey, means, primary, alpha) + for frame in (means, tukey, letters): + frame.insert(0, stratum_col, label) + means_frames.append(means) + tukey_frames.append(tukey) + letter_frames.append(letters) + + if control is not None: + dunnett = dunnett_vs_control(sub, factor=primary, control=control, alpha=alpha) + dunnett.insert(0, stratum_col, label) + dunnett_frames.append(dunnett) + + empty_dunnett = pd.DataFrame( + columns=[stratum_col, "attribute", "level", "control", "meandiff", "statistic", "p_value", "reject"] + ) + return ComparisonResult( + anova=anova, + tukey=pd.concat(tukey_frames, ignore_index=True) if tukey_frames else pd.DataFrame(), + dunnett=pd.concat(dunnett_frames, ignore_index=True) if dunnett_frames else empty_dunnett, + letters=pd.concat(letter_frames, ignore_index=True) if letter_frames else pd.DataFrame(), + means=pd.concat(means_frames, ignore_index=True) if means_frames else pd.DataFrame(), + config={ + "factors": list(factors), + "block": block, + "primary": primary, + "within": within, + "control": control, + "interactions": interactions, + "alpha": alpha, + "conf_level": conf_level, + }, + ) diff --git a/tests/test_sensory_designed.py b/tests/test_sensory_designed.py new file mode 100644 index 0000000..5d8b997 --- /dev/null +++ b/tests/test_sensory_designed.py @@ -0,0 +1,200 @@ +"""Tests for the designed-mode comparison module (:mod:`process_improve.sensory.designed`). + +A single synthetic randomized-complete-block scenario drives most tests: seven +panelists score five formulations (one named ``Control``) at three aging +conditions (``REF``, ``RT``, ``HB``) on two attributes. On attribute ``A`` the +formulation ``T4`` is planted well above the four others (which are identical), +aging lowers every score, and a formulation-by-condition interaction is planted +so that ``T2`` collapses only under ``HB``. Attribute ``B`` carries no +formulation signal (a negative control). Panelist offsets make the block matter. +Because the truth is known, the ANOVA F-tests, Tukey groupings, +Dunnett-vs-control flags and letter display can each be checked against it. +""" + +from __future__ import annotations + +import numpy as np +import pandas as pd +import pytest + +from process_improve.sensory import ( + ComparisonResult, + compare_products, + dunnett_vs_control, + factorial_anova, + tukey_hsd, +) +from process_improve.sensory.designed import _compact_letter_display + +PANELISTS = [f"P{i}" for i in range(7)] +FORMULATIONS = ["Control", "T1", "T2", "T3", "T4"] +CONDITIONS = ["REF", "RT", "HB"] +_FORM_EFFECT = {"Control": 0.0, "T1": 0.0, "T2": 0.0, "T3": 0.0, "T4": 3.0} +_COND_EFFECT = {"REF": 0.0, "RT": -0.5, "HB": -1.0} +_NOISE_SD = 0.35 + + +def make_rcbd_panel(seed: int = 0) -> pd.DataFrame: + """Return a descriptive-long panel with the planted effects described above.""" + rng = np.random.default_rng(seed) + panelist_offset = {pid: rng.normal(0.0, 0.5) for pid in PANELISTS} + rows: list[dict[str, object]] = [] + for pid in PANELISTS: + for form in FORMULATIONS: + for cond in CONDITIONS: + interaction = -2.0 if (form == "T2" and cond == "HB") else 0.0 + mean_a = 5.0 + _FORM_EFFECT[form] + _COND_EFFECT[cond] + interaction + panelist_offset[pid] + mean_b = 5.0 + _COND_EFFECT[cond] + panelist_offset[pid] # no formulation signal + for attribute, center in (("A", mean_a), ("B", mean_b)): + rows.append( + { + "panelist_id": pid, + "product": f"{form}-{cond}", + "formulation": form, + "condition": cond, + "attribute": attribute, + "replicate": 1, + "score": center + rng.normal(0.0, _NOISE_SD), + } + ) + return pd.DataFrame(rows) + + +@pytest.fixture(scope="module") +def panel() -> pd.DataFrame: + return make_rcbd_panel(seed=0) + + +def test_factorial_anova_detects_planted_effects(panel: pd.DataFrame) -> None: + table = factorial_anova(panel, factors=["formulation", "condition"]) + attr_a = table[table["attribute"] == "A"].set_index("source") + # Every planted term is strongly significant on attribute A. + assert attr_a.loc["formulation", "p_value"] < 1e-6 + assert attr_a.loc["condition", "p_value"] < 1e-3 + assert attr_a.loc["formulation:condition", "p_value"] < 1e-3 + assert attr_a.loc["panelist_id", "p_value"] < 1e-3 + # Degrees of freedom: 5 formulations, 3 conditions, 7 panelists. + assert attr_a.loc["formulation", "df"] == pytest.approx(4.0) + assert attr_a.loc["condition", "df"] == pytest.approx(2.0) + assert attr_a.loc["formulation:condition", "df"] == pytest.approx(8.0) + + +def test_factorial_anova_null_attribute_has_no_formulation_effect(panel: pd.DataFrame) -> None: + table = factorial_anova(panel, factors=["formulation", "condition"]) + attr_b = table[table["attribute"] == "B"].set_index("source") + # Attribute B carries no formulation signal, so its effect should not be significant. + assert attr_b.loc["formulation", "p_value"] > 0.05 + assert attr_b.loc["formulation:condition", "p_value"] > 0.05 + # Aging is still present on B. + assert attr_b.loc["condition", "p_value"] < 1e-3 + + +def test_factorial_anova_missing_column_raises(panel: pd.DataFrame) -> None: + with pytest.raises(KeyError, match="missing required column"): + factorial_anova(panel, factors=["not_a_column"]) + + +def test_tukey_only_high_formulation_separates_at_ref(panel: pd.DataFrame) -> None: + ref = panel[panel["condition"] == "REF"] + tukey = tukey_hsd(ref, factor="formulation") + attr_a = tukey[tukey["attribute"] == "A"] + + def rejected(a: str, b: str) -> bool: + row = attr_a[ + ((attr_a["group1"] == a) & (attr_a["group2"] == b)) | ((attr_a["group1"] == b) & (attr_a["group2"] == a)) + ] + return bool(row["reject"].iloc[0]) + + # T4 separates from all four others; none of the low four separate from each other. + for other in ["Control", "T1", "T2", "T3"]: + assert rejected("T4", other) + for a, b in [("Control", "T1"), ("Control", "T2"), ("T1", "T3"), ("T2", "T3")]: + assert not rejected(a, b) + + +def test_tukey_matches_scipy_on_balanced_oneway() -> None: + """Without a block, the studentized-range maths must match scipy.stats.tukey_hsd.""" + scipy_stats = pytest.importorskip("scipy.stats") + rng = np.random.default_rng(3) + groups = {g: rng.normal(loc, 1.0, 8) for g, loc in zip("WXYZ", [0.0, 0.5, 2.0, 2.2], strict=True)} + rows = [ + {"panelist_id": f"P{i}", "attribute": "A", "score": val, "grp": g, "replicate": 1} + for g, arr in groups.items() + for i, val in enumerate(arr) + ] + long = pd.DataFrame(rows) + ours = tukey_hsd(long, factor="grp", block=None).set_index(["group1", "group2"]) + reference = scipy_stats.tukey_hsd(*groups.values()) + labels = list(groups.keys()) + for i, a in enumerate(labels): + for j, b in enumerate(labels): + if i < j: + key = (a, b) if (a, b) in ours.index else (b, a) + assert ours.loc[key, "p_value"] == pytest.approx(reference.pvalue[i, j], abs=1e-6) + + +def test_dunnett_flags_only_the_true_treatment(panel: pd.DataFrame) -> None: + ref = panel[panel["condition"] == "REF"] + dunnett = dunnett_vs_control(ref, factor="formulation", control="Control") + attr_a = dunnett[dunnett["attribute"] == "A"].set_index("level") + assert bool(attr_a.loc["T4", "reject"]) + assert attr_a.loc["T4", "meandiff"] > 2.0 + for other in ["T1", "T2", "T3"]: + assert not bool(attr_a.loc[other, "reject"]) + + +def test_compare_products_simple_effects_capture_interaction(panel: pd.DataFrame) -> None: + result = compare_products(panel, factors=["formulation", "condition"], within="condition", control="Control") + assert isinstance(result, ComparisonResult) + assert set(result.letters["condition"]) == set(CONDITIONS) + + def letter(condition: str, formulation: str) -> str: + row = result.letters.query("attribute == 'A' and condition == @condition and formulation == @formulation") + return str(row["letters"].iloc[0]) + + # Under HB, T2 collapses to its own group, distinct from the Control cluster. + assert letter("HB", "T4") != letter("HB", "Control") + assert letter("HB", "T2") != letter("HB", "Control") + # Under REF, T2 sits with the Control cluster (no collapse there). + assert letter("REF", "T2") == letter("REF", "Control") + + +def test_compare_products_without_control_has_empty_dunnett(panel: pd.DataFrame) -> None: + result = compare_products(panel, factors=["formulation", "condition"], within="condition") + assert result.dunnett.empty + assert not result.tukey.empty + assert not result.means.empty + + +def test_unbalanced_panel_is_handled(panel: pd.DataFrame) -> None: + # Drop a whole panelist-by-formulation cell and a scatter of individual scores. + unbalanced = panel[~((panel["panelist_id"] == "P0") & (panel["formulation"] == "T3"))].copy() + rng = np.random.default_rng(11) + drop_idx = rng.choice(unbalanced.index, size=15, replace=False) + unbalanced = unbalanced.drop(index=drop_idx) + + table = factorial_anova(unbalanced, factors=["formulation", "condition"]) + assert not table.empty + # Type III still recovers the dominant formulation effect on A despite the imbalance. + attr_a = table[table["attribute"] == "A"].set_index("source") + assert attr_a.loc["formulation", "p_value"] < 1e-4 + + result = compare_products(unbalanced, factors=["formulation", "condition"], within="condition", control="Control") + assert not result.tukey.empty + assert not result.letters.empty + + +def test_compact_letter_display_handles_non_interval_pattern() -> None: + # A can equal B and C, but B and C differ: A must carry both letters. + levels = ["A", "B", "C"] # ordered by mean + sig = {frozenset(("B", "C"))} + mapping = _compact_letter_display(levels, sig) + assert mapping["B"] != mapping["C"] + # A shares a letter with B and a letter with C. + assert set(mapping["A"]) & set(mapping["B"]) + assert set(mapping["A"]) & set(mapping["C"]) + + +def test_compact_letter_display_all_same_when_nothing_significant() -> None: + mapping = _compact_letter_display(["A", "B", "C"], set()) + assert mapping["A"] == mapping["B"] == mapping["C"] == "a" From a9647d397bb03f9e7a35c637d210e2048f729e3f Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 24 Jul 2026 06:47:40 +0000 Subject: [PATCH 2/3] docs: document designed-mode comparison with a worked Tukey/Dunnett example Add a "Comparing designed treatments" section to the sensory panel user guide covering compare_products / factorial_anova / tukey_hsd / dunnett_vs_control, the randomized-complete-block model, and a worked numerical example that shows Tukey and Dunnett reaching different (both correct) verdicts on the same treatment-vs-control comparison. Cross-linked from the designed-versus- observational section. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_018Ra2KAdJeGtnUzi68yDFaR --- docs/user_guide/sensory_panel.rst | 144 ++++++++++++++++++++++++++++++ 1 file changed, 144 insertions(+) diff --git a/docs/user_guide/sensory_panel.rst b/docs/user_guide/sensory_panel.rst index 0a5b4ca..abe5920 100644 --- a/docs/user_guide/sensory_panel.rst +++ b/docs/user_guide/sensory_panel.rst @@ -68,6 +68,15 @@ Keep the interpretation in mind: an observational analysis supports only "this descriptor is associated with this attribute", whereas the future designed analysis will support "increasing this factor raises this attribute". +A related but separate question does not need the covariate table at all: when +the products *are* the controlled treatments (for example five formulations, or +a formulation crossed with an aging condition) and you simply want to know which +treatments differ on each attribute, use the direct factorial-comparison API in +:func:`~process_improve.sensory.compare_products` (see +:ref:`comparing-designed-treatments` below). That path is implemented now; the +``mode="designed"`` relate switch above is the still-planned integration of the +same idea into the covariate-driven pipeline. + Step 0: get the data into long form ----------------------------------- @@ -262,6 +271,141 @@ Rescaling does not remove genuine disagreement, so a panelist who truly ranks the products differently is better handled by dropping (``drop_panelists`` or ``correction="drop"``); align and drop can be combined. +.. _comparing-designed-treatments: + +Comparing designed treatments: ANOVA, Tukey HSD and Dunnett +----------------------------------------------------------- + +When the products are controlled treatments rather than measured market samples +(five formulations, a formulation crossed with an aging condition, a process +setting), and the same panelists score every treatment, the design is a +**randomized complete block** with panelist as the block. The question is which +treatments differ, and by how much, on each attribute. +:func:`~process_improve.sensory.compare_products` answers it with a per-attribute +factorial ANOVA followed by the matching post-hoc multiple comparisons; the +pieces are also usable on their own +(:func:`~process_improve.sensory.factorial_anova`, +:func:`~process_improve.sensory.tukey_hsd`, +:func:`~process_improve.sensory.dunnett_vs_control`). + +Per attribute the model is a Type III factorial ANOVA + +.. math:: + + \text{score} \sim C(\text{factor}_1) \times C(\text{factor}_2) \times \dots + C(\text{block}) + +Type III sums of squares keep the effects correct when the grid is unbalanced (a +panelist who skipped a sample, a missing cell), and the interaction terms test +whether one factor's effect depends on another, for example whether aging shifts +some formulations more than others. + +.. code-block:: python + + from process_improve.sensory import compare_products + + # panel: descriptive_long, with extra factor columns "formulation" and "condition". + result = compare_products( + panel, + factors=["formulation", "condition"], + block="panelist_id", + within="condition", # run the post-hoc tests within each condition (simple effects) + control="Control", # the reference level for Dunnett + ) + + result.anova # Type III ANOVA table, one row per (attribute, source) + result.tukey # all-pairwise Tukey HSD contrasts + result.dunnett # each formulation vs the Control + result.letters # compact-letter display: shared letter => not separable + result.means # per-level mean with a confidence interval + +Set ``within="condition"`` to run the post-hoc tests as *simple effects* +separately within each aging condition, which is the right follow-up once the +``formulation`` x ``condition`` interaction is significant; leave it ``None`` to +compare the primary factor pooled over the others. + +Two post-hoc procedures, two questions +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Both procedures control the family-wise error rate (the chance of *any* false +claim across the whole family of comparisons), but over different families: + +- **Tukey HSD** compares **every pair** of treatments. Its critical difference + uses the blocked-model error mean square :math:`\text{MSE}` and the + studentized-range distribution :math:`q`, so the panelist-block variance is + removed from the yardstick: + + .. math:: + + |\bar{y}_i - \bar{y}_j| \;>\; q_{\alpha,\,g,\,\nu}\,\sqrt{\text{MSE}/n} + + where :math:`g` is the number of treatments and :math:`\nu` the error degrees + of freedom. The compact-letter display in ``result.letters`` summarises it: + treatments that share a letter are not separable. + +- **Dunnett** compares **each treatment to one control**. Because those + comparisons all share the control mean they are correlated, so the critical + value comes from the multivariate-*t* distribution rather than the studentized + range. Guarding fewer comparisons makes Dunnett more powerful than Tukey for + the treatment-vs-control question, but it cannot compare two treatments to each + other. + +A rule of thumb: use **Dunnett** when there is a genuine reference you are +benchmarking against (a Control formulation, or REF for the aging axis), and +**Tukey** when the treatments are peers and you want the full pairwise map. +Running both is common; report them as answering different questions. + +Worked numbers +~~~~~~~~~~~~~~~ + +Take one attribute scored by five formulations, seven panelists each, with an +error mean square of :math:`\text{MSE}=1.0` on :math:`\nu=30` degrees of freedom. +The level means are: + +========= ======= ==== ==== ==== ==== +Statistic Control T1 T2 T3 T4 +========= ======= ==== ==== ==== ==== +mean 3.0 3.3 3.5 4.5 6.2 +========= ======= ==== ==== ==== ==== + +The standard error of a group mean is :math:`\sqrt{\text{MSE}/n}=\sqrt{1/7}=0.378` +and of a difference is :math:`\sqrt{2\,\text{MSE}/n}=0.535`. + +*Tukey.* With :math:`q_{0.05,\,5,\,30}=4.102`, the critical difference is +:math:`4.102 \times 0.378 = 1.55`. Only the gaps involving T4 exceed it, so the +letter display is T4 = ``a`` and Control, T1, T2, T3 = ``b``. + +*Dunnett.* The two-sided critical value for four treatments versus one control at +:math:`\nu=30` is about :math:`2.18`, giving a critical difference of +:math:`2.18 \times 0.535 = 1.16`. + +The two cutoffs disagree on exactly one comparison: + +=============== ========== ===================== ======================== +Comparison Difference Tukey (cutoff 1.55) Dunnett (cutoff 1.16) +=============== ========== ===================== ======================== +T1 - Control 0.30 not significant not significant +T2 - Control 0.50 not significant not significant +T3 - Control 1.50 not significant significant (p=0.03) +T4 - Control 3.20 significant significant +=============== ========== ===================== ======================== + +T3 sits 1.50 above the Control. Tukey, which must protect all ten pairwise +comparisons, raises its bar to 1.55 and calls it non-significant; Dunnett, which +protects only the four treatment-vs-control comparisons, has a lower bar (1.16) +and flags it. Same data, both correct: the procedures answer different questions +with correctly calibrated thresholds. + +.. rubric:: Post-hoc references + +- Tukey, J. W. (1949). Comparing individual means in the analysis of variance. + *Biometrics*, 5(2), 99-114. +- Dunnett, C. W. (1955). A multiple comparison procedure for comparing several + treatments with a control. *Journal of the American Statistical Association*, + 50(272), 1096-1121. +- Piepho, H.-P. (2004). An algorithm for a letter-based representation of + all-pairwise comparisons. *Journal of Computational and Graphical Statistics*, + 13(2), 456-466. + Worked example -------------- From 941eae32820447d5740bd3f489e3edbbd086f3e7 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 24 Jul 2026 06:55:41 +0000 Subject: [PATCH 3/3] Pin ruff below 0.16 to restore the lint gate ruff 0.16 promoted several rules into the select=["ALL"] set (CPY001, PLR0917, and others) that fire across the existing code base, turning the lint job red on every branch. Cap the dev dependency at ruff>=0.11.0,<0.16 (resolves to 0.15.22) so `ruff check .` passes again; the new rules can be triaged separately before lifting the cap. Version bump 1.59.0 -> 1.59.1 with matching CITATION.cff and CHANGELOG entry. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_018Ra2KAdJeGtnUzi68yDFaR --- CHANGELOG.md | 12 +++++++++++- CITATION.cff | 2 +- pyproject.toml | 12 +++++++++--- 3 files changed, 21 insertions(+), 5 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index adc6174..f748450 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,15 @@ those changes. ## [Unreleased] +## [1.59.1] - 2026-07-23 + +### Changed + +- Pin the dev-tooling `ruff` below `0.16` (`ruff>=0.11.0,<0.16`). ruff 0.16 + promoted several rules into the `select=["ALL"]` set (`CPY001`, `PLR0917`, + ...) that fire across the existing code base and broke the lint gate. The cap + keeps CI green until those rules are triaged. + ## [1.59.0] - 2026-07-23 ### Added @@ -2624,7 +2633,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.59.0...HEAD +[Unreleased]: https://github.com/kgdunn/process-improve/compare/v1.59.1...HEAD +[1.59.1]: https://github.com/kgdunn/process-improve/compare/v1.59.0...v1.59.1 [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 diff --git a/CITATION.cff b/CITATION.cff index fca9fa5..5a9f21d 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.59.0 +version: 1.59.1 date-released: "2026-07-23" keywords: - chemometrics diff --git a/pyproject.toml b/pyproject.toml index 7fc5fc3..64a32d0 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "process-improve" -version = "1.59.0" +version = "1.59.1" description = 'Designed Experiments; Latent Variables (PCA, PLS, multivariate methods with missing data); Process Monitoring; Batch data analysis.' readme = "README.md" license = "MIT" @@ -116,7 +116,10 @@ dev = [ "pytest-cov>=7.0.0", "pytest-xdist>=3.8.0", "pydata-sphinx-theme>=0.16.1", - "ruff>=0.11.0", + # Upper-bounded below 0.16: ruff 0.16 promoted several rules into the + # select=["ALL"] set (CPY001, PLR0917, ...) that fire across the existing + # code base. Revisit and triage those rules before lifting the cap. + "ruff>=0.11.0,<0.16", "sphinx>=8.1.3", "nbsphinx>=0.9.5", "ipykernel>=6.29", @@ -158,7 +161,10 @@ dev = [ "pytest-cov>=7.0.0", "pytest-xdist>=3.8.0", "pydata-sphinx-theme>=0.16.1", - "ruff>=0.11.0", + # Upper-bounded below 0.16: ruff 0.16 promoted several rules into the + # select=["ALL"] set (CPY001, PLR0917, ...) that fire across the existing + # code base. Revisit and triage those rules before lifting the cap. + "ruff>=0.11.0,<0.16", "sphinx>=8.1.3", "nbsphinx>=0.9.5", "ipykernel>=6.29",