From 960480af875cf3a1b18b39e72ae51960aba7ad20 Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 25 Jul 2026 18:37:08 +0000 Subject: [PATCH 1/5] Compute score contributions from the data, not from the score score_contributions() back-projected a score-space difference through the loadings: contributions = (t_end - t_start) @ P. That expression never receives the observation's data. With one component it reduces exactly to -t_1 * p_1, so it returns the loading vector rescaled by a constant, and every observation in a data set produces the same ranking of variables. On the food-texture data all 50 observations gave one ranking; the correct calculation gives 34. The variable blamed for the movement differed for 36 of those 50 observations, and for 47 of the 54 LDPE observations. That is precisely the failure the diagnostic exists to prevent. Miller, Swanson and Heckler (1994), who introduced contribution plots, show a batch whose loadings point at the silver concentration variables while its contributions point at the salt and silver pressure variables; the pressures were the real fault. A loading describes the whole data set, a contribution describes one observation, and a variable with a large loading contributes nothing when that observation sits at its mean. A contribution is the term-by-term breakdown of the sum that forms the score, so the terms must add up to the score being decomposed: c_ik = x_ik * R_ka sum_k c_ik = t_ia with R the score-generating matrix (loadings for PCA, direct_weights_ for PLS, following the generalisation to any latent-variable model in Westerhuis, Gurden and Smilde, 2000). For the multiblock models the super score at component a is formed from the block data deflated through the earlier components, so the decomposition starts from that same deflated data; summing over every block and variable returns the super score. score_contributions() now takes X and returns one row per observation, matching t2_contributions() and spe_contributions() next to it. Those two were already correct: t2_contributions multiplies by X_values and telescopes to T-squared, which is the pattern this function was missing. The group and two-period forms from the same paper, where x_ik is replaced by a mean or by a difference of means, are exposed as group_contributions(). The paper's two presentation scalings are available via scaling=. The old signature cannot be supported alongside the new one, since a score vector does not carry the information needed. Calls in the old form raise a TypeError naming the replacement rather than returning a different number silently. The previous tests asserted only that the output equalled dt @ P.T, which re-derived the implementation and so held whatever it computed. They now assert the defining property, that the contributions sum to the score, and that the ranking varies across observations. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01D58KmvE1L75DDE7pqfJPeL --- src/process_improve/multivariate/__init__.py | 4 + src/process_improve/multivariate/_base.py | 8 + src/process_improve/multivariate/_common.py | 25 ++ .../multivariate/_diagnostics.py | 252 ++++++++++++++++++ src/process_improve/multivariate/_mbpca.py | 133 ++++++--- src/process_improve/multivariate/_mbpls.py | 171 ++++++++---- src/process_improve/multivariate/_pca.py | 86 ------ src/process_improve/multivariate/_pls.py | 71 ----- src/process_improve/multivariate/methods.py | 4 + tests/test_multiblock_reference.py | 75 +++--- tests/test_multivariate.py | 196 ++++++++------ 11 files changed, 673 insertions(+), 352 deletions(-) diff --git a/src/process_improve/multivariate/__init__.py b/src/process_improve/multivariate/__init__.py index 0ba9fbd6..65edf243 100644 --- a/src/process_improve/multivariate/__init__.py +++ b/src/process_improve/multivariate/__init__.py @@ -10,11 +10,13 @@ MCUVScaler, center, eigenvalue_summary, + group_contributions, observation_contributions, project_variables, rv2_coefficient, rv_coefficient, scale, + score_contributions, spe_contributions, squared_cosine, t2_contributions, @@ -44,6 +46,7 @@ "correlation_loadings_plot", "eigenvalue_summary", "explained_variance_plot", + "group_contributions", "loading_plot", "observation_contributions", "predictions_vs_observed_plot", @@ -51,6 +54,7 @@ "rv2_coefficient", "rv_coefficient", "scale", + "score_contributions", "score_plot", "spe_contributions", "spe_plot", diff --git a/src/process_improve/multivariate/_base.py b/src/process_improve/multivariate/_base.py index 705bd2de..c60478d2 100644 --- a/src/process_improve/multivariate/_base.py +++ b/src/process_improve/multivariate/_base.py @@ -31,12 +31,18 @@ from ._diagnostics import ( eigenvalue_summary as _eigenvalue_summary, ) +from ._diagnostics import ( + group_contributions as _group_contributions, +) from ._diagnostics import ( observation_contributions as _observation_contributions, ) from ._diagnostics import ( project_variables as _project_variables, ) +from ._diagnostics import ( + score_contributions as _score_contributions, +) from ._diagnostics import ( spe_contributions as _spe_contributions, ) @@ -216,6 +222,8 @@ def __setstate__(self, state: dict) -> None: observation_contributions = _model_method(_observation_contributions) t2_contributions = _model_method(_t2_contributions) spe_contributions = _model_method(_spe_contributions) + score_contributions = _model_method(_score_contributions) + group_contributions = _model_method(_group_contributions) eigenvalue_summary = _model_method(_eigenvalue_summary) project_variables = _model_method(_project_variables) diff --git a/src/process_improve/multivariate/_common.py b/src/process_improve/multivariate/_common.py index 419295ed..dbdedc45 100644 --- a/src/process_improve/multivariate/_common.py +++ b/src/process_improve/multivariate/_common.py @@ -301,3 +301,28 @@ def method(self: object, *args, **kwargs) -> object: return fn(self, *args, **kwargs) return method + + +def _scale_block_contributions( + blocks: dict[str, np.ndarray], scaling: str +) -> dict[str, np.ndarray]: + """Apply a contribution-plot scaling jointly across every block. + + The scalings are those of Miller, Swanson and Heckler (1994). Both are taken + over the blocks together, not one block at a time: a super score pools all + the blocks, so scaling each block separately would make bars from different + blocks incomparable, which is the one thing a multi-block contribution plot + is read for. + """ + if scaling == "none": + return blocks + if scaling == "maximum": + largest = max((float(np.abs(v).max()) for v in blocks.values() if v.size), default=0.0) + scalar = largest if largest > epsqrt else 1.0 + return {name: values / scalar for name, values in blocks.items()} + if scaling == "within": + totals = sum(np.abs(values).sum(axis=1) for values in blocks.values()) + per_row = np.where(totals > epsqrt, totals, 1.0).reshape(-1, 1) + return {name: values / per_row for name, values in blocks.items()} + msg = f"scaling must be one of 'none', 'maximum' or 'within', got {scaling!r}." + raise ValueError(msg) diff --git a/src/process_improve/multivariate/_diagnostics.py b/src/process_improve/multivariate/_diagnostics.py index 1f7359b8..c4a60bc1 100644 --- a/src/process_improve/multivariate/_diagnostics.py +++ b/src/process_improve/multivariate/_diagnostics.py @@ -12,6 +12,8 @@ from __future__ import annotations +from collections.abc import Sequence + import numpy as np import pandas as pd from scipy.stats import f as f_dist @@ -632,6 +634,256 @@ def spe_contributions(model: BaseEstimator, X: DataMatrix) -> pd.DataFrame: return pd.DataFrame(residuals, index=X_df.index, columns=X_df.columns) +def _score_contribution_terms( + model: BaseEstimator, X: DataMatrix, component: int +) -> tuple[pd.DataFrame, np.ndarray]: + """Return ``(X_aligned, r_a)`` for the requested 1-based ``component``.""" + X_df, R, _ = _contribution_inputs(model, X) + A = R.shape[1] + if not (1 <= int(component) <= A): + msg = f"component must be a 1-based index within 1..{A}, got {component}." + raise ValueError(msg) + return X_df, R[:, int(component) - 1] + + +_OLD_SCORE_CONTRIBUTIONS_API = """\ +score_contributions() now takes the preprocessed data X, not a score vector. + +Until version {version} it back-projected a score-space difference through the +loadings, which does not depend on the observation's data at all: every +observation produced the same ranking of variables, namely the ranking of the +loadings. A contribution has to be x_ik * R_ka so that the per-variable terms +sum to the score. See Miller, Swanson and Heckler (1994). + + old: model.score_contributions(model.scores_.iloc[i]) + new: model.score_contributions(X) # all observations, component 1 + model.score_contributions(X).iloc[i] # just observation i + + old: model.score_contributions(t_a, t_b) # two observations + new: model.group_contributions(X, group=[a], reference=[b]) + + old: model.score_contributions(t, weighted=True) + new: model.t2_contributions(X) # already correct; unchanged +""" + + +_OLD_SCORE_CONTRIBUTIONS_KEYWORDS = frozenset({"t_end", "components", "weighted"}) + + +def _reject_old_score_contributions_api(X: object, keywords: dict) -> None: + """Raise a migration error if called with the pre-1.53.0 signature.""" + looks_like_scores = isinstance(X, (pd.Series, np.ndarray, list, tuple)) and np.ndim(X) == 1 + if looks_like_scores or (set(keywords) & _OLD_SCORE_CONTRIBUTIONS_KEYWORDS): + raise TypeError(_OLD_SCORE_CONTRIBUTIONS_API.format(version="1.53.0")) + if keywords: + unexpected = ", ".join(sorted(keywords)) + msg = f"score_contributions() got an unexpected keyword argument: {unexpected}." + raise TypeError(msg) + + +def score_contributions( + model: BaseEstimator, + X: DataMatrix, + component: int = 1, + scaling: str = "none", + **deprecated: object, +) -> pd.DataFrame: + r"""Per-variable contributions to a single score, :math:`t_a`. + + Works with fitted :class:`PCA` and :class:`PLS` models. A score is a + weighted sum of the (preprocessed) variables, so it splits exactly into one + term per variable. The contribution of variable :math:`k` to the score of + observation :math:`i` on component :math:`a` is + + .. math:: + + c_{ik}^{(a)} = x_{ik}\, R_{ka}, \qquad + \sum_{k=1}^{K} c_{ik}^{(a)} = t_{ia}, + + where :math:`R` is the score-generating matrix (``loadings_`` for PCA, + ``direct_weights_`` for PLS, so that :math:`T = XR`). This is the + contribution of Miller, Swanson and Heckler (1994); the generalisation from + PCA loadings to any latent-variable model's score-generating weights follows + Westerhuis, Gurden and Smilde (2000). + + The distinction from a loading plot is the point of the diagnostic. A + loading :math:`R_{ka}` describes the whole data set; a contribution + :math:`x_{ik} R_{ka}` describes *one observation*, and a variable with a + large loading contributes nothing when that observation sits at its mean. + Ranking variables by loading can therefore point at a different cause than + ranking them by contribution. + + Parameters + ---------- + model : PCA or PLS + A fitted PCA or PLS model. + X : array-like of shape (n_samples, n_features) + Preprocessed data, scaled the same way as the training data (for + example with :class:`MCUVScaler`). Passing the training data reproduces + the model's stored ``scores_``. + component : int, default=1 + **1-based** component index whose score is decomposed, matching the + model's column convention. + scaling : {"none", "maximum", "within"}, default="none" + Presentation scaling from Miller, Swanson and Heckler (1994). ``"none"`` + returns the raw contributions, which sum to the score. ``"maximum"`` + divides by the largest absolute contribution anywhere in ``X``, so a bar + of :math:`\pm 1` marks the most extreme variable-observation pair in the + data set. ``"within"`` divides each row by the sum of its absolute + contributions, so each row is on a common footing. Both scalings leave + the *pattern* of bars within a row unchanged; neither preserves the sum + to the score. + **deprecated + Captures ``t_end``, ``components`` and ``weighted`` from the pre-1.53.0 + calling convention, which took a score vector instead of ``X``, so that + those calls raise a :class:`TypeError` explaining the change rather than + quietly returning something else. Passing a 1-D ``X`` raises the same + error. + + Returns + ------- + pd.DataFrame + Signed contributions of shape (n_samples, n_features). With the default + ``scaling="none"``, each row sums to that observation's score on the + selected component. + + Examples + -------- + >>> pca = PCA(n_components=2).fit(X_scaled) + >>> contrib = pca.score_contributions(X_scaled, component=1) + >>> contrib.sum(axis=1) # equals pca.scores_[1] + >>> contrib.loc["33"].abs().sort_values() # what makes observation 33 extreme + + References + ---------- + Miller, P., Swanson, R.E. and Heckler, C.E. (1994). "Contribution plots: a + missing link in multivariate quality control." Applied Mathematics and + Computer Science, 8(4), 775-792. + + Westerhuis, J.A., Gurden, S.P. and Smilde, A.K. (2000). "Generalized + contribution plots in multivariate statistical process monitoring." + Chemometrics and Intelligent Laboratory Systems, 51(1), 95-114. + + See Also + -------- + group_contributions : The same decomposition for a group of observations, or + for the difference between two groups. + t2_contributions : Decomposes Hotelling's :math:`T^2`, which pools all + components rather than reading one at a time. + spe_contributions : The residual-space counterpart. + """ + _reject_old_score_contributions_api(X, deprecated) + X_df, r_a = _score_contribution_terms(model, X, component) + contributions = X_df.to_numpy(dtype=float) * r_a + + if scaling == "maximum": + largest = float(np.abs(contributions).max()) + contributions = contributions / (largest if largest > epsqrt else 1.0) + elif scaling == "within": + row_total = np.abs(contributions).sum(axis=1, keepdims=True) + contributions = contributions / np.where(row_total > epsqrt, row_total, 1.0) + elif scaling != "none": + msg = f"scaling must be one of 'none', 'maximum' or 'within', got {scaling!r}." + raise ValueError(msg) + + return pd.DataFrame(contributions, index=X_df.index, columns=X_df.columns) + + +def group_contributions( + model: BaseEstimator, + X: DataMatrix, + group: Sequence, + reference: Sequence | None = None, + component: int = 1, +) -> pd.Series: + r"""Per-variable contributions to a group's average score, or to a shift. + + The group form of :func:`score_contributions`. Averaging the data rows + before multiplying by the score-generating weights answers "what do these + observations have in common?" rather than "why is this one observation + unusual?", which is the question a cluster on a score plot, or a level shift + part-way through a data set, actually poses. + + With ``reference=None`` the contributions are + + .. math:: + + c_k = \bar{x}_{k}^{(G)} R_{ka}, \qquad \sum_k c_k = \bar{t}_a^{(G)}, + + comparing the group mean against the model centre. Given a ``reference`` + group the contributions are :math:`(\bar{x}_k^{(G)} - \bar{x}_k^{(H)}) + R_{ka}`, which sum to the difference in average score between the two + groups. Both forms are the linear-combination generalisation given by + Miller, Swanson and Heckler (1994). + + Parameters + ---------- + model : PCA or PLS + A fitted PCA or PLS model. + X : array-like of shape (n_samples, n_features) + Preprocessed data, scaled the same way as the training data. + group : sequence + Index labels (or a boolean mask, or positional integers) selecting the + observations of interest. + reference : sequence, optional + Index labels selecting the observations to compare against. ``None`` + (default) compares the group against the model centre. + component : int, default=1 + **1-based** component index whose score is decomposed. + + Returns + ------- + pd.Series + Signed contributions, one per variable. Sums to the group's average + score, or to the difference in average score between the two groups. + + Examples + -------- + >>> pca = PCA(n_components=2).fit(X_scaled) + >>> # Five batches that cluster together on the score plot: + >>> pca.group_contributions(X_scaled, group=[31, 142, 147, 220, 221]) + >>> # What shifted at batch 74? + >>> pca.group_contributions(X_scaled, group=range(64, 74), reference=range(74, 84)) + + See Also + -------- + score_contributions : The single-observation form. + """ + X_df, r_a = _score_contribution_terms(model, X, component) + + def _mean_row(selector: Sequence, name: str) -> np.ndarray: + selected = _select_rows(X_df, selector, name) + if len(selected) == 0: + msg = f"{name} selected no observations from X." + raise ValueError(msg) + return selected.to_numpy(dtype=float).mean(axis=0) + + deviation = _mean_row(group, "group") + if reference is not None: + deviation = deviation - _mean_row(reference, "reference") + + return pd.Series(deviation * r_a, index=X_df.columns, name="group_contributions") + + +def _select_rows(X_df: pd.DataFrame, selector: Sequence, name: str) -> pd.DataFrame: + """Select rows of ``X_df`` by label, boolean mask, or positional index.""" + values = list(selector) + if values and all(isinstance(v, (bool, np.bool_)) for v in values): + if len(values) != len(X_df): + msg = f"{name} is a boolean mask of length {len(values)}, but X has {len(X_df)} rows." + raise ValueError(msg) + return X_df.loc[np.asarray(values, dtype=bool)] + missing = [v for v in values if v not in X_df.index] + if not missing: + return X_df.loc[values] + # Fall back to positional selection, but only if every entry is a valid + # position; a partial match means the caller mixed labels and positions. + if all(isinstance(v, (int, np.integer)) and 0 <= int(v) < len(X_df) for v in values): + return X_df.iloc[[int(v) for v in values]] + msg = f"{name} contains entries that are neither index labels of X nor valid row positions: {missing}." + raise ValueError(msg) + + def eigenvalue_summary(model: BaseEstimator) -> pd.DataFrame: """Summarize the variance captured by each component as a tidy table. diff --git a/src/process_improve/multivariate/_mbpca.py b/src/process_improve/multivariate/_mbpca.py index fb15e259..394ce4c5 100644 --- a/src/process_improve/multivariate/_mbpca.py +++ b/src/process_improve/multivariate/_mbpca.py @@ -10,6 +10,7 @@ import time import typing import warnings +from collections.abc import Sequence import numpy as np import pandas as pd @@ -18,7 +19,8 @@ from sklearn.utils.validation import check_is_fitted from ._base import _HotellingsT2LimitMixin -from ._common import SpecificationWarning, _nz, epsqrt +from ._common import SpecificationWarning, _nz, _scale_block_contributions, epsqrt +from ._diagnostics import _select_rows from ._limits import spe_calculation from ._nipals import quick_regress, ssq from ._preprocessing import MCUVScaler @@ -585,49 +587,120 @@ def spe_contributions(self, X: dict[str, pd.DataFrame]) -> dict[str, pd.DataFram out[name] = pd.DataFrame(residuals_sq, index=sample_index, columns=self._block_columns[name]) return out + def _deflated_blocks( + self, X: dict[str, pd.DataFrame], component: int + ) -> tuple[dict[str, np.ndarray], pd.Index]: + """Preprocessed blocks, deflated through the first ``component - 1`` components. + + The super score at component *a* is formed from the data that remain + after the earlier components have been removed, so a decomposition of + that score has to start from the same deflated data. + """ + if not isinstance(X, dict): + raise TypeError("X must be a dict[str, pd.DataFrame].") + missing = set(self.block_names_) - set(X) + if missing: + raise ValueError(f"Missing X-blocks: {sorted(missing)}.") + a_max = int(self.n_components) + if not (1 <= int(component) <= a_max): + msg = f"component must be a 1-based index within 1..{a_max}, got {component}." + raise ValueError(msg) + + sample_index: pd.Index | None = None + x_def: dict[str, np.ndarray] = {} + for name in self.block_names_: + block = X[name] + if not isinstance(block, pd.DataFrame): + block = pd.DataFrame(block, columns=self._block_columns[name]) + if sample_index is None: + sample_index = block.index + x_def[name] = self.preproc_[name].transform(block).values.astype(float) + + sqrt_kb = {name: float(np.sqrt(self.block_widths_[name])) for name in self.block_names_} + for a in range(int(component) - 1): + t_b_row = np.column_stack([ + x_def[name] + @ self.block_loadings_[name].values[:, a] + / _nz(float(self.block_loadings_[name].values[:, a] @ self.block_loadings_[name].values[:, a])) + / sqrt_kb[name] + for name in self.block_names_ + ]) + p_s = self.super_loadings_.values[:, a] + t_super = t_b_row @ p_s / _nz(float(p_s @ p_s)) + for b_idx, name in enumerate(self.block_names_): + p_b = self.block_loadings_[name].values[:, a] + x_def[name] = x_def[name] - np.outer(t_super, p_b * p_s[b_idx] * sqrt_kb[name]) + + assert sample_index is not None + return x_def, sample_index + def score_contributions( self, - t_super_start: np.ndarray | pd.Series, - t_super_end: np.ndarray | pd.Series | None = None, - components: list[int] | None = None, - *, - weighted: bool = False, - ) -> dict[str, pd.Series]: - r"""Per-block per-variable contributions to a super-score movement (MBPCA). + X: dict[str, pd.DataFrame], + component: int = 1, + scaling: str = "none", + ) -> dict[str, pd.DataFrame]: + r"""Per-block per-variable contributions to a super-score (MBPCA). - Decomposes a super-score-space delta into preprocessed-scale variable - contributions per X-block. The MBPCA back-projection mirrors the - deflation step used during fit: + The multi-block analogue of :meth:`PCA.score_contributions`. A super + score is a weighted sum of the (deflated, preprocessed) variables across + every block, so it splits exactly into one term per variable: .. math:: - \text{contrib}_{b,j} = \sum_a (\Delta t_\mathrm{super}[a]) \cdot - P_b[j, a] \cdot p_\mathrm{super}[b, a] \cdot \sqrt{K_b} + c_{b,ij}^{(a)} = \tilde{x}_{b,ij}^{(a)}\, + \frac{P_b[j, a]\, p_\mathrm{super}[b, a]} + {(p_b^\top p_b)(p_\mathrm{super}^\top p_\mathrm{super})\sqrt{K_b}}, + \qquad + \sum_b \sum_j c_{b,ij}^{(a)} = t_{\mathrm{super},ia}, + + where :math:`\tilde{x}^{(a)}` is the block data deflated through the + first :math:`a-1` components. See :meth:`MBPLS.score_contributions` for the parameter and return descriptions; the API is identical. """ check_is_fitted(self, "block_loadings_") - t_start = np.asarray(t_super_start, dtype=float) - t_end = np.zeros(self.n_components) if t_super_end is None else np.asarray(t_super_end, dtype=float) - idx = np.arange(self.n_components) if components is None else np.array(components) - 1 - dt = t_end[idx] - t_start[idx] - if weighted: - # ``explained_variance_[a] == 0`` for a degenerate component - # would silently produce inf/NaN weighted contributions. - # Clamp the divisor so weighting is a no-op on such components. - # SEC-21 (#270) sub-item 5. - ev = np.asarray(self.explained_variance_)[idx] - dt = dt / np.sqrt(np.where(ev > epsqrt, ev, 1.0)) + deflated, sample_index = self._deflated_blocks(X, component) + a = int(component) - 1 + p_s = self.super_loadings_.values[:, a] + super_norm = _nz(float(p_s @ p_s)) - out: dict[str, pd.Series] = {} + raw: dict[str, np.ndarray] = {} for b_idx, name in enumerate(self.block_names_): sqrt_kb = float(np.sqrt(self.block_widths_[name])) - ps = self.super_loadings_.values[b_idx, idx] # (len(idx),) - pb = self.block_loadings_[name].values[:, idx] # (K_b, len(idx)) - effective = pb * (ps * sqrt_kb) # (K_b, len(idx)) - contrib = effective @ dt # (K_b,) - out[name] = pd.Series(contrib, index=self._block_columns[name], name=f"score_contributions[{name}]") + p_b = self.block_loadings_[name].values[:, a] + weight = p_b / _nz(float(p_b @ p_b)) / sqrt_kb * (p_s[b_idx] / super_norm) + raw[name] = deflated[name] * weight + + raw = _scale_block_contributions(raw, scaling) + return { + name: pd.DataFrame(values, index=sample_index, columns=self._block_columns[name]) + for name, values in raw.items() + } + + def group_contributions( + self, + X: dict[str, pd.DataFrame], + group: Sequence, + reference: Sequence | None = None, + component: int = 1, + ) -> dict[str, pd.Series]: + """Per-block per-variable contributions to a group's average super score. + + The multi-block analogue of :meth:`PCA.group_contributions`. See that + method for the definition; the result is returned one Series per + X-block, and the sum over every block equals the group's average super + score (or the difference between the two groups' average super scores + when ``reference`` is given). + """ + per_block = self.score_contributions(X, component=component) + out: dict[str, pd.Series] = {} + for name, frame in per_block.items(): + deviation = _select_rows(frame, group, "group").mean(axis=0) + if reference is not None: + deviation = deviation - _select_rows(frame, reference, "reference").mean(axis=0) + out[name] = pd.Series(deviation, name=f"group_contributions[{name}]") return out def super_score_plot(self, pc_horiz: int = 1, pc_vert: int = 2) -> go.Figure: diff --git a/src/process_improve/multivariate/_mbpls.py b/src/process_improve/multivariate/_mbpls.py index 4b4e3c5f..5ceb92c2 100644 --- a/src/process_improve/multivariate/_mbpls.py +++ b/src/process_improve/multivariate/_mbpls.py @@ -12,6 +12,7 @@ import time import typing import warnings +from collections.abc import Sequence import numpy as np import pandas as pd @@ -21,7 +22,8 @@ from ..visualization.themes import REFERENCE_LINE_COLOR from ._base import _HotellingsT2LimitMixin -from ._common import SpecificationWarning, _nz, epsqrt +from ._common import SpecificationWarning, _nz, _scale_block_contributions, epsqrt +from ._diagnostics import _select_rows from ._limits import spe_calculation from ._nipals import quick_regress, ssq from ._preprocessing import MCUVScaler @@ -620,70 +622,139 @@ def spe_contributions(self, X: dict[str, pd.DataFrame]) -> dict[str, pd.DataFram out[name] = pd.DataFrame(residuals_sq, index=sample_index, columns=self._block_columns[name]) return out + def _deflated_blocks( + self, X: dict[str, pd.DataFrame], component: int + ) -> tuple[dict[str, np.ndarray], pd.Index]: + """Preprocessed blocks, deflated through the first ``component - 1`` components. + + The super score at component *a* is formed from the data that remain + after the earlier components have been removed, so a decomposition of + that score has to start from the same deflated data. + """ + if not isinstance(X, dict): + raise TypeError("X must be a dict[str, pd.DataFrame].") + missing = set(self.block_names_) - set(X) + if missing: + raise ValueError(f"Missing X-blocks: {sorted(missing)}.") + a_max = int(self.n_components) + if not (1 <= int(component) <= a_max): + msg = f"component must be a 1-based index within 1..{a_max}, got {component}." + raise ValueError(msg) + + sample_index: pd.Index | None = None + x_def: dict[str, np.ndarray] = {} + for name in self.block_names_: + block = X[name] + if not isinstance(block, pd.DataFrame): + block = pd.DataFrame(block, columns=self._block_columns[name]) + if sample_index is None: + sample_index = block.index + x_def[name] = self.preproc_[name].transform(block).values.astype(float) + + sqrt_kb = {name: float(np.sqrt(self.block_widths_[name])) for name in self.block_names_} + for a in range(int(component) - 1): + t_b_row = np.column_stack([ + x_def[name] @ self.block_weights_[name].values[:, a] / sqrt_kb[name] + for name in self.block_names_ + ]) + t_super = t_b_row @ self.super_weights_.values[:, a] + for name in self.block_names_: + p_b = self.block_loadings_[name].values[:, a] + x_def[name] = x_def[name] - np.outer(t_super, p_b) + + assert sample_index is not None + return x_def, sample_index + def score_contributions( self, - t_super_start: np.ndarray | pd.Series, - t_super_end: np.ndarray | pd.Series | None = None, - components: list[int] | None = None, - *, - weighted: bool = False, - ) -> dict[str, pd.Series]: - r"""Per-block per-variable contributions to a super-score movement. + X: dict[str, pd.DataFrame], + component: int = 1, + scaling: str = "none", + ) -> dict[str, pd.DataFrame]: + r"""Per-block per-variable contributions to a super-score. + + The multi-block analogue of :meth:`PLS.score_contributions`. A super + score is a weighted sum of the (deflated, preprocessed) variables across + every block, so it splits exactly into one term per variable: + + .. math:: - The multi-block analogue of :meth:`PCA.score_contributions`. Decomposes - a super-score-space delta back into preprocessed-scale variable-space - deltas, one per X-block. + c_{b,ij}^{(a)} = \tilde{x}_{b,ij}^{(a)}\, + \frac{w_b[j, a]\, w_\mathrm{super}[b, a]}{\sqrt{K_b}}, + \qquad + \sum_b \sum_j c_{b,ij}^{(a)} = t_{\mathrm{super},ia}, - For MBPLS, the super-score at component *a* is built from the per-block - scores (which themselves are weighted regressions of the block on - ``w_b``), so the back-projection through ``w_super[b,a] * w_b[:,a] / - sqrt(K_b)`` gives the variable contribution. + where :math:`\tilde{x}^{(a)}` is the block data deflated through the + first :math:`a-1` components, which is what the super score at component + :math:`a` is actually formed from. Parameters ---------- - t_super_start : array-like, shape (n_components,) - Super-score row of the observation of interest. Typically a row - from ``self.super_scores_`` or from ``predict(X_new).super_scores``. - t_super_end : array-like, optional - Reference point in super-score space. Defaults to the model - centre (zeros). - components : list of int, optional - **1-based** component indices to decompose over. ``None`` (default) - uses all components - appropriate for Hotelling's T² contributions. - weighted : bool, default=False - If ``True``, divide the super-score delta by - ``sqrt(explained_variance_)`` per component before back-projecting, - giving contributions to the T² statistic instead of the Euclidean - super-score distance. + X : dict[str, pd.DataFrame] + Raw (un-preprocessed) X-blocks, keyed by block name, exactly as + passed to :meth:`fit`. The stored per-block preprocessing is applied + internally. + component : int, default=1 + **1-based** component index whose super score is decomposed. + scaling : {"none", "maximum", "within"}, default="none" + Presentation scaling, as for :meth:`PLS.score_contributions`. Under + ``"none"`` the contributions sum across all blocks to the super + score. ``"maximum"`` divides by the largest absolute contribution + over every block; ``"within"`` divides each observation by the total + absolute contribution it accumulates across every block, so both + scalings are taken over the blocks jointly rather than one block at + a time. Returns ------- - dict[str, pd.Series] - One Series per X-block (length ``K_b``), indexed by variable - (column) name. + dict[str, pd.DataFrame] + One frame per X-block, of shape (n_samples, K_b). + + Examples + -------- + >>> mbpls = MBPLS(n_components=2).fit(blocks, Y) + >>> contrib = mbpls.score_contributions(blocks, component=1) + >>> sum(frame.sum(axis=1) for frame in contrib.values()) # super score 1 """ check_is_fitted(self, "block_weights_") - t_start = np.asarray(t_super_start, dtype=float) - t_end = np.zeros(self.n_components) if t_super_end is None else np.asarray(t_super_end, dtype=float) - idx = np.arange(self.n_components) if components is None else np.array(components) - 1 - dt = t_end[idx] - t_start[idx] # (len(idx),) - if weighted: - # ``explained_variance_[a] == 0`` for a degenerate component - # would silently produce inf/NaN weighted contributions. - # Clamp the divisor so weighting is a no-op on such components. - # SEC-21 (#270) sub-item 5. - ev = np.asarray(self.explained_variance_)[idx] - dt = dt / np.sqrt(np.where(ev > epsqrt, ev, 1.0)) + deflated, sample_index = self._deflated_blocks(X, component) + a = int(component) - 1 - out: dict[str, pd.Series] = {} + raw: dict[str, np.ndarray] = {} for b_idx, name in enumerate(self.block_names_): sqrt_kb = float(np.sqrt(self.block_widths_[name])) - ws = self.super_weights_.values[b_idx, idx] # (len(idx),) - wb = self.block_weights_[name].values[:, idx] # (K_b, len(idx)) - # Effective per-component contribution per variable: w_super[b] * w_b[:,a] / sqrt(K_b) - effective = wb * (ws / sqrt_kb) # (K_b, len(idx)) - contrib = effective @ dt # (K_b,) - out[name] = pd.Series(contrib, index=self._block_columns[name], name=f"score_contributions[{name}]") + w_b = self.block_weights_[name].values[:, a] + w_s = float(self.super_weights_.values[b_idx, a]) + raw[name] = deflated[name] * (w_b * w_s / sqrt_kb) + + raw = _scale_block_contributions(raw, scaling) + return { + name: pd.DataFrame(values, index=sample_index, columns=self._block_columns[name]) + for name, values in raw.items() + } + + def group_contributions( + self, + X: dict[str, pd.DataFrame], + group: Sequence, + reference: Sequence | None = None, + component: int = 1, + ) -> dict[str, pd.Series]: + """Per-block per-variable contributions to a group's average super score. + + The multi-block analogue of :meth:`PLS.group_contributions`. See that + method for the definition; the only difference is that the result is + returned one Series per X-block, and the sum over every block equals the + group's average super score (or the difference between the two groups' + average super scores when ``reference`` is given). + """ + per_block = self.score_contributions(X, component=component) + out: dict[str, pd.Series] = {} + for name, frame in per_block.items(): + deviation = _select_rows(frame, group, "group").mean(axis=0) + if reference is not None: + deviation = deviation - _select_rows(frame, reference, "reference").mean(axis=0) + out[name] = pd.Series(deviation, name=f"group_contributions[{name}]") return out def super_score_plot(self, pc_horiz: int = 1, pc_vert: int = 2) -> go.Figure: diff --git a/src/process_improve/multivariate/_pca.py b/src/process_improve/multivariate/_pca.py index 69189e45..ca6cb4b2 100644 --- a/src/process_improve/multivariate/_pca.py +++ b/src/process_improve/multivariate/_pca.py @@ -1281,92 +1281,6 @@ def select_n_components( # noqa: PLR0913, PLR0915, C901 **consensus_fields, ) - def score_contributions( - self, - t_start: np.ndarray | pd.Series, - t_end: np.ndarray | pd.Series | None = None, - components: list[int] | None = None, - *, - weighted: bool = False, - ) -> pd.Series: - """Contribution of each variable to a score-space movement. - - Decomposes the difference (t_end - t_start) back into variable space - via the loadings for the selected components. - - Parameters - ---------- - t_start : array-like of shape (n_components,) - Score vector of the observation of interest. Can be a row from - ``self.scores_`` or from ``predict(X_new).scores``. - t_end : array-like of shape (n_components,), optional - Reference point in score space. Default is the model center - (a vector of zeros), which is the most common choice. - components : list of int, optional - **1-based** component indices to decompose over, matching the - model's column convention. Examples: ``[2, 3]`` for a PC2-vs-PC3 - score plot, or ``None`` (default) for all components - appropriate - for Hotelling's T² contributions. - weighted : bool, default False - If True, scale the score difference by 1/sqrt(explained_variance) - per component before back-projecting. This gives contributions to - the T² statistic rather than to the Euclidean score distance. - - Returns - ------- - contributions : pd.Series of shape (n_features,) - One value per variable; sign indicates direction. Index contains - the variable (feature) names. - - Examples - -------- - >>> pca = PCA(n_components=3).fit(X_scaled) - >>> # Why does observation 0 differ from the model center? - >>> contrib = pca.score_contributions(pca.scores_.iloc[0].values) - >>> # Which variables drive the difference between obs 0 and obs 5? - >>> contrib = pca.score_contributions( - ... pca.scores_.iloc[0].values, pca.scores_.iloc[5].values - ... ) - - See Also - -------- - observation_contributions : The per-observation counterpart. Despite - the similar name, the two answer different questions and are not - interchangeable. ``score_contributions`` is *per-variable*: it - decomposes a single observation's movement in score space back - onto the original variables, answering "which **variables** - explain why this observation sits where it does?". It returns one - signed value per variable. ``observation_contributions`` is - *per-observation*: it reports each observation's share of a - component's total inertia, answering "which **observations** most - strongly shape this component?". It returns a non-negative - sample-by-component table whose columns each sum to 1. The two are - orthogonal views of the same score matrix - one decomposes across - variables, the other across observations. - """ - check_is_fitted(self, "loadings_") - t_start = np.asarray(t_start, dtype=float) - t_end = np.zeros(self.n_components) if t_end is None else np.asarray(t_end, dtype=float) - - idx = ( - np.arange(self.n_components) if components is None else np.array(components) - 1 - ) # convert 1-based to 0-based - - dt = t_end[idx] - t_start[idx] - - if weighted: - # ``explained_variance_[a] == 0`` for a degenerate component - # would silently produce inf/NaN weighted contributions. - # Clamp the divisor so weighting is a no-op on such components. - # SEC-21 (#270) sub-item 5. - ev = np.asarray(self.explained_variance_)[idx] - dt = dt / np.sqrt(np.where(ev > epsqrt, ev, 1.0)) - - P = self._loadings[:, idx].T # (len(idx), n_features) - contributions = dt @ P # (n_features,) - - return pd.Series(contributions, index=self._feature_names, name="score_contributions") - def detect_outliers(self, conf_level: float = 0.95) -> list[dict]: """Detect outlier observations using SPE and Hotelling's T² diagnostics. diff --git a/src/process_improve/multivariate/_pls.py b/src/process_improve/multivariate/_pls.py index f6045d4c..9cd41b36 100644 --- a/src/process_improve/multivariate/_pls.py +++ b/src/process_improve/multivariate/_pls.py @@ -1838,77 +1838,6 @@ def nested_cv( # noqa: PLR0913, PLR0915 selected_components_distribution=distribution, ) - def score_contributions( - self, - t_start: np.ndarray | pd.Series, - t_end: np.ndarray | pd.Series | None = None, - components: list[int] | None = None, - *, - weighted: bool = False, - ) -> pd.Series: - """Contribution of each X variable to a score-space movement. - - Identical to ``PCA.score_contributions`` but uses the PLS X-loadings - (P matrix) for back-projection. - - Parameters - ---------- - t_start : array-like of shape (n_components,) - Score vector of the observation of interest. - t_end : array-like of shape (n_components,), optional - Reference point in score space. Default is the model center (zeros). - components : list of int, optional - **1-based** component indices. Default: all components. - weighted : bool, default False - If True, scale by 1/sqrt(explained_variance) for T² contributions. - - Returns - ------- - contributions : pd.Series of shape (n_features,) - - Examples - -------- - >>> pls = PLS(n_components=3).fit(X_scaled, Y_scaled) - >>> contrib = pls.score_contributions(pls.scores_.iloc[0].values) - >>> contrib.abs().sort_values(ascending=False).head() # top contributors - - See Also - -------- - observation_contributions : The per-observation counterpart. Despite - the similar name, the two answer different questions and are not - interchangeable. ``score_contributions`` is *per-variable*: it - decomposes a single observation's movement in score space back - onto the original X variables, answering "which **variables** - explain why this observation sits where it does?". It returns one - signed value per variable. ``observation_contributions`` is - *per-observation*: it reports each observation's share of a - component's total inertia, answering "which **observations** most - strongly shape this component?". It returns a non-negative - sample-by-component table whose columns each sum to 1. The two are - orthogonal views of the same score matrix - one decomposes across - variables, the other across observations. - """ - check_is_fitted(self, "x_loadings_") - t_start = np.asarray(t_start, dtype=float) - t_end = np.zeros(self.n_components) if t_end is None else np.asarray(t_end, dtype=float) - - idx = np.arange(self.n_components) if components is None else np.array(components) - 1 - - dt = t_end[idx] - t_start[idx] - - if weighted: - # ``explained_variance_[a] == 0`` for a degenerate component - # would silently produce inf/NaN weighted contributions. - # Clamp the divisor so weighting is a no-op on such components. - # SEC-21 (#270) sub-item 5. - ev = np.asarray(self.explained_variance_)[idx] - dt = dt / np.sqrt(np.where(ev > epsqrt, ev, 1.0)) - - P = self._x_loadings[:, idx].T - contributions = dt @ P - - return pd.Series(contributions, index=self._feature_names, name="score_contributions") - def detect_outliers(self, conf_level: float = 0.95) -> list[dict]: """Detect outlier observations using SPE and Hotelling's T² diagnostics. diff --git a/src/process_improve/multivariate/methods.py b/src/process_improve/multivariate/methods.py index d90dfcc8..1f069874 100644 --- a/src/process_improve/multivariate/methods.py +++ b/src/process_improve/multivariate/methods.py @@ -20,10 +20,12 @@ from ._common import NotEnoughVarianceError, SpecificationWarning, epsqrt from ._diagnostics import ( eigenvalue_summary, + group_contributions, observation_contributions, project_variables, rv2_coefficient, rv_coefficient, + score_contributions, selectivity_ratio, spe_contributions, squared_cosine, @@ -91,6 +93,7 @@ "ellipse_coordinates", "epsqrt", "explained_variance_plot", + "group_contributions", "hotellings_t2_limit", "internal_pls_nipals_fit_one_pc", "loading_plot", @@ -105,6 +108,7 @@ "rv_coefficient", "safe_inverse", "scale", + "score_contributions", "score_limit", "score_plot", "selectivity_ratio", diff --git a/tests/test_multiblock_reference.py b/tests/test_multiblock_reference.py index b99e2fb6..17d689b0 100644 --- a/tests/test_multiblock_reference.py +++ b/tests/test_multiblock_reference.py @@ -379,25 +379,33 @@ def test_spe_contributions_sum_to_block_spe_squared(self, synthetic_two_block) - block_spe_sq = model.block_spe_[name].iloc[:, -1].values ** 2 np.testing.assert_array_almost_equal(row_sums, block_spe_sq, decimal=8) - def test_score_contributions_to_origin_back_project_correctly(self, synthetic_two_block) -> None: - """For ``t_super_end=0`` (the model centre), the per-block score - contributions must equal the per-block sums of effective-loading - projections of the (negated) super-score row. + def test_score_contributions_sum_to_the_super_score(self, synthetic_two_block) -> None: + """Summed over every block and variable, the contributions give the super score. + + The super score at component ``a`` is a weighted sum of the block data + deflated through the first ``a - 1`` components, so the decomposition + has to start from that same deflated data. """ from process_improve.multivariate.methods import MBPCA x_blocks, _ = synthetic_two_block model = MBPCA(n_components=2).fit(x_blocks) - t_row = model.super_scores_.iloc[0].values - contribs = model.score_contributions(t_row) - # Effective per-block reconstruction at the row level: compare manually - for b_idx, name in enumerate(model.block_names_): - sqrt_kb = float(np.sqrt(model.block_widths_[name])) - pb = model.block_loadings_[name].values - ps = model.super_loadings_.values[b_idx, :] - effective = pb * (ps * sqrt_kb) - expected = -t_row @ effective.T # back-projected variable contributions - np.testing.assert_array_almost_equal(contribs[name].values, expected, decimal=10) + for a in (1, 2): + contribs = model.score_contributions(x_blocks, component=a) + total = sum(frame.sum(axis=1) for frame in contribs.values()) + np.testing.assert_array_almost_equal( + total.values, model.super_scores_[a].values, decimal=10 + ) + + def test_group_contributions_sum_to_the_average_super_score(self, synthetic_two_block) -> None: + from process_improve.multivariate.methods import MBPCA + + x_blocks, _ = synthetic_two_block + model = MBPCA(n_components=2).fit(x_blocks) + group = list(range(5)) + contribs = model.group_contributions(x_blocks, group=group) + total = sum(float(series.sum()) for series in contribs.values()) + assert total == pytest.approx(model.super_scores_[1].values[group].mean(), abs=1e-10) def test_super_score_and_loadings_plots_return_figures(self, synthetic_two_block) -> None: import plotly.graph_objects as go @@ -1125,35 +1133,34 @@ def test_spe_contributions_sum_to_block_spe_squared(self, ldpe) -> None: block_spe_sq = model.block_spe_[name].iloc[:, -1].values ** 2 np.testing.assert_array_almost_equal(row_sums, block_spe_sq, decimal=8) - def test_score_contributions_to_origin_back_project_correctly(self, ldpe) -> None: - """For each per-block contribution Series, manual back-projection - through ``w_super[b] * w_b / sqrt(K_b)`` reproduces the result. - """ + def test_score_contributions_sum_to_the_super_score(self, ldpe) -> None: + """On the LDPE blocks, contributions decompose the super score exactly.""" from process_improve.multivariate.methods import MBPLS x_blocks, y_df = ldpe model = MBPLS(n_components=3).fit(x_blocks, y_df) - t_row = model.super_scores_.iloc[0].values - contribs = model.score_contributions(t_row) - for b_idx, name in enumerate(model.block_names_): - sqrt_kb = float(np.sqrt(model.block_widths_[name])) - wb = model.block_weights_[name].values - ws = model.super_weights_.values[b_idx, :] - effective = wb * (ws / sqrt_kb) - expected = -t_row @ effective.T - np.testing.assert_array_almost_equal(contribs[name].values, expected, decimal=10) - - def test_score_contributions_subset_of_components(self, ldpe) -> None: + for a in (1, 2, 3): + contribs = model.score_contributions(x_blocks, component=a) + total = sum(frame.sum(axis=1) for frame in contribs.values()) + np.testing.assert_array_almost_equal( + total.values, model.super_scores_[a].values, decimal=9 + ) + + def test_score_contributions_differ_between_components(self, ldpe) -> None: + """Each component gets its own decomposition; they are not interchangeable.""" from process_improve.multivariate.methods import MBPLS x_blocks, y_df = ldpe model = MBPLS(n_components=3).fit(x_blocks, y_df) - t_row = model.super_scores_.iloc[0].values - # Components 1 + 2 should differ from full (1+2+3) decomposition - partial = model.score_contributions(t_row, components=[1, 2]) - full = model.score_contributions(t_row) + first = model.score_contributions(x_blocks, component=1) + second = model.score_contributions(x_blocks, component=2) for name in model.block_names_: - assert not np.allclose(partial[name].values, full[name].values) + assert not np.allclose(first[name].values, second[name].values) + + with pytest.raises(ValueError, match="1-based index"): + model.score_contributions(x_blocks, component=4) + with pytest.raises(ValueError, match="Missing X-blocks"): + model.score_contributions({model.block_names_[0]: x_blocks[model.block_names_[0]]}) def test_super_score_plot_returns_plotly_figure(self, ldpe) -> None: import plotly.graph_objects as go diff --git a/tests/test_multivariate.py b/tests/test_multivariate.py index ece114ec..cd7c93f0 100644 --- a/tests/test_multivariate.py +++ b/tests/test_multivariate.py @@ -921,53 +921,102 @@ def test_pca_select_n_components_consensus_mode() -> None: def test_pca_score_contributions() -> None: - """Test score_contributions method on a simple dataset.""" + """Contributions decompose the score exactly, one term per variable.""" rng = np.random.default_rng(42) X = pd.DataFrame(rng.standard_normal((30, 5)), columns=[f"V{i}" for i in range(1, 6)]) X = MCUVScaler().fit_transform(X) pca = PCA(n_components=3).fit(X) - # --- Basic shape and index --- - obs = pca.scores_.iloc[0] - contrib = pca.score_contributions(obs) - assert isinstance(contrib, pd.Series) - assert len(contrib) == 5 - assert list(contrib.index) == [f"V{i}" for i in range(1, 6)] - - # --- Conservation: sum of contributions equals the projected reconstruction --- - # For unweighted, all-component case: contributions = (0 - t) @ P.T row-wise - # and sum(contributions * p_k) should recover dt for each component - P = pca.loadings_.values # (K, A) - dt = -obs.values # t_end(0) - t_start - expected = dt @ P.T - assert contrib.values == pytest.approx(expected, abs=1e-12) - - # --- Specific components (1-based) --- - contrib_23 = pca.score_contributions(obs, components=[2, 3]) - dt_23 = -obs.values[[1, 2]] # 0-based indices 1, 2 - P_23 = P[:, [1, 2]].T - expected_23 = dt_23 @ P_23 - assert contrib_23.values == pytest.approx(expected_23, abs=1e-12) - - # --- Custom t_end --- - obs2 = pca.scores_.iloc[1] - contrib_pair = pca.score_contributions(obs, t_end=obs2) - dt_pair = obs2.values - obs.values - expected_pair = dt_pair @ P.T - assert contrib_pair.values == pytest.approx(expected_pair, abs=1e-12) - - # --- Weighted mode (T² contributions) --- - contrib_w = pca.score_contributions(obs, weighted=True) - dt_w = -obs.values / np.sqrt(pca.explained_variance_) - expected_w = dt_w @ P.T - assert contrib_w.values == pytest.approx(expected_w, abs=1e-12) - - # --- Weighted + specific components --- - contrib_w2 = pca.score_contributions(obs, components=[1], weighted=True) - dt_w1 = -obs.values[0] / np.sqrt(pca.explained_variance_[0]) - expected_w1 = dt_w1 * P[:, 0] - assert contrib_w2.values == pytest.approx(expected_w1, abs=1e-12) + contrib = pca.score_contributions(X, component=1) + assert isinstance(contrib, pd.DataFrame) + assert contrib.shape == (30, 5) + assert list(contrib.columns) == [f"V{i}" for i in range(1, 6)] + assert list(contrib.index) == list(X.index) + + # --- The defining property (Miller, Swanson and Heckler, 1994, eq. 3): + # the per-variable terms sum to the score they decompose. + for a in (1, 2, 3): + terms = pca.score_contributions(X, component=a) + assert terms.sum(axis=1).to_numpy() == pytest.approx(pca.scores_[a].to_numpy(), abs=1e-12) + # ... and each term is the data value times its loading. + expected = X.to_numpy() * pca.loadings_.to_numpy()[:, a - 1] + assert terms.to_numpy() == pytest.approx(expected, abs=1e-12) + + # --- The diagnosis must depend on the observation, not only on the model. + # A back-projection of the score through the loadings (what this method did + # before 1.53.0) gives every observation the same ranking of variables, + # which is the loading ranking; that is the defect this test guards. + rankings = {tuple(np.argsort(-contrib.iloc[i].abs().to_numpy())) for i in range(len(X))} + assert len(rankings) > 1 + + # --- A variable sitting at its mean contributes nothing, whatever its loading. + at_mean = X.copy() + at_mean.iloc[0, 2] = 0.0 + assert pca.score_contributions(at_mean, component=1).iloc[0, 2] == pytest.approx(0.0, abs=1e-12) + + # --- Scaling changes the units, not the pattern within a row. + for scaling in ("maximum", "within"): + scaled = pca.score_contributions(X, component=1, scaling=scaling) + ratio = scaled.iloc[0].to_numpy() / contrib.iloc[0].to_numpy() + assert ratio == pytest.approx(np.full(5, ratio[0]), rel=1e-10) + assert np.abs(pca.score_contributions(X, component=1, scaling="maximum").to_numpy()).max() == pytest.approx(1.0) + + with pytest.raises(ValueError, match="scaling must be one of"): + pca.score_contributions(X, scaling="loud") + with pytest.raises(ValueError, match="1-based index"): + pca.score_contributions(X, component=4) + + +def test_score_contributions_rejects_the_old_api() -> None: + """The pre-1.53.0 signature took a score vector and must not silently work.""" + rng = np.random.default_rng(3) + X = MCUVScaler().fit_transform( + pd.DataFrame(rng.standard_normal((25, 4)), columns=list("abcd")) + ) + pca = PCA(n_components=2).fit(X) + + with pytest.raises(TypeError, match="not a score vector"): + pca.score_contributions(pca.scores_.iloc[0]) + with pytest.raises(TypeError, match="not a score vector"): + pca.score_contributions(pca.scores_.iloc[0].to_numpy()) + with pytest.raises(TypeError, match="t2_contributions"): + pca.score_contributions(X, weighted=True) + with pytest.raises(TypeError, match="group_contributions"): + pca.score_contributions(X, t_end=pca.scores_.iloc[1]) + + +def test_pca_group_contributions() -> None: + """Group contributions sum to the average score, or to the shift between groups.""" + rng = np.random.default_rng(11) + X = MCUVScaler().fit_transform( + pd.DataFrame(rng.standard_normal((40, 5)), columns=[f"V{i}" for i in range(1, 6)]) + ) + pca = PCA(n_components=2).fit(X) + early, late = list(X.index[:8]), list(X.index[20:30]) + + against_centre = pca.group_contributions(X, group=early) + assert isinstance(against_centre, pd.Series) + assert list(against_centre.index) == list(X.columns) + assert against_centre.sum() == pytest.approx(pca.scores_.loc[early, 1].mean(), abs=1e-12) + + shift = pca.group_contributions(X, group=early, reference=late, component=2) + expected = pca.scores_.loc[early, 2].mean() - pca.scores_.loc[late, 2].mean() + assert shift.sum() == pytest.approx(expected, abs=1e-12) + + # A single-observation group is the corresponding row of score_contributions. + one = pca.group_contributions(X, group=[X.index[3]]) + assert one.to_numpy() == pytest.approx( + pca.score_contributions(X, component=1).iloc[3].to_numpy(), abs=1e-12 + ) + + # Positional selection and boolean masks are both accepted. + mask = [i < 8 for i in range(len(X))] + assert pca.group_contributions(X, group=mask).to_numpy() == pytest.approx( + against_centre.to_numpy(), abs=1e-12 + ) + with pytest.raises(ValueError, match="neither index labels"): + pca.group_contributions(X, group=["not-a-label"]) def test_pca_t2_spe_contributions() -> None: @@ -2471,7 +2520,7 @@ def test_vip_formula_correctness() -> None: def test_pls_score_contributions() -> None: - """Test PLS score_contributions method.""" + """PLS contributions use the direct weights, and sum to the X-score.""" rng = np.random.default_rng(42) X = pd.DataFrame(rng.standard_normal((40, 6)), columns=[f"X{i}" for i in range(1, 7)]) Y = pd.DataFrame(X.values @ rng.standard_normal((6, 2)) + rng.standard_normal((40, 2)) * 0.3, columns=["Y1", "Y2"]) @@ -2481,31 +2530,22 @@ def test_pls_score_contributions() -> None: plsmodel = PLS(n_components=3) plsmodel.fit(X_scaled, Y_scaled) - # Contribution from first observation to model center - t_obs = plsmodel.scores_.iloc[0].values - contrib = plsmodel.score_contributions(t_obs) - assert isinstance(contrib, pd.Series) - assert len(contrib) == X_scaled.shape[1] - assert contrib.index.tolist() == X_scaled.columns.tolist() - - # Conservation: contributions = dt @ P.T - P = plsmodel.x_loadings_.values # (K, A) - dt = -t_obs # t_end(0) - t_start - expected = dt @ P.T - assert contrib.values == pytest.approx(expected, abs=1e-12) + contrib = plsmodel.score_contributions(X_scaled, component=1) + assert isinstance(contrib, pd.DataFrame) + assert contrib.shape == X_scaled.shape + assert contrib.columns.tolist() == X_scaled.columns.tolist() - # Contribution between two observations - t_obs2 = plsmodel.scores_.iloc[1].values - contrib2 = plsmodel.score_contributions(t_obs, t_obs2) - assert len(contrib2) == X_scaled.shape[1] + # For PLS the score-generating matrix is direct_weights_ (T = X R), not the + # X-loadings, so that the terms sum to the score. + for a in (1, 2, 3): + terms = plsmodel.score_contributions(X_scaled, component=a) + assert terms.sum(axis=1).to_numpy() == pytest.approx(plsmodel.scores_[a].to_numpy(), abs=1e-12) + expected = X_scaled.to_numpy() * plsmodel.direct_weights_.to_numpy()[:, a - 1] + assert terms.to_numpy() == pytest.approx(expected, abs=1e-12) - # Weighted contributions (for T²) - contrib_w = plsmodel.score_contributions(t_obs, weighted=True) - assert len(contrib_w) == X_scaled.shape[1] - - # Subset of components (1-based) - contrib_sub = plsmodel.score_contributions(t_obs, components=[1, 2]) - assert len(contrib_sub) == X_scaled.shape[1] + shift = plsmodel.group_contributions(X_scaled, group=[0, 1, 2], reference=[3, 4, 5]) + expected_shift = plsmodel.scores_[1].to_numpy()[:3].mean() - plsmodel.scores_[1].to_numpy()[3:6].mean() + assert shift.sum() == pytest.approx(expected_shift, abs=1e-12) def test_pls_detect_outliers() -> None: @@ -2538,7 +2578,7 @@ def test_pls_detect_outliers() -> None: def test_pls_score_contributions_ldpe( fixture_pls_ldpe_example: dict[str, pd.DataFrame | np.ndarray | float | int], ) -> None: - """Test PLS score_contributions on LDPE real dataset.""" + """PLS score contributions on the LDPE data: the identity holds on real data.""" data = fixture_pls_ldpe_example X = pd.DataFrame(data["X"]) Y = pd.DataFrame(data["Y"]) @@ -2548,21 +2588,15 @@ def test_pls_score_contributions_ldpe( plsmodel = PLS(n_components=3) plsmodel.fit(X_scaled, Y_scaled) - # Contribution from the last observation (process fault) to model center - t_obs = plsmodel.scores_.iloc[-1].values - contrib = plsmodel.score_contributions(t_obs) - assert isinstance(contrib, pd.Series) - assert len(contrib) == X_scaled.shape[1] - - # Conservation check: contributions = dt @ P.T - P = plsmodel.x_loadings_.values - dt = -t_obs - expected = dt @ P.T - assert contrib.values == pytest.approx(expected, abs=1e-12) - - # Weighted contributions - contrib_w = plsmodel.score_contributions(t_obs, weighted=True) - assert len(contrib_w) == X_scaled.shape[1] + contrib = plsmodel.score_contributions(X_scaled, component=1) + assert contrib.shape == X_scaled.shape + assert contrib.sum(axis=1).to_numpy() == pytest.approx(plsmodel.scores_[1].to_numpy(), abs=1e-10) + + # The known process fault sits in the last observations; the variable blamed + # there must not simply be the variable with the largest loading, or the + # diagnostic would carry no per-observation information at all. + per_observation = {contrib.iloc[i].abs().idxmax() for i in range(len(contrib))} + assert len(per_observation) > 1 def test_pls_detect_outliers_ldpe( From 64df3d26003ea142fe3b2f3ac99144000f9708c6 Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 25 Jul 2026 18:49:07 +0000 Subject: [PATCH 2/5] Document contributions as the decomposition they are, and bump to 1.61.0 The user guide already carried the correct formula, x_ik * p_ka, next to code that called the old back-projection: the prose and the implementation had drifted apart. The guide now states the property that makes the formula the right one, namely that the terms sum to the score, says plainly that a contribution is not a loading and why that distinction is the reason the diagnostic exists, and covers the group and two-period forms. Both case-study notebooks print the sum of the contributions alongside the score, so a reader can see the identity hold rather than take it on trust: -1.1107 for the flagged food-texture sample, -53.4725 for the flagged tablet. Minor rather than patch: score_contributions() changes shape and its old calling convention now raises. CITATION.cff tracks the version and the release date, per the repository convention. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01D58KmvE1L75DDE7pqfJPeL --- CHANGELOG.md | 50 +++++++++++++- CITATION.cff | 4 +- docs/api/multivariate.rst | 18 +++-- docs/quickstart.rst | 4 +- .../pca-food-texture.ipynb | 16 +++-- .../pca-spectral-data.ipynb | 10 +-- docs/user_guide/pca.rst | 66 ++++++++++++++----- docs/user_guide/pls.rst | 2 +- pyproject.toml | 2 +- .../multivariate/_diagnostics.py | 6 +- tests/test_multivariate.py | 4 +- 11 files changed, 137 insertions(+), 45 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 7e561028..a2cc1b41 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,53 @@ those changes. ## [Unreleased] +## [1.61.0] - 2026-07-25 + +### Fixed + +- `score_contributions()` on `PCA`, `PLS`, `MBPCA` and `MBPLS` computed the + contribution of each variable by back-projecting a score-space difference + through the loadings, `(t_end - t_start) @ P`. That expression never receives + the observation's data. With a single component it reduces exactly to + `-t_1 * p_1`, so it returned the loading vector rescaled by a constant, and + every observation in a data set produced the same ranking of variables. On the + food-texture data all 50 observations gave one ranking; the correct + calculation gives 34. The variable blamed for the movement differed for 36 of + those 50 observations, and for 47 of the 54 LDPE observations. + + A contribution is the term-by-term breakdown of the sum that forms the score, + so the terms have to add up to the score being decomposed: + `c_ik = x_ik * R_ka` with `sum_k c_ik = t_ia`, where `R` is the + score-generating matrix (`loadings_` for PCA, `direct_weights_` for PLS). See + Miller, Swanson and Heckler (1994), who introduced contribution plots, and + Westerhuis, Gurden and Smilde (2000) for the generalisation to any + latent-variable model. `t2_contributions()` and `spe_contributions()` were + already correct and are unchanged. + +### Added + +- `group_contributions()` on `PCA`, `PLS`, `MBPCA` and `MBPLS`, for the group + and two-period forms in the same paper: the data rows are averaged (or + differenced between two sets of rows) before being multiplied by the weights. + This answers "what do these observations have in common?" and "what changed + between these two periods?", which is what a cluster on a score plot, or a + level shift part-way through a data set, actually poses. +- `scaling=` on `score_contributions()`, offering the maximum-contribution and + within-observation presentation scalings from Miller, Swanson and Heckler + (1994). The default leaves the contributions unscaled, so they sum to the + score. + +### Changed + +- **Breaking.** `score_contributions()` now takes the preprocessed data `X` and + returns one row per observation, matching `t2_contributions()` and + `spe_contributions()` beside it, rather than taking a score vector and + returning a single Series. The old signature cannot be supported alongside the + new one, because a score vector does not carry the information the + calculation needs. Calls in the old form raise a `TypeError` naming the + replacement rather than silently returning a different number. Use + `t2_contributions(X)` in place of the former `weighted=True`. + ## [1.60.0] - 2026-07-23 ### Added @@ -2652,7 +2699,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.60.0...HEAD +[Unreleased]: https://github.com/kgdunn/process-improve/compare/v1.61.0...HEAD +[1.61.0]: https://github.com/kgdunn/process-improve/compare/v1.60.0...v1.61.0 [1.60.0]: https://github.com/kgdunn/process-improve/compare/v1.59.0...v1.60.0 [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 diff --git a/CITATION.cff b/CITATION.cff index 3b94d0e9..8990021e 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -12,8 +12,8 @@ authors: repository-code: "https://github.com/kgdunn/process-improve" url: "https://kgdunn.github.io/process-improve/" license: MIT -version: 1.60.0 -date-released: "2026-07-23" +version: 1.61.0 +date-released: "2026-07-25" keywords: - chemometrics - multivariate analysis diff --git a/docs/api/multivariate.rst b/docs/api/multivariate.rst index bf4199af..6411b8cd 100644 --- a/docs/api/multivariate.rst +++ b/docs/api/multivariate.rst @@ -10,7 +10,7 @@ PCA ~~~ .. autoclass:: PCA - :members: fit, transform, fit_transform, predict, score, select_n_components, score_contributions, detect_outliers + :members: fit, transform, fit_transform, predict, score, select_n_components, score_contributions, group_contributions, detect_outliers :undoc-members: :show-inheritance: @@ -18,7 +18,7 @@ PLS ~~~ .. autoclass:: PLS - :members: fit, transform, fit_transform, predict, score, select_n_components, score_contributions, detect_outliers, cross_validate + :members: fit, transform, fit_transform, predict, score, select_n_components, score_contributions, group_contributions, detect_outliers, cross_validate :undoc-members: :show-inheritance: @@ -94,11 +94,11 @@ also bound as a convenience method on the model after :meth:`fit`. answer different questions about the same fitted score matrix. * :meth:`PCA.score_contributions` (and :meth:`PLS.score_contributions`) is - *per-variable* and signed. It decomposes a single observation's movement - in score space back onto the original variables, answering "which - **variables** explain why this observation sits where it does?". It - returns one signed value per variable, and it takes an observation's - score vector as input. + *per-variable* and signed. It splits each score into the K terms + :math:`x_{ik} R_{ka}` that form it, answering "which **variables** + explain why this observation sits where it does?". It takes the + preprocessed data and returns a sample-by-variable table whose rows sum + to the score being decomposed. * :func:`observation_contributions` is *per-observation* and non-negative. It reports each observation's share of a component's total inertia @@ -116,6 +116,10 @@ also bound as a convenience method on the model after :meth:`fit`. .. autofunction:: observation_contributions +.. autofunction:: score_contributions + +.. autofunction:: group_contributions + .. autofunction:: eigenvalue_summary .. autofunction:: project_variables diff --git a/docs/quickstart.rst b/docs/quickstart.rst index cf87f626..61f52d1c 100644 --- a/docs/quickstart.rst +++ b/docs/quickstart.rst @@ -31,7 +31,7 @@ PCA Example # Diagnostics pca.detect_outliers() - pca.score_contributions(pca.scores_.iloc[0].values) + pca.score_contributions(X_scaled, component=1) # per-variable, sums to t1 # Plots pca.score_plot() @@ -66,7 +66,7 @@ PLS Example # Diagnostics pls.detect_outliers() - pls.score_contributions(pls.scores_.iloc[0].values) + pls.score_contributions(X_scaled, component=1) # per-variable, sums to t1 Component Selection ------------------- diff --git a/docs/user_guide/case_studies/latent-variable-modelling/pca-food-texture.ipynb b/docs/user_guide/case_studies/latent-variable-modelling/pca-food-texture.ipynb index ad0cda00..d7e866fa 100644 --- a/docs/user_guide/case_studies/latent-variable-modelling/pca-food-texture.ipynb +++ b/docs/user_guide/case_studies/latent-variable-modelling/pca-food-texture.ipynb @@ -145,7 +145,12 @@ "source": [ "## Outlier explanation\n", "\n", - "Each flagged sample can be explained back in terms of the original five variables with `score_contributions`." + "Each flagged sample can be explained back in terms of the original five\n", + "variables with `score_contributions`. A score is a weighted sum of the\n", + "variables, so it splits into one term per variable, `x_ik * p_ka`, and those\n", + "terms add back up to the score. Note that this is not the loading plot: a\n", + "variable with a large loading contributes nothing to a sample sitting at its\n", + "mean, so the ranking here is specific to this sample." ] }, { @@ -168,15 +173,18 @@ "source": [ "if outliers:\n", " worst = outliers[0][\"observation\"]\n", - " contrib = model.score_contributions(model.scores_.loc[worst].values)\n", + " contrib = model.score_contributions(X, component=1).loc[worst]\n", "\n", " fig, ax = plt.subplots(figsize=(6, 3))\n", " contrib.plot.bar(ax=ax, color=\"steelblue\")\n", " ax.axhline(0, color=\"k\", linewidth=0.5)\n", - " ax.set_ylabel(\"Contribution\")\n", + " ax.set_ylabel(\"Contribution to $t_1$\")\n", " ax.set_title(f\"Score contributions for sample {worst}\")\n", " plt.tight_layout()\n", - " plt.show()" + " plt.show()\n", + "\n", + " print(f\"contributions sum to {contrib.sum():.4f}\")\n", + " print(f\"score t1 for this sample is {model.scores_.loc[worst, 1]:.4f}\")" ] }, { diff --git a/docs/user_guide/case_studies/latent-variable-modelling/pca-spectral-data.ipynb b/docs/user_guide/case_studies/latent-variable-modelling/pca-spectral-data.ipynb index 1412505c..a80c2900 100644 --- a/docs/user_guide/case_studies/latent-variable-modelling/pca-spectral-data.ipynb +++ b/docs/user_guide/case_studies/latent-variable-modelling/pca-spectral-data.ipynb @@ -261,17 +261,19 @@ "outputs": [], "source": [ "worst = outliers[0][\"observation\"]\n", - "scores_worst = model.scores_.loc[worst].values\n", - "contrib = model.score_contributions(scores_worst)\n", + "contrib = model.score_contributions(X, component=1).loc[worst]\n", "\n", "fig, ax = plt.subplots(figsize=(8, 3))\n", "ax.plot(range(len(contrib)), contrib.values, linewidth=0.7)\n", "ax.axhline(0, color=\"k\", linewidth=0.5)\n", "ax.set_xlabel(\"Wavelength index\")\n", - "ax.set_ylabel(\"Contribution\")\n", + "ax.set_ylabel(\"Contribution to $t_1$\")\n", "ax.set_title(f\"Score contributions for tablet {worst}\")\n", "plt.tight_layout()\n", - "plt.show()" + "plt.show()\n", + "\n", + "print(f\"contributions sum to {contrib.sum():.4f}\")\n", + "print(f\"score t1 for this tablet is {model.scores_.loc[worst, 1]:.4f}\")" ] }, { diff --git a/docs/user_guide/pca.rst b/docs/user_guide/pca.rst index ddd64a16..40a47462 100644 --- a/docs/user_guide/pca.rst +++ b/docs/user_guide/pca.rst @@ -167,31 +167,61 @@ component accounts for. Score Contributions ------------------- -Score contributions decompose a score-space movement back into the original -variable space, answering: *"Which variables caused this observation to score -where it did?"* +Score contributions decompose a score back into the original variable space, +answering: *"Which variables caused this observation to score where it did?"* -Each score value :math:`t_{i,a}` can be written as a sum of K contributions -:math:`x_{i,k} \, p_{k,a}`, one per original variable. Plotting these as a -bar chart reveals the dominant drivers. +A score is a weighted sum of the variables, so it splits exactly into K terms, +one per variable: -The reference point matters: contributions always measure the difference -*from* one point *to* another. The default reference is the model center -(the origin after centering), but you can compare any two points or groups: +.. math:: + + t_{i,a} = \sum_{k=1}^{K} x_{i,k} \, p_{k,a} + \qquad\Longrightarrow\qquad + c_{i,k}^{(a)} = x_{i,k} \, p_{k,a} + +The terms sum back to the score, which is what makes a bar chart of them +readable: each bar is the amount of :math:`t_{i,a}` that variable *k* supplied. +For PLS the weights are ``direct_weights_`` rather than the loadings, since +those are what generate the scores. + +A contribution is not a loading. A loading :math:`p_{k,a}` describes the whole +data set; a contribution :math:`x_{i,k} p_{k,a}` describes one observation, and +a variable with a large loading contributes nothing to an observation sitting +at its mean. Ranking variables by loading can therefore point at a different +cause than ranking them by contribution: that difference is the reason the +diagnostic exists (Miller, Swanson and Heckler, 1994). .. code-block:: python - # Why does observation 5 differ from the model center? - contrib = model.score_contributions(model.scores_.iloc[5].values) + # Why does each observation score where it does on component 1? + contrib = model.score_contributions(X_scaled, component=1) + contrib.sum(axis=1) # equals model.scores_[1] + contrib.iloc[5].sort_values() # the drivers for observation 5 - # Why do observations 5 and 10 differ from each other? - contrib = model.score_contributions( - model.scores_.iloc[5].values, - model.scores_.iloc[10].values, - ) +Pass the same preprocessed data used to fit the model. To compare *groups* of +observations rather than one at a time, average the rows first with +``group_contributions``: this is the question a cluster on a score plot, or a +level shift part-way through a data set, actually poses. + +.. code-block:: python + + # What do these five observations have in common? + model.group_contributions(X_scaled, group=[31, 42, 47, 20, 21]) + + # What changed between the first and second halves of the campaign? + model.group_contributions(X_scaled, group=range(64, 74), reference=range(74, 84)) + +For contributions to Hotelling's :math:`T^2`, which pools every component +rather than reading one at a time, use ``model.t2_contributions(X_scaled)``. +Those sum to the observation's :math:`T^2`. + +.. note:: - # T²-weighted contributions (scale by 1/sqrt(eigenvalue)) - contrib = model.score_contributions(model.scores_.iloc[5].values, weighted=True) + Before version 1.61.0 ``score_contributions`` took a *score vector* and + back-projected it through the loadings. That expression never saw the + observation's data, so it returned the loading vector rescaled by a + constant and gave every observation the same ranking of variables. Calls in + the old form now raise a ``TypeError`` naming the replacement. Observation Contributions ------------------------- diff --git a/docs/user_guide/pls.rst b/docs/user_guide/pls.rst index d710b808..3e7f21a8 100644 --- a/docs/user_guide/pls.rst +++ b/docs/user_guide/pls.rst @@ -105,7 +105,7 @@ These diagnostics work identically to PCA (see :doc:`pca`): correlation structure. - ``model.hotellings_t2_`` and ``model.hotellings_t2_limit()`` - extremity within the model. -- ``model.score_contributions()`` - decompose scores back to original +- ``model.score_contributions(X)`` - decompose a score back to original variables. - ``model.detect_outliers()`` - combined statistical + robust ESD detection. - ``model.squared_cosine()`` - quality of representation of each observation, diff --git a/pyproject.toml b/pyproject.toml index deed236b..cdc0efac 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "process-improve" -version = "1.60.0" +version = "1.61.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/multivariate/_diagnostics.py b/src/process_improve/multivariate/_diagnostics.py index c4a60bc1..4b5f7544 100644 --- a/src/process_improve/multivariate/_diagnostics.py +++ b/src/process_improve/multivariate/_diagnostics.py @@ -671,10 +671,10 @@ def _score_contribution_terms( def _reject_old_score_contributions_api(X: object, keywords: dict) -> None: - """Raise a migration error if called with the pre-1.53.0 signature.""" + """Raise a migration error if called with the pre-1.61.0 signature.""" looks_like_scores = isinstance(X, (pd.Series, np.ndarray, list, tuple)) and np.ndim(X) == 1 if looks_like_scores or (set(keywords) & _OLD_SCORE_CONTRIBUTIONS_KEYWORDS): - raise TypeError(_OLD_SCORE_CONTRIBUTIONS_API.format(version="1.53.0")) + raise TypeError(_OLD_SCORE_CONTRIBUTIONS_API.format(version="1.61.0")) if keywords: unexpected = ", ".join(sorted(keywords)) msg = f"score_contributions() got an unexpected keyword argument: {unexpected}." @@ -734,7 +734,7 @@ def score_contributions( the *pattern* of bars within a row unchanged; neither preserves the sum to the score. **deprecated - Captures ``t_end``, ``components`` and ``weighted`` from the pre-1.53.0 + Captures ``t_end``, ``components`` and ``weighted`` from the pre-1.61.0 calling convention, which took a score vector instead of ``X``, so that those calls raise a :class:`TypeError` explaining the change rather than quietly returning something else. Passing a 1-D ``X`` raises the same diff --git a/tests/test_multivariate.py b/tests/test_multivariate.py index cd7c93f0..8f834355 100644 --- a/tests/test_multivariate.py +++ b/tests/test_multivariate.py @@ -945,7 +945,7 @@ def test_pca_score_contributions() -> None: # --- The diagnosis must depend on the observation, not only on the model. # A back-projection of the score through the loadings (what this method did - # before 1.53.0) gives every observation the same ranking of variables, + # before 1.61.0) gives every observation the same ranking of variables, # which is the loading ranking; that is the defect this test guards. rankings = {tuple(np.argsort(-contrib.iloc[i].abs().to_numpy())) for i in range(len(X))} assert len(rankings) > 1 @@ -969,7 +969,7 @@ def test_pca_score_contributions() -> None: def test_score_contributions_rejects_the_old_api() -> None: - """The pre-1.53.0 signature took a score vector and must not silently work.""" + """The pre-1.61.0 signature took a score vector and must not silently work.""" rng = np.random.default_rng(3) X = MCUVScaler().fit_transform( pd.DataFrame(rng.standard_normal((25, 4)), columns=list("abcd")) From af70daca61429f72c31f3b2e2453d8cd6c08be23 Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 25 Jul 2026 19:03:57 +0000 Subject: [PATCH 3/5] Select contribution groups by label only, and cover the new branches Codecov flagged 45 uncovered lines on the patch. Writing tests for them turned up a real defect in group_contributions() rather than only a coverage gap. _select_rows fell back to positional selection when a label was missing, which made the meaning depend on the data. On the LDPE blocks, indexed 1..54, group=[0..9] selected positions (0 is not a label) while reference=[10..19] selected labels, so the two halves of one comparison were resolved by different rules and the result was silently off by one row. Selection is now by index label or boolean mask only; callers wanting positions pass X.index[...], which the error message says. A rule that quietly changes with the index is worse than one that refuses. Coverage of the new code: _scale_block_contributions is now fully exercised, including both of the paper's scalings and the joint-over-blocks behaviour that distinguishes them from scaling each block separately. _diagnostics reaches 97% under the contribution tests alone, with the remaining misses in unrelated functions. The MBPCA and MBPLS guards were each tested on one class only; they are now tested on both, along with blocks supplied as plain arrays. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01D58KmvE1L75DDE7pqfJPeL --- CHANGELOG.md | 4 +- docs/user_guide/pca.rst | 7 +- .../multivariate/_diagnostics.py | 36 +++++--- tests/test_multiblock_reference.py | 88 ++++++++++++++++++- tests/test_multivariate.py | 48 +++++++++- 5 files changed, 162 insertions(+), 21 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index a2cc1b41..ad0fe456 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -41,7 +41,9 @@ those changes. differenced between two sets of rows) before being multiplied by the weights. This answers "what do these observations have in common?" and "what changed between these two periods?", which is what a cluster on a score plot, or a - level shift part-way through a data set, actually poses. + level shift part-way through a data set, actually poses. Observations are + selected by index label or boolean mask; pass `X.index[...]` to select by + position. - `scaling=` on `score_contributions()`, offering the maximum-contribution and within-observation presentation scalings from Miller, Swanson and Heckler (1994). The default leaves the contributions unscaled, so they sum to the diff --git a/docs/user_guide/pca.rst b/docs/user_guide/pca.rst index 40a47462..937ab3c6 100644 --- a/docs/user_guide/pca.rst +++ b/docs/user_guide/pca.rst @@ -209,7 +209,12 @@ level shift part-way through a data set, actually poses. model.group_contributions(X_scaled, group=[31, 42, 47, 20, 21]) # What changed between the first and second halves of the campaign? - model.group_contributions(X_scaled, group=range(64, 74), reference=range(74, 84)) + model.group_contributions( + X_scaled, group=X_scaled.index[64:74], reference=X_scaled.index[74:84] + ) + +Observations are selected by index label, or by a boolean mask the same length +as ``X``. To select by position, pass ``X.index[...]`` as above. For contributions to Hotelling's :math:`T^2`, which pools every component rather than reading one at a time, use ``model.t2_contributions(X_scaled)``. diff --git a/src/process_improve/multivariate/_diagnostics.py b/src/process_improve/multivariate/_diagnostics.py index 4b5f7544..d424e475 100644 --- a/src/process_improve/multivariate/_diagnostics.py +++ b/src/process_improve/multivariate/_diagnostics.py @@ -823,8 +823,9 @@ def group_contributions( X : array-like of shape (n_samples, n_features) Preprocessed data, scaled the same way as the training data. group : sequence - Index labels (or a boolean mask, or positional integers) selecting the - observations of interest. + Index labels of the observations of interest, or a boolean mask the + same length as ``X``. Selection is by label, not by position; pass + ``X.index[...]`` to select positionally. reference : sequence, optional Index labels selecting the observations to compare against. ``None`` (default) compares the group against the model centre. @@ -842,8 +843,10 @@ def group_contributions( >>> pca = PCA(n_components=2).fit(X_scaled) >>> # Five batches that cluster together on the score plot: >>> pca.group_contributions(X_scaled, group=[31, 142, 147, 220, 221]) - >>> # What shifted at batch 74? - >>> pca.group_contributions(X_scaled, group=range(64, 74), reference=range(74, 84)) + >>> # What shifted at batch 74? (ten batches either side, by position) + >>> pca.group_contributions( + ... X_scaled, group=X_scaled.index[64:74], reference=X_scaled.index[74:84] + ... ) See Also -------- @@ -866,7 +869,14 @@ def _mean_row(selector: Sequence, name: str) -> np.ndarray: def _select_rows(X_df: pd.DataFrame, selector: Sequence, name: str) -> pd.DataFrame: - """Select rows of ``X_df`` by label, boolean mask, or positional index.""" + """Select rows of ``X_df`` by index label or by boolean mask. + + Deliberately label-based only. Falling back to positional selection when a + label happens to be missing would make the meaning depend on the data: on a + frame indexed 1..54, ``[0, 1, 2]`` would be positions (0 is not a label) + while ``[10, 11, 12]`` would be labels, silently selecting a different set + of rows. Callers who want positions pass ``X.index[...]``. + """ values = list(selector) if values and all(isinstance(v, (bool, np.bool_)) for v in values): if len(values) != len(X_df): @@ -874,14 +884,14 @@ def _select_rows(X_df: pd.DataFrame, selector: Sequence, name: str) -> pd.DataFr raise ValueError(msg) return X_df.loc[np.asarray(values, dtype=bool)] missing = [v for v in values if v not in X_df.index] - if not missing: - return X_df.loc[values] - # Fall back to positional selection, but only if every entry is a valid - # position; a partial match means the caller mixed labels and positions. - if all(isinstance(v, (int, np.integer)) and 0 <= int(v) < len(X_df) for v in values): - return X_df.iloc[[int(v) for v in values]] - msg = f"{name} contains entries that are neither index labels of X nor valid row positions: {missing}." - raise ValueError(msg) + if missing: + msg = ( + f"{name} contains entries that are not index labels of X: {missing}. " + "Selection is by index label or boolean mask; to select by position, " + "pass X.index[...] instead." + ) + raise ValueError(msg) + return X_df.loc[values] def eigenvalue_summary(model: BaseEstimator) -> pd.DataFrame: diff --git a/tests/test_multiblock_reference.py b/tests/test_multiblock_reference.py index 17d689b0..82a3b835 100644 --- a/tests/test_multiblock_reference.py +++ b/tests/test_multiblock_reference.py @@ -402,10 +402,72 @@ def test_group_contributions_sum_to_the_average_super_score(self, synthetic_two_ x_blocks, _ = synthetic_two_block model = MBPCA(n_components=2).fit(x_blocks) - group = list(range(5)) - contribs = model.group_contributions(x_blocks, group=group) + index = next(iter(x_blocks.values())).index + contribs = model.group_contributions(x_blocks, group=index[:5]) total = sum(float(series.sum()) for series in contribs.values()) - assert total == pytest.approx(model.super_scores_[1].values[group].mean(), abs=1e-10) + assert total == pytest.approx(model.super_scores_[1].values[:5].mean(), abs=1e-10) + + def test_contribution_guards_and_plain_array_blocks(self, synthetic_two_block) -> None: + """Guards mirror MBPLS, and blocks may arrive as plain arrays.""" + from process_improve.multivariate.methods import MBPCA + + x_blocks, _ = synthetic_two_block + model = MBPCA(n_components=2).fit(x_blocks) + + with pytest.raises(ValueError, match="Missing X-blocks"): + model.score_contributions({model.block_names_[0]: x_blocks[model.block_names_[0]]}) + with pytest.raises(ValueError, match="1-based index"): + model.score_contributions(x_blocks, component=3) + + # A block given as a bare ndarray is labelled from the fitted columns. + as_arrays = {name: block.to_numpy() for name, block in x_blocks.items()} + from_arrays = model.score_contributions(as_arrays, component=1) + from_frames = model.score_contributions(x_blocks, component=1) + for name in model.block_names_: + np.testing.assert_array_almost_equal( + from_arrays[name].values, from_frames[name].values, decimal=12 + ) + + # The two-period form, against a reference set of observations. + index = next(iter(x_blocks.values())).index + shift = model.group_contributions(x_blocks, group=index[:5], reference=index[5:10]) + total = sum(float(series.sum()) for series in shift.values()) + scores = model.super_scores_[1].values + assert total == pytest.approx(scores[:5].mean() - scores[5:10].mean(), abs=1e-10) + + def test_score_contributions_scalings_span_the_blocks(self, synthetic_two_block) -> None: + """Both scalings are taken over the blocks jointly, not block by block. + + A super score pools every block, so scaling each block against its own + maximum would make bars from different blocks incomparable, which is the + one thing a multi-block contribution plot is read for. + """ + from process_improve.multivariate.methods import MBPCA + + x_blocks, _ = synthetic_two_block + model = MBPCA(n_components=2).fit(x_blocks) + raw = model.score_contributions(x_blocks, component=1) + + largest = model.score_contributions(x_blocks, component=1, scaling="maximum") + peak = max(float(np.abs(f.values).max()) for f in largest.values()) + assert peak == pytest.approx(1.0) + # Exactly one block holds the peak: the divisor was shared. + peaks = [float(np.abs(f.values).max()) for f in largest.values()] + assert sum(p == pytest.approx(1.0) for p in peaks) == 1 + + within = model.score_contributions(x_blocks, component=1, scaling="within") + row_total = sum(np.abs(f.values).sum(axis=1) for f in within.values()) + np.testing.assert_array_almost_equal(row_total, np.ones(len(row_total)), decimal=10) + + # Scaling changes the units, not the pattern within a row. + for name in model.block_names_: + ratio = within[name].values[0] / raw[name].values[0] + np.testing.assert_array_almost_equal(ratio, np.full_like(ratio, ratio[0]), decimal=10) + + with pytest.raises(ValueError, match="scaling must be one of"): + model.score_contributions(x_blocks, scaling="loud") + with pytest.raises(TypeError, match="dict"): + model.score_contributions(next(iter(x_blocks.values()))) def test_super_score_and_loadings_plots_return_figures(self, synthetic_two_block) -> None: import plotly.graph_objects as go @@ -1162,6 +1224,26 @@ def test_score_contributions_differ_between_components(self, ldpe) -> None: with pytest.raises(ValueError, match="Missing X-blocks"): model.score_contributions({model.block_names_[0]: x_blocks[model.block_names_[0]]}) + def test_group_contributions_sum_to_the_shift_in_super_score(self, ldpe) -> None: + """The two-period form: contributions sum to the change in average super score.""" + from process_improve.multivariate.methods import MBPLS + + x_blocks, y_df = ldpe + model = MBPLS(n_components=2).fit(x_blocks, y_df) + index = next(iter(x_blocks.values())).index + early, late = index[:10], index[10:20] + + contribs = model.group_contributions(x_blocks, group=early, reference=late) + total = sum(float(series.sum()) for series in contribs.values()) + scores = model.super_scores_[1].values + assert total == pytest.approx(scores[:10].mean() - scores[10:20].mean(), abs=1e-9) + + # Against the model centre instead of another period. + centred = model.group_contributions(x_blocks, group=early) + assert sum(float(s.sum()) for s in centred.values()) == pytest.approx( + scores[:10].mean(), abs=1e-9 + ) + def test_super_score_plot_returns_plotly_figure(self, ldpe) -> None: import plotly.graph_objects as go diff --git a/tests/test_multivariate.py b/tests/test_multivariate.py index 8f834355..30e42893 100644 --- a/tests/test_multivariate.py +++ b/tests/test_multivariate.py @@ -986,6 +986,46 @@ def test_score_contributions_rejects_the_old_api() -> None: pca.score_contributions(X, t_end=pca.scores_.iloc[1]) +def test_score_contributions_selector_and_argument_errors() -> None: + """The selector and argument guards each report what was actually wrong.""" + rng = np.random.default_rng(5) + X = MCUVScaler().fit_transform( + pd.DataFrame(rng.standard_normal((12, 3)), columns=list("abc")) + ) + pca = PCA(n_components=2).fit(X) + + # A keyword that was never part of any signature is reported as such, + # rather than being mistaken for the pre-1.61.0 calling convention. + with pytest.raises(TypeError, match="unexpected keyword argument: colour"): + pca.score_contributions(X, colour="red") + + # An empty group has no mean to take. + with pytest.raises(ValueError, match="group selected no observations"): + pca.group_contributions(X, group=[]) + with pytest.raises(ValueError, match="reference selected no observations"): + pca.group_contributions(X, group=[0], reference=[]) + + # A boolean mask has to line up with the rows it is masking. + with pytest.raises(ValueError, match="boolean mask of length 3, but X has 12 rows"): + pca.group_contributions(X, group=[True, False, True]) + + # Selection is by label. Positions go through X.index[...], and an entry + # that is not a label is refused rather than quietly read as a position: + # on a frame indexed 1..N that fallback would make [0, 1, 2] mean positions + # and [10, 11, 12] mean labels. + assert pca.group_contributions(X, group=X.index[:2]).to_numpy() == pytest.approx( + pca.score_contributions(X, component=1).iloc[:2].mean(axis=0).to_numpy(), abs=1e-12 + ) + with pytest.raises(ValueError, match="not index labels of X"): + pca.group_contributions(X, group=[999]) + + shifted = X.copy() + shifted.index = range(1, len(X) + 1) + pca_shifted = PCA(n_components=2).fit(shifted) + with pytest.raises(ValueError, match=r"not index labels of X: \[0\]"): + pca_shifted.group_contributions(shifted, group=[0, 1, 2]) + + def test_pca_group_contributions() -> None: """Group contributions sum to the average score, or to the shift between groups.""" rng = np.random.default_rng(11) @@ -1010,12 +1050,12 @@ def test_pca_group_contributions() -> None: pca.score_contributions(X, component=1).iloc[3].to_numpy(), abs=1e-12 ) - # Positional selection and boolean masks are both accepted. + # A boolean mask is accepted alongside labels. mask = [i < 8 for i in range(len(X))] assert pca.group_contributions(X, group=mask).to_numpy() == pytest.approx( against_centre.to_numpy(), abs=1e-12 ) - with pytest.raises(ValueError, match="neither index labels"): + with pytest.raises(ValueError, match="not index labels of X"): pca.group_contributions(X, group=["not-a-label"]) @@ -2543,7 +2583,9 @@ def test_pls_score_contributions() -> None: expected = X_scaled.to_numpy() * plsmodel.direct_weights_.to_numpy()[:, a - 1] assert terms.to_numpy() == pytest.approx(expected, abs=1e-12) - shift = plsmodel.group_contributions(X_scaled, group=[0, 1, 2], reference=[3, 4, 5]) + shift = plsmodel.group_contributions( + X_scaled, group=X_scaled.index[:3], reference=X_scaled.index[3:6] + ) expected_shift = plsmodel.scores_[1].to_numpy()[:3].mean() - plsmodel.scores_[1].to_numpy()[3:6].mean() assert shift.sum() == pytest.approx(expected_shift, abs=1e-12) From 5af109efa3badd5a37538d6bd5aea75d17279205 Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 25 Jul 2026 19:16:21 +0000 Subject: [PATCH 4/5] Test contribution plots against the paper that defines them tests/test_multivariate_contribution_plots.py encodes Miller, Swanson and Heckler (1994), the paper the diagnostic comes from, one test per equation, page or figure, each citing what it encodes: https://literature.learnche.org/item/78/contribution-plots-a-missing-link-in-multivariate-quality-control Equations 1 and 2 (T-squared from the covariance form and from the scores, and that they agree when no dimension is dropped); equation 3 and its per-term form, for PCA, for PLS, and for the multiblock super score; equation 4 and the squared, non-negative Q contributions; the group average of page 14 and its limiting case where averaging over every batch gives zero on centred data; the plus and minus one-tenth level-shift weights of page 18 written out exactly as the paper prints them; the orthogonal polynomial for a drift from the same page; and both scalings of pages 20 and 21, including the proportion property and the statement that scaling leaves the pattern of bar heights unchanged. The tests are written against the paper rather than against the code. The old tests asserted the output equalled dt @ P.T, which re-derived the implementation and so held whatever it computed; that is why the defect survived. To check these are not vacuous in the same way, the old back-projection was reintroduced temporarily: 9 of the 21 PCA and PLS tests failed, including every load-bearing one. Two things came out of writing them. The paper's general form is an arbitrary linear combination of the rows, not only a mean or a difference of means, and page 18 gives a case the group/reference pair cannot express: a run of batches drifting rather than stepping, weighted by the first order orthogonal polynomial. group_ contributions() now takes weights= for that, with group and reference as sugar over it. The remaining Codecov gaps were input guards on the multiblock spe_contributions, the multiblock form of equation 4, so they are covered here too. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01D58KmvE1L75DDE7pqfJPeL --- CHANGELOG.md | 6 + .../multivariate/_diagnostics.py | 79 +- tests/test_multiblock_reference.py | 31 + tests/test_multivariate_contribution_plots.py | 691 ++++++++++++++++++ 4 files changed, 784 insertions(+), 23 deletions(-) create mode 100644 tests/test_multivariate_contribution_plots.py diff --git a/CHANGELOG.md b/CHANGELOG.md index ad0fe456..e78e309a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -44,6 +44,12 @@ those changes. level shift part-way through a data set, actually poses. Observations are selected by index label or boolean mask; pass `X.index[...]` to select by position. +- `weights=` on `group_contributions()`, giving the paper's general form + directly: contributions to any linear combination of the rows, + `(sum_i w_i x_ik) R_ka`, which sum to `sum_i w_i t_ia`. `group` and + `reference` are sugar over it. The paper suggests the first-order orthogonal + polynomial as the weights when a run of batches is drifting rather than + stepping, which the group/reference pair alone cannot express. - `scaling=` on `score_contributions()`, offering the maximum-contribution and within-observation presentation scalings from Miller, Swanson and Heckler (1994). The default leaves the contributions unscaled, so they sum to the diff --git a/src/process_improve/multivariate/_diagnostics.py b/src/process_improve/multivariate/_diagnostics.py index d424e475..b6ae97c1 100644 --- a/src/process_improve/multivariate/_diagnostics.py +++ b/src/process_improve/multivariate/_diagnostics.py @@ -789,32 +789,37 @@ def score_contributions( return pd.DataFrame(contributions, index=X_df.index, columns=X_df.columns) -def group_contributions( +def group_contributions( # noqa: PLR0913 - group/reference are sugar for weights model: BaseEstimator, X: DataMatrix, - group: Sequence, + group: Sequence | None = None, reference: Sequence | None = None, component: int = 1, + weights: Sequence | None = None, ) -> pd.Series: r"""Per-variable contributions to a group's average score, or to a shift. - The group form of :func:`score_contributions`. Averaging the data rows + The group form of :func:`score_contributions`. Combining the data rows before multiplying by the score-generating weights answers "what do these observations have in common?" rather than "why is this one observation unusual?", which is the question a cluster on a score plot, or a level shift part-way through a data set, actually poses. - With ``reference=None`` the contributions are + In general any linear combination of the rows may be used + (Miller, Swanson and Heckler, 1994): .. math:: - c_k = \bar{x}_{k}^{(G)} R_{ka}, \qquad \sum_k c_k = \bar{t}_a^{(G)}, + c_k = \Bigl(\sum_i w_i x_{ik}\Bigr) R_{ka}, \qquad + \sum_k c_k = \sum_i w_i t_{ia}. - comparing the group mean against the model centre. Given a ``reference`` - group the contributions are :math:`(\bar{x}_k^{(G)} - \bar{x}_k^{(H)}) - R_{ka}`, which sum to the difference in average score between the two - groups. Both forms are the linear-combination generalisation given by - Miller, Swanson and Heckler (1994). + The common cases have their own arguments. With ``group`` alone the weights + are :math:`1/n_G` over the group, comparing its mean against the model + centre. With ``group`` and ``reference`` they are :math:`+1/n_G` and + :math:`-1/n_H`, so the contributions sum to the difference in average score, + which is the level-shift diagnostic of the paper's Figure 9. Pass ``weights`` + directly for anything else: the paper suggests the first-order orthogonal + polynomial when a run of batches is drifting rather than stepping. Parameters ---------- @@ -822,21 +827,26 @@ def group_contributions( A fitted PCA or PLS model. X : array-like of shape (n_samples, n_features) Preprocessed data, scaled the same way as the training data. - group : sequence + group : sequence, optional Index labels of the observations of interest, or a boolean mask the same length as ``X``. Selection is by label, not by position; pass - ``X.index[...]`` to select positionally. + ``X.index[...]`` to select positionally. Required unless ``weights`` + is given. reference : sequence, optional Index labels selecting the observations to compare against. ``None`` (default) compares the group against the model centre. component : int, default=1 **1-based** component index whose score is decomposed. + weights : sequence, optional + One weight per row of ``X``, giving the linear combination directly. + Mutually exclusive with ``group`` / ``reference``. Returns ------- pd.Series - Signed contributions, one per variable. Sums to the group's average - score, or to the difference in average score between the two groups. + Signed contributions, one per variable. Sums to the weighted combination + of the scores: the group's average score, the difference in average + score between the two groups, or :math:`\sum_i w_i t_{ia}`. Examples -------- @@ -847,6 +857,11 @@ def group_contributions( >>> pca.group_contributions( ... X_scaled, group=X_scaled.index[64:74], reference=X_scaled.index[74:84] ... ) + >>> # A run of batches drifting rather than stepping: weight by a + >>> # first-order orthogonal polynomial over the run. + >>> slope = np.zeros(len(X_scaled)) + >>> slope[40:60] = np.arange(20) - 9.5 + >>> pca.group_contributions(X_scaled, weights=slope, component=3) See Also -------- @@ -854,20 +869,38 @@ def group_contributions( """ X_df, r_a = _score_contribution_terms(model, X, component) - def _mean_row(selector: Sequence, name: str) -> np.ndarray: - selected = _select_rows(X_df, selector, name) - if len(selected) == 0: - msg = f"{name} selected no observations from X." + if weights is not None: + if group is not None or reference is not None: + msg = "Pass either weights, or group/reference, not both." raise ValueError(msg) - return selected.to_numpy(dtype=float).mean(axis=0) - - deviation = _mean_row(group, "group") - if reference is not None: - deviation = deviation - _mean_row(reference, "reference") + w = np.asarray(weights, dtype=float) + if w.ndim != 1 or w.size != len(X_df): + msg = f"weights must have one entry per row of X ({len(X_df)}), got shape {w.shape}." + raise ValueError(msg) + else: + if group is None: + msg = "Pass group (optionally with reference), or weights." + raise ValueError(msg) + w = np.zeros(len(X_df), dtype=float) + w += _mean_weights(X_df, group, "group") + if reference is not None: + w -= _mean_weights(X_df, reference, "reference") + deviation = w @ X_df.to_numpy(dtype=float) return pd.Series(deviation * r_a, index=X_df.columns, name="group_contributions") +def _mean_weights(X_df: pd.DataFrame, selector: Sequence, name: str) -> np.ndarray: + """Row weights that average the selected rows: ``1/n`` on each, 0 elsewhere.""" + selected = _select_rows(X_df, selector, name) + if len(selected) == 0: + msg = f"{name} selected no observations from X." + raise ValueError(msg) + w = np.zeros(len(X_df), dtype=float) + w[X_df.index.get_indexer_for(selected.index)] = 1.0 / len(selected) + return w + + def _select_rows(X_df: pd.DataFrame, selector: Sequence, name: str) -> pd.DataFrame: """Select rows of ``X_df`` by index label or by boolean mask. diff --git a/tests/test_multiblock_reference.py b/tests/test_multiblock_reference.py index 82a3b835..bf86a909 100644 --- a/tests/test_multiblock_reference.py +++ b/tests/test_multiblock_reference.py @@ -1244,6 +1244,37 @@ def test_group_contributions_sum_to_the_shift_in_super_score(self, ldpe) -> None scores[:10].mean(), abs=1e-9 ) + def test_contribution_guards_and_plain_array_blocks(self, ldpe) -> None: + """The MBPLS guards mirror MBPCA, and blocks may arrive as plain arrays.""" + from process_improve.multivariate.methods import MBPLS + + x_blocks, y_df = ldpe + model = MBPLS(n_components=2).fit(x_blocks, y_df) + + with pytest.raises(TypeError, match="dict"): + model.score_contributions(next(iter(x_blocks.values()))) + with pytest.raises(ValueError, match="Missing X-blocks"): + model.score_contributions({model.block_names_[0]: x_blocks[model.block_names_[0]]}) + with pytest.raises(ValueError, match="1-based index"): + model.score_contributions(x_blocks, component=5) + with pytest.raises(ValueError, match="scaling must be one of"): + model.score_contributions(x_blocks, scaling="loud") + + as_arrays = {name: block.to_numpy() for name, block in x_blocks.items()} + from_arrays = model.score_contributions(as_arrays, component=1) + from_frames = model.score_contributions(x_blocks, component=1) + for name in model.block_names_: + np.testing.assert_array_almost_equal( + from_arrays[name].values, from_frames[name].values, decimal=12 + ) + + # Both scalings, taken jointly over the blocks. + largest = model.score_contributions(x_blocks, component=1, scaling="maximum") + assert max(float(np.abs(f.values).max()) for f in largest.values()) == pytest.approx(1.0) + within = model.score_contributions(x_blocks, component=1, scaling="within") + row_total = sum(np.abs(f.values).sum(axis=1) for f in within.values()) + np.testing.assert_array_almost_equal(row_total, np.ones(len(row_total)), decimal=10) + def test_super_score_plot_returns_plotly_figure(self, ldpe) -> None: import plotly.graph_objects as go diff --git a/tests/test_multivariate_contribution_plots.py b/tests/test_multivariate_contribution_plots.py new file mode 100644 index 00000000..ae919f3d --- /dev/null +++ b/tests/test_multivariate_contribution_plots.py @@ -0,0 +1,691 @@ +# (c) Kevin Dunn, 2010-2026. MIT License. Based on own private work over the years. +"""Contribution plots, checked against the paper that defines them. + +Source +------ +Miller, P., Swanson, R.E. and Heckler, C.E. (1994). "Contribution plots: a +missing link in multivariate quality control." Presented at the ASQC/ASA Fall +Technical Conference; later in *Applied Mathematics and Computer Science*, +8(4), 775-792. + +https://literature.learnche.org/item/78/contribution-plots-a-missing-link-in-multivariate-quality-control + +Every test below names the equation, page or figure of that paper it encodes. +Page numbers refer to the article as paginated in the PDF at the link above +("Contribution Plots Article - Page N" in its footer). + +Why this file exists +-------------------- +Until version 1.61.0 ``score_contributions`` back-projected a score-space +difference through the loadings, ``(t_end - t_start) @ P``. That expression +never receives the observation's data: with one component it reduces exactly to +``-t_1 * p_1``, so it returned the loading vector rescaled by a constant and +gave every observation in a data set the same ranking of variables. The tests +that existed asserted only that the output equalled ``dt @ P.T``, which +re-derives the implementation and therefore held whatever it computed. + +These tests are written against the paper instead of against the code, so they +constrain the behaviour rather than describe it. The load-bearing ones are +:func:`test_eq3_contributions_sum_to_the_score` (the defining property) and +:func:`test_contributions_are_not_the_loadings`, which is the failure the paper +was written about and the regression this file exists to prevent. +""" + +from __future__ import annotations + +import numpy as np +import pandas as pd +import pytest + +from process_improve.multivariate.methods import PCA, PLS, MCUVScaler + +PAPER = "Miller, Swanson and Heckler (1994)" + + +# --------------------------------------------------------------------------- +# Fixtures +# --------------------------------------------------------------------------- + + +@pytest.fixture +def emulsion() -> pd.DataFrame: + """Correlated process data in the shape of the paper's application. + + The paper monitors a photographic emulsion process: 230 batches of 27 + correlated measurements (flows, pressures, temperatures, pH, mixer speeds), + of which only a handful of independent phenomena drive the variation + (pages 3 and 5). This is that structure at a size a test can run: 230 rows + and 27 variables generated from six latent factors plus noise. + """ + rng = np.random.default_rng(208) + n_rows, n_vars, n_factors = 230, 27, 6 + factors = rng.standard_normal((n_rows, n_factors)) + weights = rng.standard_normal((n_factors, n_vars)) + values = factors @ weights + 0.35 * rng.standard_normal((n_rows, n_vars)) + frame = pd.DataFrame(values, columns=[f"v{k:02d}" for k in range(1, n_vars + 1)]) + # Batch numbers, so that a test using labels cannot silently be reading + # positions: the paper refers to batches 101, 208 and 74 by number. + frame.index = range(1, n_rows + 1) + return frame + + +@pytest.fixture +def emulsion_model(emulsion: pd.DataFrame) -> tuple[PCA, pd.DataFrame]: + """Fit a 10-component PCA on autoscaled data, as in Figure 1 of the paper.""" + scaled = MCUVScaler().fit_transform(emulsion) + return PCA(n_components=10).fit(scaled), scaled + + +# --------------------------------------------------------------------------- +# Equations 1 and 2, page 4: T-squared, from the covariance form and the scores +# --------------------------------------------------------------------------- + + +def test_eq1_and_eq2_agree_when_every_component_is_kept(emulsion: pd.DataFrame) -> None: + """T-squared from the covariance form equals T-squared from the scores. + + Equation (1), ``T2_i = x_i S^-1 x_i'``, and equation (2), + ``T2_i = sum_a t_ia^2 / lambda_a``, are the same quantity when the PCA + keeps every dimension. The paper notes this equivalence ("It is well known + that T2 can be computed from the PCA scores") before choosing to truncate + at ``a`` dimensions. + """ + scaled = MCUVScaler().fit_transform(emulsion) + n_vars = scaled.shape[1] + model = PCA(n_components=n_vars).fit(scaled) + + values = scaled.to_numpy() + covariance = np.cov(values, rowvar=False) + from_covariance = np.einsum("ij,jk,ik->i", values, np.linalg.inv(covariance), values) + + scores = model.scores_.to_numpy() + eigenvalues = scores.var(axis=0, ddof=1) + from_scores = ((scores**2) / eigenvalues).sum(axis=1) + + assert from_scores == pytest.approx(from_covariance, rel=1e-8) + + +def test_t2_contributions_sum_to_t2(emulsion_model: tuple[PCA, pd.DataFrame]) -> None: + """Contributions to T-squared add up to T-squared. + + The T-squared decomposition is the pooled-over-components counterpart of + equation (3); ``t2_contributions`` already satisfied it before 1.61.0 and + must keep doing so. + """ + model, scaled = emulsion_model + contributions = model.t2_contributions(scaled) + assert contributions.sum(axis=1).to_numpy() == pytest.approx( + model.hotellings_t2_.iloc[:, -1].to_numpy(), rel=1e-9 + ) + + +# --------------------------------------------------------------------------- +# Equation 3, pages 6 and 7: the contribution to a score +# --------------------------------------------------------------------------- + + +def test_eq3_contributions_sum_to_the_score(emulsion_model: tuple[PCA, pd.DataFrame]) -> None: + """The k contributions to a score add up to that score. + + Equation (3) writes the score as a weighted sum of the data, + ``t_id = sum_j x_ij p_jd``, and page 7 decomposes it into "k terms + ``x_ij p_jd`` for j = 1,...,k. These k terms are the contributions to the + score ``t_id``." Adding up to the score is what makes them contributions + rather than merely a set of numbers per variable, and it is the property + the pre-1.61.0 implementation did not have. + """ + model, scaled = emulsion_model + for component in range(1, model.n_components + 1): + contributions = model.score_contributions(scaled, component=component) + assert contributions.sum(axis=1).to_numpy() == pytest.approx( + model.scores_[component].to_numpy(), abs=1e-10 + ) + + +def test_eq3_each_term_is_the_datum_times_its_weight( + emulsion_model: tuple[PCA, pd.DataFrame], +) -> None: + """Each term is literally ``x_ij * p_jd``, not just something that sums right.""" + model, scaled = emulsion_model + for component in (1, 4, 10): + contributions = model.score_contributions(scaled, component=component) + expected = scaled.to_numpy() * model.loadings_.to_numpy()[:, component - 1] + assert contributions.to_numpy() == pytest.approx(expected, abs=1e-12) + + +def test_eq3_holds_for_pls_using_the_score_generating_weights( + emulsion: pd.DataFrame, +) -> None: + """For PLS the weights are the ones that generate the scores. + + Page 21 notes that contributions apply to "other dimension reduction + techniques", naming PLS. Westerhuis, Gurden and Smilde (2000) give the + generalisation: the decomposition uses whatever matrix ``R`` satisfies + ``T = XR``, which for PLS is ``direct_weights_`` and not the X-loadings. + Using the loadings here would break the sum in equation (3). + """ + scaled = MCUVScaler().fit_transform(emulsion) + y = pd.DataFrame( + {"quality": scaled.to_numpy() @ np.linspace(1.0, -1.0, scaled.shape[1])} + ) + model = PLS(n_components=3).fit(scaled, MCUVScaler().fit_transform(y)) + + for component in (1, 2, 3): + contributions = model.score_contributions(scaled, component=component) + assert contributions.sum(axis=1).to_numpy() == pytest.approx( + model.scores_[component].to_numpy(), abs=1e-10 + ) + expected = scaled.to_numpy() * model.direct_weights_.to_numpy()[:, component - 1] + assert contributions.to_numpy() == pytest.approx(expected, abs=1e-12) + + +# --------------------------------------------------------------------------- +# Pages 7 and 8, Figure 2: a contribution is not a loading +# --------------------------------------------------------------------------- + + +def test_contributions_are_not_the_loadings() -> None: + """A large loading contributes nothing when the observation sits at its mean. + + This is the point of the paper. Page 7: "A practical difference between + contributions and loadings occurs when some of the process variables have a + value close to zero, even though those same variables may have large + loadings." Page 8 gives the case: for batch 208 the loadings pick out + variables 3, 5, 19 and 25, while the contributions pick out 12, 13, 16 and + 17; the pressures were the real fault, confirmed by univariate charts and + traced to a replaced valve. "Thus interpreting the loadings would + potentially detect a different process problem for this batch than actually + occurred." + + Constructed to make that unambiguous: ``big_loading`` dominates the + component, but the observation under test sits exactly at its mean, so it + can have contributed nothing to that observation's score. + """ + rng = np.random.default_rng(0) + n_rows = 120 + driver = rng.standard_normal(n_rows) + correlated = ["press_a", "press_b", "press_c"] + frame = pd.DataFrame( + { + # Three variables moving together: after autoscaling, a variable + # earns a large loading on component 1 by correlating with others. + "press_a": driver + 0.05 * rng.standard_normal(n_rows), + "press_b": driver + 0.05 * rng.standard_normal(n_rows), + "press_c": driver + 0.05 * rng.standard_normal(n_rows), + # Nearly independent of them, so its loading is small. + "silver": rng.standard_normal(n_rows), + } + ) + scaled = MCUVScaler().fit_transform(frame) + + # The batch under investigation: every large-loading variable sits exactly + # at its mean, and the small-loading one is the only thing that moved. + culprit = 7 + scaled.loc[culprit, correlated] = 0.0 + scaled.loc[culprit, "silver"] = 6.0 + + model = PCA(n_components=2).fit(scaled) + loadings = model.loadings_[1].abs() + contributions = model.score_contributions(scaled, component=1).loc[culprit] + + # Reading the loadings points at the pressures ... + assert loadings.idxmax() in correlated + assert loadings["silver"] < loadings[correlated].min() + # ... but they contributed nothing to this batch, and the contributions + # point at the variable that actually moved. + assert contributions[correlated].to_numpy() == pytest.approx(np.zeros(3), abs=1e-12) + assert contributions.abs().idxmax() == "silver" + assert loadings.idxmax() != contributions.abs().idxmax() + + +def test_contributions_distinguish_observations( + emulsion_model: tuple[PCA, pd.DataFrame], +) -> None: + """Different observations get different answers. + + The regression guard for the pre-1.61.0 defect. Back-projecting the score + through the loadings gives a result proportional to the loading vector, so + the ranking of variables is the same for every observation and the plot + carries no per-observation information at all. Contributions are per-batch + by construction (page 7: "Contributions represent the particular process + variables that were unusual *for a given batch*"). + """ + model, scaled = emulsion_model + contributions = model.score_contributions(scaled, component=1) + + rankings = { + tuple(np.argsort(-contributions.iloc[i].abs().to_numpy())) for i in range(len(scaled)) + } + assert len(rankings) > len(scaled) // 2 + + blamed = {contributions.iloc[i].abs().idxmax() for i in range(len(scaled))} + assert len(blamed) > 1 + + # Not proportional to the loadings: that is the defect's signature. + loadings = model.loadings_.to_numpy()[:, 0] + for i in range(0, len(scaled), 40): + row = contributions.iloc[i].to_numpy() + if np.allclose(row, 0.0): + continue + ratios = row / np.where(np.abs(loadings) > 1e-12, loadings, np.nan) + assert not np.allclose(ratios, ratios[0], equal_nan=True) + + +# --------------------------------------------------------------------------- +# Equation 4, pages 10 and 12: contributions to Q (the residual, SPE) +# --------------------------------------------------------------------------- + + +def test_eq4_q_contributions_are_squared_residuals_that_sum_to_q( + emulsion_model: tuple[PCA, pd.DataFrame], +) -> None: + """Q splits into the squares of the k residual elements. + + Equation (4) defines ``Q = (x_i - x_hat_i)(x_i - x_hat_i)'`` and page 12 + says: "There are k elements in ``(x_i - x_hat_i)``, and the squares of these + k values are plotted as bars in a contribution plot for Q." Squares, so + unlike the score contributions these are non-negative, as in Figure 4 where + every bar points upwards. + """ + model, scaled = emulsion_model + residuals = model.spe_contributions(scaled) + + fitted = model.scores_.to_numpy() @ model.loadings_.to_numpy().T + assert residuals.to_numpy() == pytest.approx(scaled.to_numpy() - fitted, abs=1e-10) + + q_contributions = residuals**2 + assert (q_contributions.to_numpy() >= 0).all() + q_statistic = (scaled.to_numpy() - fitted) ** 2 + assert q_contributions.sum(axis=1).to_numpy() == pytest.approx( + q_statistic.sum(axis=1), rel=1e-10 + ) + assert q_contributions.sum(axis=1).to_numpy() == pytest.approx( + model.spe_.iloc[:, -1].to_numpy() ** 2, rel=1e-8 + ) + + +# --------------------------------------------------------------------------- +# Page 14: contributions for a group of batches +# --------------------------------------------------------------------------- + + +def test_page14_group_contributions_use_the_group_average( + emulsion_model: tuple[PCA, pd.DataFrame], +) -> None: + """A cluster of batches is summarised by averaging the data, then weighting. + + Page 14: "When we compute the contributions, we can replace ``x_ij`` with a + suitably chosen average, or other linear combination of the data. In this + case, we would use the average value ``x_bar_.j``, where the averaging is + over the batches of interest." The paper's example is a cluster of five + non-sequential batches (31, 142, 147, 220, 221) seen on a score plot. + """ + model, scaled = emulsion_model + cluster = [31, 142, 147, 220, 221] + + contributions = model.group_contributions(scaled, group=cluster, component=1) + + mean_row = scaled.loc[cluster].to_numpy().mean(axis=0) + expected = mean_row * model.loadings_.to_numpy()[:, 0] + assert contributions.to_numpy() == pytest.approx(expected, abs=1e-12) + assert contributions.sum() == pytest.approx(model.scores_.loc[cluster, 1].mean(), abs=1e-10) + + +def test_page14_the_comparison_is_against_the_centre_which_is_zero( + emulsion_model: tuple[PCA, pd.DataFrame], +) -> None: + """Averaging over every batch gives nothing, because the data are centred. + + Page 14: "This in effect compares the average value of the batches of + interest to the mean of all the process variables, which is zero for each + (mean centered) process variable." Taking the whole data set as the group + is that statement's limiting case. + """ + model, scaled = emulsion_model + everything = model.group_contributions(scaled, group=list(scaled.index), component=1) + assert everything.to_numpy() == pytest.approx(np.zeros(scaled.shape[1]), abs=1e-10) + + +def test_single_observation_group_is_the_single_observation_contribution( + emulsion_model: tuple[PCA, pd.DataFrame], +) -> None: + """A group of one reduces to equation (3) for that batch.""" + model, scaled = emulsion_model + one = model.group_contributions(scaled, group=[208], component=2) + direct = model.score_contributions(scaled, component=2).loc[208] + assert one.to_numpy() == pytest.approx(direct.to_numpy(), abs=1e-12) + + +# --------------------------------------------------------------------------- +# Page 18, Figures 8 and 9: a level shift between two periods +# --------------------------------------------------------------------------- + + +def test_page18_level_shift_uses_plus_and_minus_one_tenth_weights( + emulsion_model: tuple[PCA, pd.DataFrame], +) -> None: + """The shift diagnostic: ten batches before against ten batches after. + + Page 18 works the case where score 3 steps at batch 74: "the linear + combination of the data has weights of +0.1 for batches 64 through 73, -0.1 + for batches 74 through 83 and 0 elsewhere. The contributions are the k + values of ``(sum_{i=64..73} x_ij - sum_{i=74..83} x_ij) p_j3 / 10``." + + Both spellings must agree: the explicit weight vector, and the + group/reference form that is sugar for it. + """ + model, scaled = emulsion_model + before, after = list(range(64, 74)), list(range(74, 84)) + + weights = np.zeros(len(scaled)) + positions = scaled.index.get_indexer_for(before) + weights[positions] = 0.1 + weights[scaled.index.get_indexer_for(after)] = -0.1 + + by_weights = model.group_contributions(scaled, weights=weights, component=3) + by_groups = model.group_contributions(scaled, group=before, reference=after, component=3) + assert by_weights.to_numpy() == pytest.approx(by_groups.to_numpy(), abs=1e-12) + + # Written out exactly as the paper prints it. + literal = ( + scaled.loc[before].to_numpy().sum(axis=0) - scaled.loc[after].to_numpy().sum(axis=0) + ) * model.loadings_.to_numpy()[:, 2] / 10.0 + assert by_weights.to_numpy() == pytest.approx(literal, abs=1e-12) + + scores = model.scores_[3] + assert by_weights.sum() == pytest.approx( + scores.loc[before].mean() - scores.loc[after].mean(), abs=1e-10 + ) + + +def test_page18_a_drift_can_be_weighted_by_an_orthogonal_polynomial( + emulsion_model: tuple[PCA, pd.DataFrame], +) -> None: + """Any linear combination is allowed, including a slope. + + Page 18: "Sometimes, a drift upwards or downwards will be seen in the time + sequence plot of the scores. In that situation, a linear combination of the + data that estimates a slope will be useful for determining which process + variables were drifting... we could use the first order orthogonal + polynomial for n data points as the weights." + """ + model, scaled = emulsion_model + run = list(range(100, 120)) + weights = np.zeros(len(scaled)) + # First-order orthogonal polynomial over n points: centred, evenly spaced. + weights[scaled.index.get_indexer_for(run)] = np.arange(len(run)) - (len(run) - 1) / 2.0 + assert weights.sum() == pytest.approx(0.0) + + contributions = model.group_contributions(scaled, weights=weights, component=2) + + expected = (weights @ scaled.to_numpy()) * model.loadings_.to_numpy()[:, 1] + assert contributions.to_numpy() == pytest.approx(expected, abs=1e-12) + assert contributions.sum() == pytest.approx( + float(weights @ model.scores_[2].to_numpy()), abs=1e-9 + ) + + +def test_weights_and_groups_are_mutually_exclusive( + emulsion_model: tuple[PCA, pd.DataFrame], +) -> None: + """The two spellings of the same idea cannot be combined.""" + model, scaled = emulsion_model + with pytest.raises(ValueError, match="not both"): + model.group_contributions(scaled, group=[1, 2], weights=np.zeros(len(scaled))) + with pytest.raises(ValueError, match="or weights"): + model.group_contributions(scaled) + with pytest.raises(ValueError, match="one entry per row"): + model.group_contributions(scaled, weights=np.zeros(3)) + + +# --------------------------------------------------------------------------- +# Pages 20 and 21: the two scalings +# --------------------------------------------------------------------------- + + +def test_page20_maximum_contribution_scaling( + emulsion_model: tuple[PCA, pd.DataFrame], +) -> None: + """Method 1: divide by the largest absolute contribution in the data set. + + Page 20: "for dimension d, we plot ``x_ij p_jd / max_ij |x_ij p_jd|``... we + compare the contributions for batch i to the maximum, in absolute value, of + the contributions for all of the batches. If the contribution for batch i is + +/-1, then this represents the worst deviation from the mean of all of the + batches over all variables." + """ + model, scaled = emulsion_model + raw = model.score_contributions(scaled, component=1) + scaled_bars = model.score_contributions(scaled, component=1, scaling="maximum") + + assert np.abs(scaled_bars.to_numpy()).max() == pytest.approx(1.0) + assert scaled_bars.to_numpy() == pytest.approx( + raw.to_numpy() / np.abs(raw.to_numpy()).max(), abs=1e-12 + ) + # Exactly one variable-batch pair attains it: the worst deviation. + assert int((np.abs(scaled_bars.to_numpy()) > 1.0 - 1e-12).sum()) == 1 + + +def test_page21_within_batch_scaling(emulsion_model: tuple[PCA, pd.DataFrame]) -> None: + """Method 2: divide by the total absolute contribution for that batch. + + Page 21: "for dimension d, we plot ``x_ij p_jd / sum_j |x_ij p_jd|``... The + biggest bars in this method are truly the ones which contribute most to the + score for this particular batch and the height of the bar is roughly the + proportion of the variable's contribution." + """ + model, scaled = emulsion_model + scaled_bars = model.score_contributions(scaled, component=1, scaling="within") + totals = np.abs(scaled_bars.to_numpy()).sum(axis=1) + assert totals == pytest.approx(np.ones(len(scaled)), abs=1e-10) + + +def test_page21_within_batch_scaling_is_the_exact_proportion_when_signs_agree() -> None: + """Method 2 gives the exact proportion when every contribution shares a sign. + + Page 21 qualifies method 2 with that parenthesis. Constructed so the + qualification is met: with every ``x_ij p_j1`` of one sign, each scaled bar + is exactly that variable's share of the score. + """ + rng = np.random.default_rng(3) + frame = pd.DataFrame( + rng.uniform(1.0, 2.0, size=(40, 4)) + np.arange(4), columns=list("abcd") + ) + scaled = MCUVScaler().fit_transform(frame) + model = PCA(n_components=2).fit(scaled) + + raw = model.score_contributions(scaled, component=1) + within = model.score_contributions(scaled, component=1, scaling="within") + + same_sign = raw.index[ + [bool(np.all(np.sign(row) == np.sign(row[0])) and row[0] != 0) for row in raw.to_numpy()] + ] + assert len(same_sign) > 0, "fixture no longer exercises the same-sign case" + + for label in same_sign: + row = raw.loc[label].to_numpy() + share = row / row.sum() + # "The height of the bar is roughly the proportion of the variable's + # contribution": heights are magnitudes. The scaled value carries the + # sign of the contribution, so it equals the proportion up to the sign + # of the score, which is what the division by a sum of absolute values + # does when the common sign is negative. + assert np.abs(within.loc[label].to_numpy()) == pytest.approx(np.abs(share), abs=1e-10) + assert np.sign(within.loc[label].to_numpy()) == pytest.approx(np.sign(row)) + + +def test_page20_scalings_leave_the_pattern_of_bar_heights_unchanged( + emulsion_model: tuple[PCA, pd.DataFrame], +) -> None: + """Both scalings zoom, they do not reshape. + + Page 20: "The two methods either 'zoom in' or 'zoom out' on the plot, but + leave the pattern of bar heights unchanged." So within a batch every bar is + multiplied by the same positive number: signs are preserved and the ranking + of variables is untouched. + """ + model, scaled = emulsion_model + raw = model.score_contributions(scaled, component=1) + + for scaling in ("maximum", "within"): + bars = model.score_contributions(scaled, component=1, scaling=scaling) + for i in (0, 73, 207): + row_raw = raw.iloc[i].to_numpy() + row_scaled = bars.iloc[i].to_numpy() + ratio = row_scaled / row_raw + assert ratio == pytest.approx(np.full(len(ratio), ratio[0]), rel=1e-10) + assert ratio[0] > 0 + assert np.array_equal(np.argsort(-np.abs(row_scaled)), np.argsort(-np.abs(row_raw))) + + +def test_default_scaling_preserves_the_sum_to_the_score( + emulsion_model: tuple[PCA, pd.DataFrame], +) -> None: + """The scalings are for presentation; unscaled is the quantity in equation (3). + + Page 20 introduces both methods as ways "to make the variables with the + biggest contributions stand out visually". Neither preserves the sum, so + neither is the default here. + """ + model, scaled = emulsion_model + unscaled = model.score_contributions(scaled, component=1) + assert unscaled.sum(axis=1).to_numpy() == pytest.approx( + model.scores_[1].to_numpy(), abs=1e-10 + ) + for scaling in ("maximum", "within"): + bars = model.score_contributions(scaled, component=1, scaling=scaling) + assert not np.allclose(bars.sum(axis=1).to_numpy(), model.scores_[1].to_numpy()) + + +# --------------------------------------------------------------------------- +# The pre-1.61.0 API must not come back +# --------------------------------------------------------------------------- + + +def test_the_back_projection_api_is_gone(emulsion_model: tuple[PCA, pd.DataFrame]) -> None: + """Calling the old way raises rather than returning a different number. + + A score vector cannot carry the information equation (3) needs, so the old + signature cannot be supported alongside the correct one. + """ + model, scaled = emulsion_model + with pytest.raises(TypeError, match="not a score vector"): + model.score_contributions(model.scores_.iloc[0]) + with pytest.raises(TypeError, match="t2_contributions"): + model.score_contributions(scaled, weighted=True) + + +def test_contributions_are_not_proportional_to_the_loadings_on_real_data() -> None: + """The same guard, on the LDPE data set rather than a generated one. + + Real process data, so the test cannot be satisfied by a quirk of the + generator. Under the pre-1.61.0 calculation every one of these observations + produced the identical ranking of variables. + """ + import pathlib + + folder = ( + pathlib.Path(__file__).parents[1] + / "src" + / "process_improve" + / "datasets" + / "multivariate" + / "LDPE" + ) + values = pd.read_csv(folder / "LDPE.csv", index_col=0).select_dtypes("number") + scaled = MCUVScaler().fit_transform(values) + model = PCA(n_components=3).fit(scaled) + + contributions = model.score_contributions(scaled, component=1) + assert contributions.sum(axis=1).to_numpy() == pytest.approx( + model.scores_[1].to_numpy(), abs=1e-9 + ) + + blamed = {contributions.iloc[i].abs().idxmax() for i in range(len(scaled))} + assert len(blamed) > 1, f"{PAPER}: contributions must vary between observations" + + +# --------------------------------------------------------------------------- +# Page 21: "contributions can be used with other dimension reduction +# techniques" - here the multi-block models +# --------------------------------------------------------------------------- + + +@pytest.fixture +def two_blocks(emulsion: pd.DataFrame) -> dict[str, pd.DataFrame]: + """Split the emulsion variables into two blocks, as a plant would group them.""" + return {"phase1": emulsion.iloc[:, :14], "phase2": emulsion.iloc[:, 14:]} + + +def test_eq3_holds_for_the_super_score_of_a_multiblock_model( + two_blocks: dict[str, pd.DataFrame], +) -> None: + """Summed over every block and variable, the contributions give the super score. + + Equation (3) applied to a super score. The wrinkle is that the super score + at component ``a`` is formed from the block data deflated through the + earlier components, so the decomposition has to start from that same + deflated data rather than from the raw preprocessed blocks. + """ + from process_improve.multivariate.methods import MBPCA + + model = MBPCA(n_components=3).fit(two_blocks) + for component in (1, 2, 3): + contributions = model.score_contributions(two_blocks, component=component) + total = sum(frame.sum(axis=1) for frame in contributions.values()) + assert total.to_numpy() == pytest.approx( + model.super_scores_[component].to_numpy(), abs=1e-9 + ) + + +def test_page14_group_contributions_for_a_multiblock_model( + two_blocks: dict[str, pd.DataFrame], +) -> None: + """The group and two-period forms carry over to the super score.""" + from process_improve.multivariate.methods import MBPCA + + model = MBPCA(n_components=2).fit(two_blocks) + index = two_blocks["phase1"].index + before, after = index[63:73], index[73:83] + + shift = model.group_contributions(two_blocks, group=before, reference=after) + total = sum(float(series.sum()) for series in shift.values()) + scores = model.super_scores_[1] + assert total == pytest.approx( + scores.loc[before].mean() - scores.loc[after].mean(), abs=1e-9 + ) + + +def test_eq4_q_contributions_for_a_multiblock_model( + two_blocks: dict[str, pd.DataFrame], +) -> None: + """Equation (4) per block: squared residuals summing to the block SPE. + + The multi-block form of the Q contribution plot of page 12 and Figure 4. + Squares, so every bar is non-negative, and they sum across a block's + variables to that block's squared residual distance. + """ + from process_improve.multivariate.methods import MBPLS + + y = pd.DataFrame( + {"quality": two_blocks["phase1"].to_numpy() @ np.linspace(1.0, -1.0, 14)}, + index=two_blocks["phase1"].index, + ) + model = MBPLS(n_components=2).fit(two_blocks, y) + + contributions = model.spe_contributions(two_blocks) + for name, frame in contributions.items(): + assert (frame.to_numpy() >= 0).all() + assert frame.sum(axis=1).to_numpy() == pytest.approx( + model.block_spe_[name].iloc[:, -1].to_numpy() ** 2, rel=1e-8 + ) + + # The same input guards as the score-contribution side. + with pytest.raises(TypeError, match="dict"): + model.spe_contributions(two_blocks["phase1"]) + with pytest.raises(ValueError, match="Missing X-blocks"): + model.spe_contributions({"phase1": two_blocks["phase1"]}) + as_arrays = {name: block.to_numpy() for name, block in two_blocks.items()} + from_arrays = model.spe_contributions(as_arrays) + for name, frame in from_arrays.items(): + assert frame.to_numpy() == pytest.approx(contributions[name].to_numpy(), abs=1e-12) From 389b26ce6a17a039cca0015bc2a00a4345016052 Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 25 Jul 2026 19:35:26 +0000 Subject: [PATCH 5/5] Drop the version-history notes from the code, docs and tests The user guide, the test module and several docstrings carried a note describing what score_contributions() did before this release. That history belongs in the changelog, not in prose a reader meets while trying to use the function: it puts a defunct calculation in front of someone who never called it, and it dates the moment the version is no longer current. The error message for a score-vector call now states the requirement and the alternatives without narrating what changed. It is shorter and reads as usage guidance rather than a migration note. The internal helper is named for what it checks, _reject_score_vector_call, rather than for the API it used to reject. Behaviour is unchanged: the same calls raise, with the same three replacements named, and the tests that match on the message text still pass. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01D58KmvE1L75DDE7pqfJPeL --- docs/user_guide/pca.rst | 8 --- .../multivariate/_diagnostics.py | 48 ++++++++--------- tests/test_multivariate.py | 13 +++-- tests/test_multivariate_contribution_plots.py | 52 ++++++++----------- 4 files changed, 51 insertions(+), 70 deletions(-) diff --git a/docs/user_guide/pca.rst b/docs/user_guide/pca.rst index 937ab3c6..bd5472b0 100644 --- a/docs/user_guide/pca.rst +++ b/docs/user_guide/pca.rst @@ -220,14 +220,6 @@ For contributions to Hotelling's :math:`T^2`, which pools every component rather than reading one at a time, use ``model.t2_contributions(X_scaled)``. Those sum to the observation's :math:`T^2`. -.. note:: - - Before version 1.61.0 ``score_contributions`` took a *score vector* and - back-projected it through the loadings. That expression never saw the - observation's data, so it returned the loading vector rescaled by a - constant and gave every observation the same ranking of variables. Calls in - the old form now raise a ``TypeError`` naming the replacement. - Observation Contributions ------------------------- diff --git a/src/process_improve/multivariate/_diagnostics.py b/src/process_improve/multivariate/_diagnostics.py index b6ae97c1..7e36fd3e 100644 --- a/src/process_improve/multivariate/_diagnostics.py +++ b/src/process_improve/multivariate/_diagnostics.py @@ -646,35 +646,34 @@ def _score_contribution_terms( return X_df, R[:, int(component) - 1] -_OLD_SCORE_CONTRIBUTIONS_API = """\ -score_contributions() now takes the preprocessed data X, not a score vector. +_SCORE_CONTRIBUTIONS_USAGE = """\ +score_contributions() takes the preprocessed data X, not a score vector. -Until version {version} it back-projected a score-space difference through the -loadings, which does not depend on the observation's data at all: every -observation produced the same ranking of variables, namely the ranking of the -loadings. A contribution has to be x_ik * R_ka so that the per-variable terms -sum to the score. See Miller, Swanson and Heckler (1994). +A contribution is x_ik * R_ka, so the per-variable terms sum to the score being +decomposed. That needs the observation's data, which a score vector does not +carry. See Miller, Swanson and Heckler (1994). - old: model.score_contributions(model.scores_.iloc[i]) - new: model.score_contributions(X) # all observations, component 1 - model.score_contributions(X).iloc[i] # just observation i + model.score_contributions(X) # all observations, component 1 + model.score_contributions(X).iloc[i] # just observation i - old: model.score_contributions(t_a, t_b) # two observations - new: model.group_contributions(X, group=[a], reference=[b]) +To compare two observations or two groups: - old: model.score_contributions(t, weighted=True) - new: model.t2_contributions(X) # already correct; unchanged + model.group_contributions(X, group=[a], reference=[b]) + +For contributions to Hotelling's T-squared, which pools every component: + + model.t2_contributions(X) """ -_OLD_SCORE_CONTRIBUTIONS_KEYWORDS = frozenset({"t_end", "components", "weighted"}) +_SCORE_VECTOR_KEYWORDS = frozenset({"t_end", "components", "weighted"}) -def _reject_old_score_contributions_api(X: object, keywords: dict) -> None: - """Raise a migration error if called with the pre-1.61.0 signature.""" +def _reject_score_vector_call(X: object, keywords: dict) -> None: + """Raise a usage error if called with a score vector instead of ``X``.""" looks_like_scores = isinstance(X, (pd.Series, np.ndarray, list, tuple)) and np.ndim(X) == 1 - if looks_like_scores or (set(keywords) & _OLD_SCORE_CONTRIBUTIONS_KEYWORDS): - raise TypeError(_OLD_SCORE_CONTRIBUTIONS_API.format(version="1.61.0")) + if looks_like_scores or (set(keywords) & _SCORE_VECTOR_KEYWORDS): + raise TypeError(_SCORE_CONTRIBUTIONS_USAGE) if keywords: unexpected = ", ".join(sorted(keywords)) msg = f"score_contributions() got an unexpected keyword argument: {unexpected}." @@ -734,11 +733,10 @@ def score_contributions( the *pattern* of bars within a row unchanged; neither preserves the sum to the score. **deprecated - Captures ``t_end``, ``components`` and ``weighted`` from the pre-1.61.0 - calling convention, which took a score vector instead of ``X``, so that - those calls raise a :class:`TypeError` explaining the change rather than - quietly returning something else. Passing a 1-D ``X`` raises the same - error. + Rejected. Captures ``t_end``, ``components`` and ``weighted`` so that a + call passing a score vector rather than ``X`` raises a + :class:`TypeError` explaining the correct usage. Passing a 1-D ``X`` + raises the same error. Returns ------- @@ -772,7 +770,7 @@ def score_contributions( components rather than reading one at a time. spe_contributions : The residual-space counterpart. """ - _reject_old_score_contributions_api(X, deprecated) + _reject_score_vector_call(X, deprecated) X_df, r_a = _score_contribution_terms(model, X, component) contributions = X_df.to_numpy(dtype=float) * r_a diff --git a/tests/test_multivariate.py b/tests/test_multivariate.py index 30e42893..1296f7ce 100644 --- a/tests/test_multivariate.py +++ b/tests/test_multivariate.py @@ -944,9 +944,8 @@ def test_pca_score_contributions() -> None: assert terms.to_numpy() == pytest.approx(expected, abs=1e-12) # --- The diagnosis must depend on the observation, not only on the model. - # A back-projection of the score through the loadings (what this method did - # before 1.61.0) gives every observation the same ranking of variables, - # which is the loading ranking; that is the defect this test guards. + # Anything proportional to the loading vector ranks the variables the same + # way for every observation, carrying no per-observation information. rankings = {tuple(np.argsort(-contrib.iloc[i].abs().to_numpy())) for i in range(len(X))} assert len(rankings) > 1 @@ -968,8 +967,8 @@ def test_pca_score_contributions() -> None: pca.score_contributions(X, component=4) -def test_score_contributions_rejects_the_old_api() -> None: - """The pre-1.61.0 signature took a score vector and must not silently work.""" +def test_score_contributions_rejects_a_score_vector() -> None: + """A score vector cannot carry the data the calculation needs.""" rng = np.random.default_rng(3) X = MCUVScaler().fit_transform( pd.DataFrame(rng.standard_normal((25, 4)), columns=list("abcd")) @@ -994,8 +993,8 @@ def test_score_contributions_selector_and_argument_errors() -> None: ) pca = PCA(n_components=2).fit(X) - # A keyword that was never part of any signature is reported as such, - # rather than being mistaken for the pre-1.61.0 calling convention. + # An unknown keyword is reported as such, rather than being mistaken for a + # call that passed a score vector. with pytest.raises(TypeError, match="unexpected keyword argument: colour"): pca.score_contributions(X, colour="red") diff --git a/tests/test_multivariate_contribution_plots.py b/tests/test_multivariate_contribution_plots.py index ae919f3d..52cb2456 100644 --- a/tests/test_multivariate_contribution_plots.py +++ b/tests/test_multivariate_contribution_plots.py @@ -16,19 +16,15 @@ Why this file exists -------------------- -Until version 1.61.0 ``score_contributions`` back-projected a score-space -difference through the loadings, ``(t_end - t_start) @ P``. That expression -never receives the observation's data: with one component it reduces exactly to -``-t_1 * p_1``, so it returned the loading vector rescaled by a constant and -gave every observation in a data set the same ranking of variables. The tests -that existed asserted only that the output equalled ``dt @ P.T``, which -re-derives the implementation and therefore held whatever it computed. - -These tests are written against the paper instead of against the code, so they -constrain the behaviour rather than describe it. The load-bearing ones are -:func:`test_eq3_contributions_sum_to_the_score` (the defining property) and -:func:`test_contributions_are_not_the_loadings`, which is the failure the paper -was written about and the regression this file exists to prevent. +These tests are written against the paper rather than against the code, so they +constrain the behaviour rather than describe it. A test that re-derives the +implementation holds whatever that implementation computes; a test that encodes +equation (3) holds only if the code is right. + +The load-bearing ones are :func:`test_eq3_contributions_sum_to_the_score`, the +property that makes a set of per-variable numbers a decomposition of the score, +and :func:`test_contributions_are_not_the_loadings`, the confusion the paper was +written to resolve. """ from __future__ import annotations @@ -109,8 +105,7 @@ def test_t2_contributions_sum_to_t2(emulsion_model: tuple[PCA, pd.DataFrame]) -> """Contributions to T-squared add up to T-squared. The T-squared decomposition is the pooled-over-components counterpart of - equation (3); ``t2_contributions`` already satisfied it before 1.61.0 and - must keep doing so. + equation (3). """ model, scaled = emulsion_model contributions = model.t2_contributions(scaled) @@ -131,8 +126,7 @@ def test_eq3_contributions_sum_to_the_score(emulsion_model: tuple[PCA, pd.DataFr ``t_id = sum_j x_ij p_jd``, and page 7 decomposes it into "k terms ``x_ij p_jd`` for j = 1,...,k. These k terms are the contributions to the score ``t_id``." Adding up to the score is what makes them contributions - rather than merely a set of numbers per variable, and it is the property - the pre-1.61.0 implementation did not have. + rather than merely a set of numbers per variable. """ model, scaled = emulsion_model for component in range(1, model.n_components + 1): @@ -243,12 +237,11 @@ def test_contributions_distinguish_observations( ) -> None: """Different observations get different answers. - The regression guard for the pre-1.61.0 defect. Back-projecting the score - through the loadings gives a result proportional to the loading vector, so - the ranking of variables is the same for every observation and the plot - carries no per-observation information at all. Contributions are per-batch - by construction (page 7: "Contributions represent the particular process - variables that were unusual *for a given batch*"). + Anything proportional to the loading vector ranks the variables identically + for every observation, so the plot would carry no per-observation + information at all. Contributions are per-batch by construction (page 7: + "Contributions represent the particular process variables that were unusual + *for a given batch*"). """ model, scaled = emulsion_model contributions = model.score_contributions(scaled, component=1) @@ -558,15 +551,15 @@ def test_default_scaling_preserves_the_sum_to_the_score( # --------------------------------------------------------------------------- -# The pre-1.61.0 API must not come back +# Usage: the data are required, a score vector is not enough # --------------------------------------------------------------------------- -def test_the_back_projection_api_is_gone(emulsion_model: tuple[PCA, pd.DataFrame]) -> None: - """Calling the old way raises rather than returning a different number. +def test_a_score_vector_is_rejected(emulsion_model: tuple[PCA, pd.DataFrame]) -> None: + """Passing a score vector raises rather than returning a different number. - A score vector cannot carry the information equation (3) needs, so the old - signature cannot be supported alongside the correct one. + Equation (3) needs ``x_ij``, which a score vector does not carry, so this + call cannot be served and must not appear to succeed. """ model, scaled = emulsion_model with pytest.raises(TypeError, match="not a score vector"): @@ -579,8 +572,7 @@ def test_contributions_are_not_proportional_to_the_loadings_on_real_data() -> No """The same guard, on the LDPE data set rather than a generated one. Real process data, so the test cannot be satisfied by a quirk of the - generator. Under the pre-1.61.0 calculation every one of these observations - produced the identical ranking of variables. + generator. """ import pathlib