kgdunn · kgdunn · Jul 25, 2026 · Jul 25, 2026 · Jul 25, 2026 · Jul 25, 2026
diff --git 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
@@ -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
@@ -10,15 +10,15 @@ 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:
 
 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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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,13 +46,15 @@
     "correlation_loadings_plot",
     "eigenvalue_summary",
     "explained_variance_plot",
+    "group_contributions",
     "loading_plot",
     "observation_contributions",
     "predictions_vs_observed_plot",
     "project_variables",
     "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
@@ -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
@@ -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)