diff --git a/CHANGELOG.md b/CHANGELOG.md index 7e561028..e78e309a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,61 @@ 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. 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 + 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 +2707,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..bd5472b0 100644 --- a/docs/user_guide/pca.rst +++ b/docs/user_guide/pca.rst @@ -167,31 +167,58 @@ 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 + +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. - # Why do observations 5 and 10 differ from each other? - contrib = model.score_contributions( - model.scores_.iloc[5].values, - model.scores_.iloc[10].values, +.. 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=X_scaled.index[64:74], reference=X_scaled.index[74:84] ) - # T²-weighted contributions (scale by 1/sqrt(eigenvalue)) - contrib = model.score_contributions(model.scores_.iloc[5].values, weighted=True) +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)``. +Those sum to the observation's :math:`T^2`. 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/__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..7e36fd3e 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,297 @@ 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] + + +_SCORE_CONTRIBUTIONS_USAGE = """\ +score_contributions() takes the preprocessed data X, not a score vector. + +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). + + model.score_contributions(X) # all observations, component 1 + model.score_contributions(X).iloc[i] # just observation i + +To compare two observations or two groups: + + model.group_contributions(X, group=[a], reference=[b]) + +For contributions to Hotelling's T-squared, which pools every component: + + model.t2_contributions(X) +""" + + +_SCORE_VECTOR_KEYWORDS = frozenset({"t_end", "components", "weighted"}) + + +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) & _SCORE_VECTOR_KEYWORDS): + raise TypeError(_SCORE_CONTRIBUTIONS_USAGE) + 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 + 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 + ------- + 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_score_vector_call(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( # noqa: PLR0913 - group/reference are sugar for weights + model: BaseEstimator, + X: DataMatrix, + 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`. 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. + + In general any linear combination of the rows may be used + (Miller, Swanson and Heckler, 1994): + + .. math:: + + c_k = \Bigl(\sum_i w_i x_{ik}\Bigr) R_{ka}, \qquad + \sum_k c_k = \sum_i w_i t_{ia}. + + 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 + ---------- + 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, 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. 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 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 + -------- + >>> 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? (ten batches either side, by position) + >>> 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 + -------- + score_contributions : The single-observation form. + """ + X_df, r_a = _score_contribution_terms(model, X, component) + + 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) + 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. + + 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): + 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 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: """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..bf86a909 100644 --- a/tests/test_multiblock_reference.py +++ b/tests/test_multiblock_reference.py @@ -379,25 +379,95 @@ 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) + 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) + 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[: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) - 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) + 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 @@ -1125,35 +1195,85 @@ 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(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_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_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_: - assert not np.allclose(partial[name].values, full[name].values) + 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.py b/tests/test_multivariate.py index ece114ec..1296f7ce 100644 --- a/tests/test_multivariate.py +++ b/tests/test_multivariate.py @@ -921,53 +921,141 @@ 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. + # 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 + + # --- 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_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")) + ) + 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_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) + + # 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") + + # 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) + 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 + ) + + # 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="not index labels of X"): + pca.group_contributions(X, group=["not-a-label"]) def test_pca_t2_spe_contributions() -> None: @@ -2471,7 +2559,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 +2569,24 @@ 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) - - # 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] - - # 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] + 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() + + # 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) + + 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) def test_pls_detect_outliers() -> None: @@ -2538,7 +2619,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 +2629,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( diff --git a/tests/test_multivariate_contribution_plots.py b/tests/test_multivariate_contribution_plots.py new file mode 100644 index 00000000..52cb2456 --- /dev/null +++ b/tests/test_multivariate_contribution_plots.py @@ -0,0 +1,683 @@ +# (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 +-------------------- +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 + +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). + """ + 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. + """ + 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. + + 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) + + 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()) + + +# --------------------------------------------------------------------------- +# Usage: the data are required, a score vector is not enough +# --------------------------------------------------------------------------- + + +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. + + 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"): + 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. + """ + 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)