Descriptive panel data: validate, check the panel, relate to the product

A descriptive panel study asks a group of assessors (panelists) to score a set of products on a set of attributes, usually with replicates and sometimes across several sessions. The recurring question is: which product characteristics drive the attribute scores, and can we trust the panel that produced them?

The process_improve.sensory subpackage answers that with a small, generic pipeline. The data is described only as panelist, product, attribute, replicate, and score; nothing about the pipeline is specific to any particular kind of product. The flow is:

0. reshape a raw wide export into the long schema, if the data does not already arrive in it;

1. validate the panel data and a product-covariate table;

2. check the panel, then either correct each panelist’s scale use (the Mixed Assessor Model) or drop anomalous panelists;

3. relate each attribute back to the product, and discriminate the descriptors that carry out-of-sample predictive signal from those that only correlate in this one sample.

The worked example near the end runs all of these on a synthetic dataset.

The data contract

Panel data is supplied in the descriptive_long schema, one row per score:

Column	Meaning
panelist_id	Who gave the score.
session	Which sitting the score came from.
product	Which product was scored.
attribute	Which attribute was scored.
replicate	Which repeat of that product/attribute for that panelist.
score	The numeric rating.

Alongside it you supply a product-covariate table: one row per product, describing what each product is. This table comes in two flavours, and the distinction decides how the relate step works.

Designed versus observational

The covariate table is one of two kinds, and you must say which with the mode argument:

• Observational (mode="observational") - implemented. You did not set the formulation (for example the products are existing market products), but you measured each one, for example by chemical or instrumental analysis. These measured descriptors are correlated covariates, not a designed matrix, so the relate step reports association, not causation: the attribute block is related to the descriptors with PLS, and per-descriptor correlations show which descriptors track which attributes.

