"""(c) Kevin Dunn, 2010-2026. MIT License.
Relate descriptive panel attributes to the product.
:func:`analyze_descriptive` runs the proof-of-concept pipeline on a validated
dataset: score and (optionally) correct the panel, then relate each sensory
attribute to the product. The relate step dispatches on the validation mode:
* **observational** (supported) - the product has measured descriptors but
unknown formulation, so the attribute block is related to the descriptors
with PLS (:class:`process_improve.multivariate.PLS` plus VIP) and
per-descriptor correlations, reported as association rather than causation.
* **designed** (stub, not implemented yet) - the product is a controlled
experimental run; the plan is to regress each attribute on the design factors
via :func:`process_improve.experiments.analyze_experiment` for factor effects.
See :func:`relate_designed`; it raises ``NotImplementedError`` for now.
The observational relate corrects across the family of tests with
Benjamini-Hochberg FDR and returns supporting product means with confidence
intervals and a PCA sensory map. Both the marginal associations and the
per-attribute predictive-descriptor search are additionally gated on a
leave-one-out jackknife, so an association or predictive coefficient that rests on a single
high-leverage observation (a predictor that is non-zero on only one product,
common in sparse, wide descriptor blocks) is demoted rather than reported. The
jackknife adds no threshold of its own: it reuses the same ``alpha`` and the
number of observations, so a genuine multi-observation driver is unaffected.
"""
from __future__ import annotations
import itertools
import warnings
from dataclasses import dataclass, field
from typing import Any
import numpy as np
import pandas as pd
from scipy.stats import pearsonr
from scipy.stats import t as t_dist
from process_improve.multivariate._common import NotEnoughVarianceError
from process_improve.multivariate.methods import PCA, PLS, MCUVScaler, selectivity_ratio, vip
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.validation import ValidationResult
from process_improve.univariate.metrics import benjamini_hochberg, confidence_interval
[docs]
@dataclass
class AnalysisResult:
"""Outcome of :func:`analyze_descriptive`.
Attributes
----------
mode : str
``"designed"`` or ``"observational"``.
panel : PanelScorecard
The per-panelist scorecard and flags.
dropped : list of str
Panelists removed before the relate step.
mam : MAMResult
Mixed Assessor Model: per-panelist scaling coefficients and the MAM vs
classical product-effect F-tests.
correction : str
The panel correction applied before relating: ``"none"``, ``"align"``,
or ``"drop"``.
relate : dict
Mode-specific relate results; see :func:`analyze_descriptive`.
product_means : pandas.DataFrame
Per product-by-attribute mean with a confidence interval.
pca : dict
Product sensory map: ``scores``, ``loadings``, ``explained_variance``.
config : dict
The options the analysis ran with.
"""
mode: str
panel: PanelScorecard
dropped: list[str]
mam: MAMResult
correction: str
relate: dict[str, Any]
product_means: pd.DataFrame
pca: dict[str, Any]
config: dict[str, Any] = field(default_factory=dict)
[docs]
def aggregate_to_product(panel: pd.DataFrame) -> pd.DataFrame:
"""Return a product-by-attribute table of mean scores.
Parameters
----------
panel : pandas.DataFrame
Validated ``descriptive_long`` panel data.
Returns
-------
pandas.DataFrame
Index ``product``, one column per attribute, values the mean score
over panelists and replicates.
"""
wide = panel.pivot_table(index="product", columns="attribute", values="score", aggfunc="mean", observed=True)
wide.index = wide.index.astype(str)
wide.columns.name = None
return wide
#: Minimum usable observations before the leave-one-out jackknife of an
#: association is defined; below this an association cannot be certified as
#: influence-robust (mirrors the ``len < 4`` guard in ``panel._mad_bands``).
_MIN_OBS_FOR_JACKKNIFE = 4
#: Correlations are clipped off +/-1 before the Fisher-z transform so ``arctanh``
#: stays finite.
_R_CLIP = 1.0 - 1e-12
def _resolve_find_predictive(find_predictive: bool, discriminator: bool | None) -> bool:
"""Accept the deprecated ``discriminator=`` spelling of ``find_predictive=``.
.. deprecated:: 1.77.0
``discriminator`` will be removed in 2.0.0.
"""
if discriminator is None:
return find_predictive
if not find_predictive:
msg = (
"Pass either 'find_predictive' or the deprecated 'discriminator', not both. "
"They set the same thing: whether to run the per-attribute predictive-descriptor search."
)
raise ValueError(msg)
warnings.warn(
"The 'discriminator' argument is deprecated since 1.77.0 and will be removed in 2.0.0; "
"use 'find_predictive' instead.",
category=DeprecationWarning,
stacklevel=3,
)
return discriminator
def _attach_fdr(records: list[dict[str, Any]], alpha: float) -> list[dict[str, Any]]:
"""Attach Benjamini-Hochberg q-values (and reject flags) to ``records``."""
pvals = [r["p_value"] for r in records]
if not pvals:
return records
bh = benjamini_hochberg(np.asarray(pvals), alpha=alpha)
for rec, q, rej in zip(records, bh.p_adjusted, bh.reject, strict=True):
rec["q_value"] = float(q)
# Harden the marginal significance: an association counts only when it also
# survives the leave-one-out jackknife, so a single high-leverage
# observation cannot manufacture a "significant" correlation.
rec["significant"] = bool(rej) and bool(rec.get("influence_robust", True))
return records
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. ``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)
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]
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 = float(pseudo.std(ddof=1) / np.sqrt(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]:
"""Fit PLS, stepping the component count down on a near-collinear (singular) block.
Near-duplicate descriptor columns (a proxy that is almost an exact function
of a driver) can make the high-order PLS deflation singular. Retry with one
fewer component until the fit succeeds, returning ``(model, components_used)``
or ``(None, 0)`` if even a single component fails.
"""
for k in range(max(1, n_components), 0, -1):
try:
with warnings.catch_warnings():
warnings.simplefilter("ignore", RuntimeWarning)
return PLS(n_components=k).fit(x, y), k
except np.linalg.LinAlgError: # noqa: PERF203 - retry on singular block; loop is a few iterations
continue
return None, 0
def _collinear_clusters(x_block: pd.DataFrame, threshold: float) -> dict[str, int]:
"""Group descriptors into clusters of mutually high absolute correlation.
Single-linkage connected components on the descriptor ``|corr|`` matrix:
two descriptors join the same cluster when their absolute Pearson
correlation is at least ``threshold``. Returns ``{descriptor: cluster_id}``
with cluster ids assigned in column order (a singleton descriptor gets its
own id). Collinear proxies therefore share an id, which is how
:func:`find_predictive_descriptors` reports that they cannot be told apart.
"""
cols = list(x_block.columns)
n = len(cols)
corr = np.nan_to_num(x_block.corr().abs().to_numpy(), nan=0.0)
parent = list(range(n))
def find(i: int) -> int:
while parent[i] != i:
parent[i] = parent[parent[i]]
i = parent[i]
return i
for i in range(n):
for j in range(i + 1, n):
if corr[i, j] >= threshold:
ri, rj = find(i), find(j)
if ri != rj:
parent[max(ri, rj)] = min(ri, rj)
roots: dict[int, int] = {}
cluster_of: dict[str, int] = {}
for i, name in enumerate(cols):
root = find(i)
if root not in roots:
roots[root] = len(roots)
cluster_of[str(name)] = roots[root]
return cluster_of
[docs]
def find_predictive_descriptors( # noqa: PLR0913, PLR0915
agg: pd.DataFrame,
covariates: pd.DataFrame,
*,
n_components: int = 2,
alpha: float = 0.05,
n_permutations: int = 199,
random_state: int = 0,
cluster_threshold: float = 0.95,
max_components_cv: int = 4,
) -> dict[str, Any]:
"""Find which descriptors carry predictive signal, per attribute.
The marginal associations (:func:`relate_observational`) flag every
descriptor that correlates with an attribute in-sample, genuine drivers and
proxies alike. This step adds out-of-sample evidence:
1. a per-attribute cross-validated Q-squared gate (is the attribute
predictable from the descriptor block at all),
2. a selectivity ratio per descriptor on the target-projected predictive
direction, with a permutation p-value corrected for multiplicity by the
max-statistic (Westfall-Young) permutation, so a descriptor that merely
correlates by chance but does not enter the predictive direction is
demoted, and
3. a collinear-cluster id per descriptor.
What it cannot do is rank descriptors *within* a collinear cluster: two
descriptors that carry the same information predict equally well out of
sample, so they share a cluster id and both stay significant. Separating
them needs an external dataset or a designed experiment.
.. important::
**The multiplicity correction is within an attribute, not across them.**
The max-statistic null is rebuilt for each attribute over its own
descriptors, so ``p_value_fwer`` controls the family-wise error rate for
*that attribute's* descriptor family only. Nothing corrects across
attributes: on a panel of ``A`` attributes at ``alpha``, roughly
``alpha * A`` attributes are expected to produce a spurious family by
chance alone. Read a single flagged descriptor on a many-attribute panel
with that in mind.
Parameters
----------
agg : pandas.DataFrame
Product-by-attribute mean table (index ``product``).
covariates : pandas.DataFrame
One row per product with the measured descriptors (plus a ``product``
column, which is dropped here).
n_components : int
Latent components for the in-sample selectivity-ratio fit.
alpha : float
Target false-discovery rate for the permutation family.
n_permutations : int
Number of label permutations for the selectivity-ratio null.
random_state : int
Seed for the permutations and the cross-validation folds.
cluster_threshold : float
Absolute-correlation threshold for the collinear clustering.
max_components_cv : int
Cap on the component count the Q-squared gate may select.
Returns
-------
dict
``per_attribute`` (the Q-squared gate per attribute), ``descriptors``
(the per attribute-descriptor selectivity ratio, raw permutation
``p_value``, family-wise-error-adjusted ``p_value_fwer``,
``jackknife_significant`` flag from the leave-one-out beta confidence
interval, ``is_predictive`` flag and ``cluster_id``), ``clusters`` (the
descriptor-to-cluster map), and the settings used. A descriptor is
``is_predictive`` only when it also survives the jackknife, so a
coefficient carried by a single product is demoted.
See Also
--------
permutation_column_null :
The block-level counterpart. It fits one multi-response PLS over the
whole attribute block and returns one record per descriptor, answering
"which descriptors matter for the panel as a whole" rather than "which
matter for this attribute". Reach for it to screen a descriptor block
before committing to per-attribute work; reach for this function when
the answer has to name the attribute.
"""
x_all = covariates.loc[agg.index]
descriptors = [c for c in x_all.columns if c != "product" and pd.api.types.is_numeric_dtype(x_all[c])]
x_block = x_all[descriptors].astype(float)
clusters = _collinear_clusters(x_block, cluster_threshold)
rng = np.random.default_rng(random_state)
per_attribute: list[dict[str, Any]] = []
records: list[dict[str, Any]] = []
for attr in agg.columns:
y = agg[attr].astype(float)
mask = y.notna()
x_attr = x_block.loc[mask]
y_attr = y.loc[mask]
n_rows = x_attr.shape[0]
cap = max(1, min(max_components_cv, x_attr.shape[1], n_rows - 2))
a = max(1, min(n_components, cap))
# Fit one PLS for this attribute, reused for the Q-squared gate and the
# selectivity ratio. ``_fit_pls_safe`` steps the component count down if
# the near-collinear descriptor block makes the fit singular.
x_scaled = MCUVScaler().fit_transform(x_attr)
y_scaled = MCUVScaler().fit_transform(y_attr.to_frame())
pls, a = _fit_pls_safe(x_scaled, y_scaled, a)
# 1. Leave-one-out cross-validated Q-squared gate: is the attribute
# predictable from the descriptor block out of sample? The same LOO
# refit also yields a jackknife confidence interval per descriptor
# coefficient (Martens' uncertainty test); it is reused below so a
# descriptor whose predictive weight rests on a single high-leverage
# product is demoted even when it survives the permutation null.
#
# `q2_cv > 0.0` is deliberately a low, uncalibrated bar. It is a
# cheap pre-screen, not the test: it decides whether the 199-refit
# permutation loop below is worth running, and it is ANDed with the
# calibrated max-statistic p-value and the jackknife flag when
# `is_predictive` is set. Passing it therefore cannot make a
# descriptor a finding; only failing it can rule one out early.
# Tightening it would remove findings that already cleared a
# family-wise-error-controlled null, at the cost of a permutation
# null per attribute on top of the one already here.
predictable = False
q2_cv = float("nan")
rmsep_cv = float("nan")
jack_significant: dict[str, bool] = {}
if pls is not None and n_rows >= 5 and y_attr.std() > 0:
with warnings.catch_warnings():
warnings.simplefilter("ignore", RuntimeWarning)
cv = pls.cross_validate(x_scaled, y_scaled, cv="loo", conf_level=1.0 - alpha, show_progress=False)
q2_cv = float(cv.q_squared.iloc[0])
rmsep_cv = float(cv.rmse_cv.iloc[0])
predictable = q2_cv > 0.0
jack_significant = {str(name): bool(flag) for name, flag in cv.significant.iloc[:, 0].items()}
per_attribute.append(
{
"attribute": str(attr),
"n_components_cv": a if pls is not None else 0,
"q2_cv": q2_cv,
"rmsep_cv": rmsep_cv,
"predictable": bool(predictable),
}
)
# 2. Selectivity ratio on the predictive direction, with a permutation
# null. Multiplicity across descriptors is controlled by the
# max-statistic (Westfall-Young) permutation: for each label
# permutation, the *largest* selectivity ratio over all descriptors
# forms the null. An observed SR above that null is significant after
# correction, so even a single genuine driver is detectable without
# the resolution loss of a per-test Benjamini-Hochberg floor. The
# permutation loop is skipped for attributes the Q-squared gate found
# unpredictable: none of their descriptors can be flagged anyway.
k = len(descriptors)
if pls is None: # degenerate block: nothing to relate for this attribute
sr_values = np.zeros(k)
p_raw = np.ones(k)
p_maxt = np.ones(k)
else:
sr_values = np.asarray(selectivity_ratio(pls, x_scaled).reindex(descriptors), dtype=float)
if predictable:
y_values = y_scaled.to_numpy().ravel()
ge_each = np.zeros(k) # per-descriptor null exceedances
ge_max = np.zeros(k) # exceedances of the family-wide max null
done = 0
for _ in range(n_permutations):
permuted = pd.DataFrame(
y_values[rng.permutation(n_rows)], index=x_scaled.index, columns=y_scaled.columns
)
pls_p, _ = _fit_pls_safe(x_scaled, permuted, a)
if pls_p is None:
continue # a degenerate permutation contributes nothing to the null
sr_p = selectivity_ratio(pls_p, x_scaled).reindex(descriptors).to_numpy()
ge_each += sr_p >= sr_values
ge_max += np.nanmax(sr_p) >= sr_values
done += 1
denom = done + 1.0
p_raw = (ge_each + 1.0) / denom
p_maxt = (ge_max + 1.0) / denom
else:
p_raw = np.ones(k)
p_maxt = np.ones(k)
for i, desc in enumerate(descriptors):
desc_robust = jack_significant.get(str(desc), False)
records.append(
{
"attribute": str(attr),
"descriptor": str(desc),
"selectivity_ratio": float(sr_values[i]),
"p_value": float(p_raw[i]),
"p_value_fwer": float(p_maxt[i]),
# Deprecated since 1.77.0, removed in 2.0.0: q_value is the
# old name for p_value_fwer. It was always a family-wise
# error rate, never an FDR q-value, which is why it moved.
"q_value": float(p_maxt[i]),
"jackknife_significant": bool(desc_robust),
"is_predictive": bool(p_maxt[i] <= alpha and predictable and desc_robust),
# Deprecated since 1.77.0, removed in 2.0.0.
"discriminator_significant": bool(p_maxt[i] <= alpha and predictable and desc_robust),
"cluster_id": clusters[str(desc)],
}
)
return {
"per_attribute": per_attribute,
"descriptors": records,
"clusters": clusters,
"alpha": alpha,
"n_permutations": n_permutations,
"cluster_threshold": cluster_threshold,
}
[docs]
def relate_designed(
agg: pd.DataFrame,
covariates: pd.DataFrame,
*,
model: str = "main_effects",
alpha: float = 0.05,
) -> dict[str, Any]:
"""Relate attributes to controlled design factors (not implemented yet).
Stub for a later release. The plan is to regress each attribute on the
design factors via :func:`process_improve.experiments.analyze_experiment`
(or ``analyze_omars`` for DSD/OMARS designs) and report factor effects with
Benjamini-Hochberg correction. For now use ``mode="observational"``.
Raises
------
NotImplementedError
Always, until the designed-mode relate step is built.
"""
del agg, covariates, model, alpha
raise NotImplementedError(
"Designed (DoE/OMARS) relate is not implemented yet; use mode='observational'. Planned for a later release."
)
[docs]
def relate_observational( # noqa: PLR0913
agg: pd.DataFrame,
covariates: pd.DataFrame,
*,
n_components: int = 2,
alpha: float = 0.05,
find_predictive: bool = True,
n_permutations: int = 199,
random_state: int = 0,
influence_deletions: int = 1,
discriminator: bool | None = None,
) -> 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
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.
"""
find_predictive = _resolve_find_predictive(find_predictive, discriminator)
x_block = covariates.loc[agg.index].astype(float)
y_block = agg.astype(float)
max_comp = max(1, min(n_components, x_block.shape[1], x_block.shape[0] - 1))
pls = PLS(n_components=max_comp).fit(x_block, y_block)
vips = vip(pls)
drivers: list[dict[str, Any]] = [{"descriptor": str(name), "vip": float(value)} for name, value in vips.items()]
drivers.sort(key=lambda r: float(r["vip"]), reverse=True)
# Per (attribute, descriptor) association, BH-corrected across the family.
assoc: list[dict[str, Any]] = []
for attr in y_block.columns:
for desc in x_block.columns:
pair = pd.concat([y_block[attr], x_block[desc]], axis=1).dropna()
if pair.shape[0] >= 3 and pair.iloc[:, 0].std() > 0 and pair.iloc[:, 1].std() > 0:
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, max_deletions=influence_deletions)
assoc.append(
{
"attribute": str(attr),
"descriptor": str(desc),
"r": float(r),
"p_value": float(p),
"jackknife_se": float(jack_se),
"influence_robust": bool(robust),
"n_supporting": int(n_support),
}
)
assoc = _attach_fdr(assoc, alpha)
result: dict[str, Any] = {
"mode": "observational",
"n_components": max_comp,
"alpha": alpha,
"vip": drivers,
"associations": assoc,
}
if find_predictive:
result["predictive_descriptors"] = find_predictive_descriptors(
agg,
covariates.loc[agg.index],
n_components=max_comp,
alpha=alpha,
n_permutations=n_permutations,
random_state=random_state,
)
# Deprecated since 1.77.0, removed in 2.0.0: the same object under its
# old key, so an existing caller reading result["discriminator"] works.
result["discriminator"] = result["predictive_descriptors"]
return result
def _knockoff_block(x_block: pd.DataFrame, k: int, rng: np.random.Generator) -> pd.DataFrame:
"""Return ``k`` null columns, each a row-permutation of a randomly chosen real column.
Sampling the source columns with replacement lets ``k`` exceed the number of real
columns (needed for the ``min_knockoffs`` floor on a narrow block); a permuted copy
keeps the source column's own marginal - its spread, sparsity, and support - so the
null is matched to the data rather than to a simulated distribution.
"""
n_rows, p = x_block.shape
source = rng.integers(0, p, size=k)
columns = {f"__null_{j}": x_block.iloc[:, src].to_numpy()[rng.permutation(n_rows)] for j, src in enumerate(source)}
return pd.DataFrame(columns, index=x_block.index)
[docs]
def permutation_column_null( # noqa: PLR0913
agg: pd.DataFrame,
covariates: pd.DataFrame,
*,
ignore: list[str] | None = None,
n_components: int = 2,
fraction: float = 0.15,
min_knockoffs: int = 7,
max_knockoffs: int | None = None,
n_iter: int = 200,
quantile: float = 0.95,
random_state: int = 0,
) -> dict[str, Any]:
"""Empirical VIP / cross-validated-beta null for the descriptor block.
Adds ``k`` permuted "knockoff" columns - each a row-shuffled copy of a real
descriptor (:func:`_knockoff_block`) - to the descriptor block, fits the PLS relate,
and reads the VIP and cross-validated beta of every column. Repeated over ``n_iter``
permutations, the knockoff columns form an empirical null band: a real descriptor is
only credible if its VIP / beta clears a high quantile of the null the permuted
columns achieve. Because even a descriptor with no real relationship earns a
non-trivial VIP in a ``p >> n`` fit, this calibrates the magnitude against the data's
own permuted columns rather than a parametric cutoff.
This is decoupled from the influence gate: it does not itself remove any descriptors.
Pass the descriptors the gate demoted (single or twin-support spikes) as ``ignore``;
they are dropped from the fit entirely - not merely skipped when building knockoffs -
so they no longer distort the scores or VIP of the survivors.
Parameters
----------
agg : pandas.DataFrame
Product-by-attribute mean table (index ``product``).
covariates : pandas.DataFrame
One row per product with the measured descriptors (index ``product``; a
``product`` column, if present, is ignored).
ignore : list of str, optional
Descriptor names to drop from the fit before building the null (default: none).
A name absent from the descriptor block raises ``ValueError`` so a typo fails
loudly instead of silently doing nothing.
n_components : int
Latent components for each PLS fit.
fraction : float
Fraction of surviving descriptors used as the knockoff count.
min_knockoffs : int
Floor on the knockoff count, so a narrow block still gets a usable null.
max_knockoffs : int, optional
Optional cap on the knockoff count (default: uncapped).
n_iter : int
Number of permutations (refits); more gives a smoother null threshold.
quantile : float
Null quantile used as the significance threshold (e.g. 0.95).
random_state : int
Seed for the permutations.
Returns
-------
dict
``descriptors`` (per surviving descriptor: ``vip`` / ``cv_beta`` and the
``*_null_threshold`` and ``*_exceeds_null`` fields), plus the settings used and
the counts (``n_descriptors``, ``n_knockoffs``, ``n_iter``, ``ignored``).
See Also
--------
find_predictive_descriptors :
The per-attribute counterpart, and the one to prefer when the answer
has to name an attribute. It reports one record per (attribute,
descriptor) pair and controls the family-wise error rate within each
attribute by a max-statistic permutation. This function instead fits a
single multi-response PLS over the whole attribute block and returns one
record per descriptor, so it answers "which descriptors matter for the
panel as a whole"; its knockoff quantile band is a calibrated screen
rather than formal error control. Use it to triage a wide descriptor
block before committing to the per-attribute work.
"""
if fraction <= 0.0:
raise ValueError(f"fraction must be positive, got {fraction}.")
if not 0.0 < quantile < 1.0:
raise ValueError(f"quantile must be in (0, 1), got {quantile}.")
if min_knockoffs < 1 or n_iter < 1:
raise ValueError("min_knockoffs and n_iter must both be at least 1.")
ignore = list(ignore or [])
x_all = covariates.loc[agg.index]
descriptors = [c for c in x_all.columns if c != "product" and pd.api.types.is_numeric_dtype(x_all[c])]
unknown = sorted(set(ignore) - set(descriptors))
if unknown:
raise ValueError(f"ignore contains descriptor name(s) not in the covariate block: {unknown}.")
kept = [c for c in descriptors if c not in set(ignore)]
p = len(kept)
config = {
"n_descriptors": p,
"n_ignored": len(ignore),
"ignored": sorted(ignore),
"fraction": fraction,
"quantile": quantile,
"n_iter_requested": n_iter,
}
if p < 2:
return {"ok": False, "reason": f"need at least 2 descriptors after ignoring; got {p}.", **config}
y_block = agg.astype(float)
x_block = x_all[kept].astype(float)
n_rows = x_block.shape[0]
k = max(int(min_knockoffs), round(fraction * p))
if max_knockoffs is not None:
k = min(k, int(max_knockoffs))
max_comp = max(1, min(n_components, p, n_rows - 1))
null_names = [f"__null_{j}" for j in range(k)]
rng = np.random.default_rng(random_state)
real_vip_runs: list[pd.Series] = []
real_beta_runs: list[pd.Series] = []
null_vip: list[float] = []
null_beta: list[float] = []
for _ in range(n_iter):
x_aug = pd.concat([x_block, _knockoff_block(x_block, k, rng)], axis=1)
try:
with warnings.catch_warnings():
warnings.simplefilter("ignore", RuntimeWarning)
pls = PLS(n_components=max_comp).fit(x_aug, y_block)
beta = pls.cross_validate(x_aug, y_block, cv="loo", show_progress=False).beta_mean.abs().max(axis=1)
except (np.linalg.LinAlgError, NotEnoughVarianceError):
continue # a degenerate (singular / no-variance) block contributes nothing to the null
vips = vip(pls)
real_vip_runs.append(vips.reindex(kept))
real_beta_runs.append(beta.reindex(kept))
null_vip.extend(float(vips[name]) for name in null_names)
null_beta.extend(float(beta[name]) for name in null_names)
if not real_vip_runs:
return {"ok": False, "reason": "every permuted fit was singular; no null could be formed.", **config}
real_vip = pd.concat(real_vip_runs, axis=1).mean(axis=1)
real_beta = pd.concat(real_beta_runs, axis=1).mean(axis=1)
vip_threshold = float(np.nanquantile(null_vip, quantile))
beta_threshold = float(np.nanquantile(null_beta, quantile))
records: list[dict[str, Any]] = [
{
"descriptor": desc,
"vip": float(real_vip[desc]),
"vip_null_threshold": vip_threshold,
"vip_exceeds_null": bool(real_vip[desc] > vip_threshold),
"cv_beta": float(real_beta[desc]),
"cv_beta_null_threshold": beta_threshold,
"cv_beta_exceeds_null": bool(real_beta[desc] > beta_threshold),
}
for desc in kept
]
records.sort(key=lambda r: r["vip"], reverse=True)
return {"ok": True, "descriptors": records, "n_knockoffs": k, "n_iter": len(real_vip_runs), **config}
[docs]
def product_means(panel: pd.DataFrame, conf_level: float = 0.95) -> pd.DataFrame:
"""Return per product-by-attribute mean with a confidence interval."""
rows: list[dict[str, Any]] = []
for (prod, attr), grp in panel.groupby(["product", "attribute"], observed=True):
scores = grp[["score"]].dropna()
center = float(scores["score"].mean())
if scores.shape[0] >= 2:
lo, hi = confidence_interval(scores, "score", conflevel=conf_level, style="regular")
else:
lo = hi = float("nan")
rows.append(
{"product": str(prod), "attribute": str(attr), "mean": center, "ci_low": float(lo), "ci_high": float(hi)}
)
return pd.DataFrame(rows)
def _pca_map(agg: pd.DataFrame, n_components: int = 2) -> dict[str, Any]:
"""Fit a PCA on the product-by-attribute means and return the map."""
filled = agg.dropna(axis=1, how="any")
max_comp = max(1, min(n_components, filled.shape[0] - 1, filled.shape[1]))
pca = PCA(n_components=max_comp).fit(filled)
return {
"scores": pca.scores_,
"loadings": pca.loadings_,
"explained_variance": [float(v) for v in pca.r2_per_component_],
}
[docs]
def analyze_descriptive( # noqa: PLR0913
validated: ValidationResult,
*,
drop_panelists: str | list[str] | None = None,
correction: str = "none",
align_method: str = "both",
model: str = "main_effects",
n_components: int = 2,
conf_level: float = 0.95,
alpha: float = 0.05,
find_predictive: bool = True,
n_permutations: int = 199,
random_state: int = 0,
influence_deletions: int = 1,
discriminator: bool | None = None,
) -> AnalysisResult:
"""Run the descriptive pipeline: panel check, correction, and relate.
Parameters
----------
validated : ValidationResult
A passing result from
:func:`process_improve.sensory.validate_descriptive`.
drop_panelists : {"auto", None} or list of str
``"auto"`` drops every flagged panelist; a list drops exactly those
ids; ``None`` keeps all panelists.
correction : {"none", "align", "drop"}
Panel correction before relating. ``"none"`` (default) leaves scores as
is; ``"align"`` applies the Mixed Assessor Model scale alignment to all
panelists (:func:`process_improve.sensory.mam.align_scores`); ``"drop"``
is a synonym for using ``drop_panelists``. Alignment and dropping
compose: panelists are aligned first, then any dropped.
align_method : {"both", "location", "scale"}
Which MAM lever to apply when ``correction="align"``.
model : str
Design model for the ``designed`` relate step (default
``"main_effects"``).
n_components : int
Components for the PLS relate step and the PCA map.
conf_level : float
Confidence level for the product-mean intervals.
alpha : float
Target false-discovery rate for the relate step.
find_predictive : bool
Whether to run the per-attribute predictive-descriptor search
(:func:`find_predictive_descriptors`) in the observational relate step.
n_permutations : int
Permutations for the selectivity-ratio null.
random_state : int
Seed for the 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
-------
AnalysisResult
See the class docstring.
Raises
------
ValueError
If ``validated`` did not pass validation.
"""
if not validated.ok or validated.normalized_df is None or validated.covariates is None:
raise ValueError(
"analyze_descriptive requires a validated dataset; "
"validate_descriptive reported errors: "
f"{validated.errors}"
)
panel = validated.normalized_df
card = panel_scorecard(panel)
mam = mixed_assessor_model(panel)
# Correction: align all panelists onto a common scale (MAM), then drop.
working = align_scores(panel, method=align_method) if correction == "align" else panel
if drop_panelists == "auto":
dropped = list(card.flagged)
elif isinstance(drop_panelists, list):
dropped = drop_panelists
else:
dropped = []
find_predictive = _resolve_find_predictive(find_predictive, discriminator)
clean = apply_correction(working, dropped)
agg = aggregate_to_product(clean)
if validated.mode == "designed":
relate = relate_designed(agg, validated.covariates, model=model, alpha=alpha)
else:
relate = relate_observational(
agg,
validated.covariates,
n_components=n_components,
alpha=alpha,
find_predictive=find_predictive,
n_permutations=n_permutations,
random_state=random_state,
influence_deletions=influence_deletions,
)
return AnalysisResult(
mode=validated.mode,
panel=card,
dropped=dropped,
mam=mam,
correction=correction,
relate=relate,
product_means=product_means(clean, conf_level=conf_level),
pca=_pca_map(agg, n_components=n_components),
config={
"model": model,
"correction": correction,
"align_method": align_method,
"n_components": n_components,
"conf_level": conf_level,
"alpha": alpha,
"find_predictive": find_predictive,
# Deprecated since 1.77.0, removed in 2.0.0.
"discriminator": find_predictive,
"n_permutations": n_permutations,
"random_state": random_state,
"influence_deletions": influence_deletions,
"content_hash": validated.content_hash,
},
)
# ---------------------------------------------------------------------------
# Deprecated aliases - removal scheduled for 2.0.0
# ---------------------------------------------------------------------------
[docs]
def discriminate_observational( # noqa: PLR0913
agg: pd.DataFrame,
covariates: pd.DataFrame,
*,
n_components: int = 2,
alpha: float = 0.05,
n_permutations: int = 199,
random_state: int = 0,
cluster_threshold: float = 0.95,
max_components_cv: int = 4,
) -> dict[str, Any]:
"""Forward to :func:`find_predictive_descriptors`; emits a :class:`DeprecationWarning`.
.. deprecated:: 1.77.0
Use :func:`find_predictive_descriptors` instead. Will be removed in
2.0.0.
"""
warnings.warn(
"process_improve.sensory.discriminate_observational is deprecated since 1.77.0 and will "
"be removed in 2.0.0; use find_predictive_descriptors instead.",
category=DeprecationWarning,
stacklevel=2,
)
return find_predictive_descriptors(
agg,
covariates,
n_components=n_components,
alpha=alpha,
n_permutations=n_permutations,
random_state=random_state,
cluster_threshold=cluster_threshold,
max_components_cv=max_components_cv,
)