• Designed (mode="designed") - planned, not implemented yet. You controlled the formulation, so each product is an experimental run and the covariates are the design factors. Because the factors were deliberately varied, the relate step will speak of effects: each attribute regressed on the factors, so a coefficient estimates how changing a factor changes the attribute. For now mode="designed" raises NotImplementedError.

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 compare_products() (see Comparing designed treatments: ANOVA, Tukey HSD and Dunnett 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

Panel data usually arrives wide (one column per attribute) rather than in the long schema. reshape_to_long() turns a parsed table into descriptive_long from an explicit column mapping, melting the attribute columns when needed. It verifies round-trip invariants (the grand mean, the mean per attribute, the mean per panelist, and the cell count are identical before and after), so a wrong mapping fails loudly instead of corrupting the analysis silently.

from process_improve.sensory import reshape_to_long

# wide_table: rows are assessor x sample x replicate, one column per attribute.
long_df, checks = reshape_to_long(
    wide_table,
    layout="wide_by_attribute",
    mapping={"panelist_id": "Assessor", "product": "Sample", "replicate": "Rep"},
)
assert checks["ok"]   # grand / per-attribute / per-panelist means preserved

Already-long data passes through with layout="long" and an attribute / score mapping. Means-only tables (no panelist column) are refused, since the Mixed Assessor Model needs panelist-level scores.

Step 1: validate

validate_descriptive() coerces the inputs to the schema and checks them: required columns and dtypes, the score range, panel balance (the fraction of the full panelist x product x attribute x replicate grid that is missing), label encoding, and mode-specific covariate checks. It returns a result whose ok flag gates the rest of the pipeline.

import pandas as pd
from process_improve.sensory import validate_descriptive, analyze_descriptive

# panel: a DataFrame in the descriptive_long schema (loaded from your data).
# descriptors: one row per product, with the measured (e.g. instrumental)
# covariates for each product.
descriptors = pd.DataFrame(
    {"product": ["A", "B", "C", "D"], "sodium": [0.2, 0.5, 0.8, 1.0], "fat": [3.1, 3.0, 2.9, 3.2]}
)

validated = validate_descriptive(panel, descriptors, mode="observational", score_min=0, score_max=10)
print(validated.ok, validated.warnings)

Step 2 and 3: check the panel and relate

analyze_descriptive() runs the rest. It first builds a per-panelist scorecard (panel_scorecard()) rating each panelist on discrimination (do they separate the products), agreement (do they rank products like the rest of the panel), scale use, and drift across sessions. A panelist is flagged only when it is both an outlier and genuinely poor on agreement or discrimination. Passing drop_panelists="auto" removes the flagged panelists before relating, so a noisy panelist does not contaminate the product conclusions.

result = analyze_descriptive(validated, drop_panelists="auto")

print(result.panel.flagged)          # panelists flagged as anomalous
print(result.dropped)                # panelists actually removed

# Observational relate: which descriptors are associated with which attributes.
drivers = pd.DataFrame(result.relate["vip"])             # PLS descriptor importance
assoc = pd.DataFrame(result.relate["associations"])      # attribute-descriptor links
print(assoc[assoc["significant"]])

The observational relate output is:

• result.relate["vip"] ranks descriptors by their PLS variable-importance.

• result.relate["associations"] gives the per-(attribute, descriptor) correlation with a raw p-value, a Benjamini-Hochberg q_value across the whole family of tests, and a significant flag.

• result.relate["discriminator"] adds the out-of-sample evidence described in the next section.

The result also carries supporting context: result.product_means (each product-by-attribute mean with a confidence interval) and result.pca (a PCA sensory map of the products over the attributes).

The designed-mode relate (factor effects rather than associations) is planned; mode="designed" raises NotImplementedError for now.

Telling genuine drivers from proxies: the discriminator

The associations table reports marginal correlations: each attribute-descriptor pair on its own. A pair can be significant for three different reasons, which a marginal correlation cannot tell apart: the descriptor genuinely drives the attribute; the descriptor only rides on a genuine driver (a proxy); or the descriptor happens to line up with the attribute in this particular set of products (a coincidence). The discriminator adds out-of-sample evidence to separate these:

• a per-attribute cross-validated \(Q^2\) (leave-one-out): is the attribute predictable from the descriptor block at all, on products the model did not see? result.relate["discriminator"]["per_attribute"] reports q2_cv and a predictable flag. A coincidence does not predict held-out products, so its attribute often fails this gate.

• a selectivity ratio per descriptor, computed on the target-projected predictive direction (the single PLS component aligned with the attribute). The selectivity ratio is the share of a descriptor’s variance that the predictive direction explains; unlike the marginal correlation it is judged given the other descriptors, so a descriptor that adds nothing beyond the real drivers scores low. Significance is assessed with a permutation test that controls for testing many descriptors at once (a max-statistic null), gated on the attribute being predictable. result.relate["discriminator"]["descriptors"] reports selectivity_ratio, a permutation q_value and a discriminator_significant flag.

• a collinear cluster id per descriptor (cluster_id): descriptors whose absolute correlation exceeds a threshold are grouped. Two descriptors in the same cluster carry the same information, so they predict equally well and the discriminator keeps them both. This is the honest limit of an observational analysis: it can report that a group of descriptors is predictive, but it cannot say which one inside a collinear cluster is the cause. Separating them needs an external dataset that breaks the collinearity, a designed experiment, or mechanistic knowledge.

So the discriminator demotes coincidences (they fail the gate or the selectivity-ratio test) and groups proxies with their driver, but it does not, and cannot, rank descriptors within a collinear cluster.

Under the hood, the two predictive-importance steps are the standalone diagnostics target_projection() and selectivity_ratio() (also available as PLS methods). Target projection rotates the fitted PLS so that one component lies along the regression vector \(b\) for the chosen attribute: the weight is \(w_{\text{TP}} = b / \lVert b \rVert\) and the scores are \(t_{\text{TP}} = X\, w_{\text{TP}}\). The selectivity ratio of descriptor \(j\) is the ratio of its explained to its residual sum of squares on that single component, \(\text{SR}_j = \text{SS}_{\text{explained},j} / \text{SS}_{\text{residual},j}\), so it is large only for descriptors aligned with the predictive direction. Two collinear descriptors carry near-identical selectivity ratios, which is why the clustering step is needed alongside it.

References

• Kvalheim, O. M. and Karstang, T. V. (1989). Interpretation of latent-variable regression models. Chemometrics and Intelligent Laboratory Systems, 7(1-2), 39-51. (target projection)

• Rajalahti, T. et al. (2009). Biomarker discovery in mass spectral profiles by means of selectivity ratio plot. Chemometrics and Intelligent Laboratory Systems, 95(1), 35-48. (selectivity ratio)

Correcting the panel: the Mixed Assessor Model

Panelists differ in how they use the scale: some compress it into a narrow range, some expand it, some sit consistently high or low. The Mixed Assessor Model (MAM) separates this scale usage from genuine disagreement about the products. For each attribute it regresses every panelist’s product scores on the panel consensus; the slope is the panelist’s scaling coefficient beta:

• beta near 1: uses the scale like the panel;

• beta < 1: compresses; beta > 1: expands.

mixed_assessor_model() returns these coefficients per panelist and attribute, plus two product-effect F-tests: the MAM one (using the leftover disagreement as the error term) and the classical one (using the raw interaction). Removing the scale-usage differences from the error makes the MAM F-test more powerful.

Instead of dropping a panelist who merely scales differently, you can align the panel: align_scores() rescales every panelist onto a common scale (a location lever removes their offset, a scale lever divides by beta), keeping their data while removing the scale-usage artefact. analyze_descriptive exposes this through correction:

from process_improve.sensory import mixed_assessor_model, align_scores

mam = mixed_assessor_model(validated.normalized_df)
print(mam.scaling.sort_values("beta").head())   # who compresses / expands
print(mam.ftests)                                # MAM vs classical F per attribute

# Align all panelists onto a common scale, then relate to the product.
result = analyze_descriptive(validated, correction="align")
print(result.correction, result.mam.scaling.head())

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: 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. 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 (factorial_anova(), tukey_hsd(), dunnett_vs_control()).

Per attribute the model is a Type III factorial ANOVA

\[\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.

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 \(\text{MSE}\) and the studentized-range distribution \(q\), so the panelist-block variance is removed from the yardstick:

  \[|\bar{y}_i - \bar{y}_j| \;>\; q_{\alpha,\,g,\,\nu}\,\sqrt{\text{MSE}/n}\]

  where \(g\) is the number of treatments and \(\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 \(\text{MSE}=1.0\) on \(\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 \(\sqrt{\text{MSE}/n}=\sqrt{1/7}=0.378\) and of a difference is \(\sqrt{2\,\text{MSE}/n}=0.535\).

Tukey. With \(q_{0.05,\,5,\,30}=4.102\), the critical difference is \(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 \(\nu=30\) is about \(2.18\), giving a critical difference of \(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.

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

This example runs the whole pipeline on a small synthetic panel. Ten assessors (J01-J10) scored eighteen products (Product A-R) on nine attributes (Aroma intensity, Sweetness, Sourness, Bitterness, Firmness, Juiciness, Colour intensity, Aftertaste, Liking), as integers on a 0-10 scale. The scores came out of the scoring software in the wide layout, with an extra site column we do not need and a few missing cells. We also have instrumental measurements per product, in realistic physical units, split on purpose into three kinds:

• genuine mechanistic correlates: brix for sweetness, titratable_acidity for sourness, polyphenols for bitterness, aroma_oav for aroma, viscosity for firmness;

• proxies that ride on a genuine driver: refractive_index and specific_gravity rise with dissolved sugar (brix); conductivity rides on the acid ions; total_dissolved_solids is an aggregate; price tracks Liking through the sample frame;

• unrelated measurement-condition nuisances with no causal link to any attribute: serving_temperature, headspace_volume, sample_mass and lab_humidity.

Three assessors were constructed to misbehave: J07 scores at random, J03 rates everything high, and J09 uses only the middle of the scale.

Step 1, reshape. The export is wide, so melt the attribute columns to the long schema and ignore site. The round-trip check confirms the grand, per-attribute, and per-assessor means are unchanged.

from process_improve.sensory import (
    reshape_to_long, validate_descriptive, panel_scorecard,
    mixed_assessor_model, analyze_descriptive,
)

long_df, checks = reshape_to_long(
    wide_table,
    layout="wide_by_attribute",
    mapping={"panelist_id": "Assessor", "product": "Product", "ignore": ["site"]},
)
assert checks["ok"]

Step 2, validate.

validated = validate_descriptive(long_df, covariates, mode="observational")

Any out-of-range or missing-cell issues are surfaced as warnings.

Step 3, check the panel. The scorecard flags J07 (low agreement with the panel: its random scores do not track the consensus). The Mixed Assessor Model reports each assessor’s scaling coefficient beta: about 1 for most, well below 1 for the compressor J09, and a large offset for the high rater J03. Once the random assessor is removed, the leftover assessor-by-product interaction is mostly scale usage, so the MAM product F-test exceeds the classical one (the disagreement error term shrinks).

card = panel_scorecard(long_df)
print(card.flagged)                       # ['J07']
mam = mixed_assessor_model(long_df)
print(mam.scaling.sort_values("beta").head())

Step 4, correct and relate. Align all assessors onto a common scale (J03 is brought down, J09’s range is stretched), then relate the attributes to the measurements.

result = analyze_descriptive(validated, correction="align")
import pandas as pd
assoc = pd.DataFrame(result.relate["associations"])
print(assoc[assoc["significant"]])

The marginal test flags sixteen pairs. The genuine correlates are there (brix with Sweetness, titratable_acidity with Sourness, price with Liking), but so are the proxies (refractive_index and specific_gravity with Sweetness, which ride on brix) and even a coincidence: lab_humidity correlates with Liking (r about -0.67) purely by chance in these eighteen products. A within-sample correlation on its own cannot tell these three cases apart.

Step 5, discriminate. The discriminator adds the out-of-sample evidence.

disc = result.relate["discriminator"]
gate = pd.DataFrame(disc["per_attribute"])      # cross-validated Q-squared per attribute
drivers = pd.DataFrame(disc["descriptors"])     # selectivity ratio, q-value, cluster id
print(drivers[drivers["discriminator_significant"]])

It narrows the sixteen marginal hits to twelve:

• The cross-validated \(Q^2\) gate keeps Sweetness (q2_cv about 0.95) and Liking (about 0.94) but rejects Juiciness, Colour intensity and Aftertaste, for which no covariate predicts held-out products.

• lab_humidity is demoted: although it correlates with Liking in-sample, it carries no predictive selectivity, so the permutation test does not keep it, while the genuine brix and price for Liking remain significant. This is the coincidence caught.

• brix, refractive_index and specific_gravity all stay significant for Sweetness and share one cluster_id. The discriminator reports them as a single inseparable group: from these data you cannot say which of the three is the cause, because each carries the same information.

This is the trap, stated precisely: a within-sample correlation is not a transferable, causal link. Cross-validation demotes coincidences, and the selectivity ratio plus clustering groups proxies with their driver, but ranking descriptors within a collinear cluster needs evidence this single panel cannot supply: an external dataset that breaks the collinearity, a designed experiment, or mechanistic knowledge.

The full runnable scenario is in the test suite (tests/test_sensory_end_to_end.py).

Using the tools from an agent

The steps are exposed as agent-callable tools (see process_improve.sensory.tools.get_sensory_tool_specs()), taking the panel and covariate tables as lists of row-records and returning JSON:

• sensory_reshape_to_long - reshape a parsed wide/long table into the descriptive_long schema with round-trip checks.

• sensory_validate_descriptive - validate the inputs.

• sensory_panel_check - panel quality from the panel alone (no covariates): the scorecard with flags, the MAM scaling coefficients and F-tests, and, with align=true, the rescaled panel.

• sensory_analyze_descriptive - the full pipeline, with a correction option ("none" / "align" / "drop"), the MAM results, and (unless discriminator is set false) the cross-validated discriminator in its output.

The analyze tool validates first and refuses to run if validation fails, so an agent cannot skip the gate.