# (c) Kevin Dunn, 2010-2026. MIT License. Based on own private work over the years.
"""Post-fit diagnostics and matrix-correlation helpers (ENG-01).
Functions that read a *fitted* PCA / PLS model and summarise it: variable
importance (VIP), squared cosine, observation contributions, the eigenvalue
summary, supplementary-variable projection, and the RV / modified-RV matrix
correlation coefficients. They consume model attributes only (via duck typing),
so they depend just on :mod:`process_improve.multivariate._common` and ``center``
from :mod:`process_improve.multivariate._preprocessing`, and annotate the model
parameter as ``BaseEstimator`` to avoid importing the estimator modules.
"""
from __future__ import annotations
from collections.abc import Callable, Sequence
import numpy as np
import pandas as pd
from scipy.cluster.hierarchy import leaves_list, linkage, optimal_leaf_ordering
from scipy.stats import f as f_dist
from sklearn.base import BaseEstimator
from sklearn.utils import Bunch
from .._linalg import safe_inverse
from ..univariate.metrics import _check_probability
from ._common import (
DataMatrix,
_align_to_fit_features,
_scale_block_contributions,
components_with_spread,
epsqrt,
)
from ._nipals import quick_regress, ssq
from ._preprocessing import center
from ._projection import project_rows
# These diagnostics operate on a *fitted* PCA or PLS model via duck typing (they
# read attributes such as ``scores_`` / ``spe_`` after guarding with
# ``hasattr``). They are annotated against ``BaseEstimator`` rather than the
# concrete ``PCA`` / ``PLS`` classes so this leaf module does not import the
# estimator modules - which import these functions in turn - and so avoids a
# module-level import cycle.
[docs]
def vip(model: BaseEstimator, n_components: int | None = None) -> pd.Series:
r"""Calculate Variable Importance in Projection (VIP) scores.
Works with fitted :class:`PCA` and :class:`PLS` models. For PCA the
principal-component loadings ``loadings_`` are used as the weight matrix;
for PLS the X-block weights ``x_weights_`` are used.
The formula is:
.. math::
\\text{VIP}_j = \\sqrt{K \\cdot
\\frac{\\sum_{a=1}^{A} r2_a \\cdot w_{ja}^2}{\\sum_{a=1}^{A} r2_a}}
where :math:`K` is the number of features, :math:`A` the number of
components, :math:`r2_a` the model's ``r2_per_component_`` at component
:math:`a`, and :math:`w_{ja}` the weight for feature :math:`j` in
component :math:`a`.
The same formula weights each component by a different block's variance on
the two models: :math:`r2_a` is the fraction of **X** variance component
:math:`a` explains for PCA (:math:`R^2X`), and the fraction of **Y**
variance for PLS (:math:`R^2Y`). A PLS VIP therefore ranks the variables by
how much they contribute to the components that predict Y, and a PCA VIP by
how much they contribute to the components that describe X; the two are
not comparable with each other.
Parameters
----------
model : PCA or PLS
A fitted PCA or PLS model.
n_components : int or None, default=None
Number of components to include. ``None`` uses all fitted components.
Returns
-------
pd.Series
VIP scores indexed by feature names, named ``"VIP"``.
Raises
------
ValueError
If the model is not fitted, if neither ``x_weights_`` nor
``loadings_`` is found, or if *n_components* is out of range.
Notes
-----
The ``K`` factor in the formula above normalises the scores so that
.. math::
\sum_{j=1}^{K} \text{VIP}_j^2 = K
*exactly*, for any model, on any data. That identity is what makes the
familiar "VIP > 1" rule a sensible relative cut-off: the mean square is 1
by construction, so a score above 1 means the variable is
above-average within this model.
It also means **the number of variables exceeding VIP 1 is not a test
statistic.** That count describes the *shape* of the VIP distribution, not
whether any relationship exists, and it barely moves when the response is
permuted: a null built on it has almost no power and will report a
false-discovery rate near 100% on data that genuinely contains signal. If
you want to ask "is there anything here at all", permute the response and
compare *out-of-sample* performance instead. See
:func:`~process_improve.multivariate.check_predictive_signal`.
See Also
--------
process_improve.multivariate.check_predictive_signal :
A permutation null on Q² that does respond to signal.
Examples
--------
>>> pls = PLS(n_components=3).fit(X_scaled, Y_scaled)
>>> pls.vip() # bound convenience method after fit()
>>> vip(pls) # or call the standalone function directly
>>> pca = PCA(n_components=3).fit(X_scaled)
>>> pca.vip(n_components=2)
"""
if not hasattr(model, "r2_per_component_"):
msg = "Model is not fitted. Call fit() before computing VIP."
raise ValueError(msg)
if hasattr(model, "x_weights_"):
weights: pd.DataFrame = model.x_weights_
elif hasattr(model, "loadings_"):
weights = model.loadings_
else:
msg = "Model must have 'x_weights_' (PLS) or 'loadings_' (PCA) to compute VIP."
raise ValueError(msg)
r2: np.ndarray = model.r2_per_component_.values
w: np.ndarray = weights.values # (n_features, total_components)
total_components = w.shape[1]
if n_components is not None and not (1 <= n_components <= total_components):
msg = f"n_components must be between 1 and {total_components}, got {n_components}."
raise ValueError(msg)
w = w[:, :n_components] # a slice up to None keeps every component
r2 = r2[:n_components]
n_features = w.shape[0]
vip_values = np.sqrt(n_features * np.sum(r2 * w**2, axis=1) / np.sum(r2))
return pd.Series(vip_values, index=weights.index, name="VIP")
def _select_response(beta: pd.DataFrame, response: str | int | None) -> str:
"""Resolve a response selector to a single column label of ``beta``."""
columns = list(beta.columns)
if response is None:
if len(columns) != 1:
msg = f"This model has several responses ({columns}); pass response=<name> to pick one."
raise ValueError(msg)
return columns[0]
if isinstance(response, (int, np.integer)) and response not in columns:
if not (0 <= int(response) < len(columns)):
msg = f"response index {response} is out of range for {len(columns)} responses."
raise ValueError(msg)
return columns[int(response)]
if response not in columns:
msg = f"response {response!r} is not one of the model responses {columns}."
raise ValueError(msg)
return response # type: ignore[return-value]
def _target_projection_arrays(
model: BaseEstimator, X: DataMatrix, response: str | int | None
) -> tuple[pd.DataFrame, np.ndarray, np.ndarray, np.ndarray, str]:
"""Return ``(X_aligned, t_tp, p_tp, w_tp, response_label)`` for the TP component.
The target-projection direction is the (unit-normalised) regression vector
``b`` for the chosen response, read from ``model.beta_coefficients_``; the
scores are ``t = X w_tp`` and the loadings ``p = X^T t / (t^T t)``. ``X``
is the data ``predict`` takes, moved into the model's scaled space first,
as :func:`t2_contributions` and :func:`spe_contributions` move it.
"""
if not hasattr(model, "beta_coefficients_"):
msg = "Model is not fitted, or is not a PLS model with 'beta_coefficients_'."
raise ValueError(msg)
beta = model.beta_coefficients_
label = _select_response(beta, response)
# The TP direction must live in the SAME space as the X the scores are
# computed from. ``beta_coefficients_`` is deliberately rescaled to raw
# X -> raw Y units in fit(); with scale=True (the default) that vector is
# NOT a direction in the model's internal (scaled) space, and normalising
# it does not repair the per-column unit mismatch. Use the scaled-space
# regression vector (direct_weights_ @ y_loadings_.T) instead, and map X
# into the same internal space via the model's own scaler when present.
b_scaled = model.direct_weights_.to_numpy() @ model.y_loadings_.to_numpy().T
b = b_scaled[:, list(beta.columns).index(label)]
norm_b = float(np.sqrt(b @ b))
if norm_b <= epsqrt:
msg = f"The regression vector for response {label!r} is ~0; it is not predicted by X."
raise ValueError(msg)
w_tp = b / norm_b
if not isinstance(X, pd.DataFrame):
X = pd.DataFrame(X)
x_frame: pd.DataFrame = _align_to_fit_features(X, beta.index)
x_scaler = getattr(model, "_x_scaler", None)
if x_scaler is not None:
x_frame = x_scaler.transform(x_frame)
X_values = x_frame.to_numpy(dtype=float)
# A missing cell made every score, loading and ratio NaN (#1114). It is skipped, as NIPALS
# skips it: a row's score is regressed on the observed part of w_tp (0 for a row with nothing
# observed), and a feature's loading on the scores of the rows where it is observed.
t_tp = quick_regress(X_values, w_tp[:, np.newaxis], per="row").ravel()
ttt = float(t_tp @ t_tp)
if ttt <= epsqrt:
msg = "The target-projected scores have ~0 variance; cannot form the TP loading."
raise ValueError(msg)
p_tp = quick_regress(X_values, t_tp[:, np.newaxis], per="column").ravel()
return x_frame, t_tp, p_tp, w_tp, label
def target_projection(model: BaseEstimator, X: DataMatrix, response: str | int | None = None) -> Bunch:
r"""Target-projected (TP) component of a fitted PLS model for one response.
Target projection (Kvalheim and Karstang, 1989) rotates the PLS solution so
that a *single* latent component carries all of the predictive information
for one response. The component points along the regression vector
:math:`b` (the column of ``beta_coefficients_`` for that response):
.. math::
w_{\text{TP}} = \frac{b}{\lVert b \rVert}, \qquad
t_{\text{TP}} = X\, w_{\text{TP}}, \qquad
p_{\text{TP}} = \frac{X^\top t_{\text{TP}}}{t_{\text{TP}}^\top t_{\text{TP}}}.
The TP component is the basis for the selectivity ratio
(:func:`selectivity_ratio`).
A missing cell (NaN) is skipped, as the NIPALS fit skips it: a row's score
is regressed on the observed part of :math:`w_{\text{TP}}`, and a
feature's loading on the scores of the rows where it is observed. A row
with nothing observed scores 0 and adds nothing, and a feature with
nothing observed has a loading of 0.
Parameters
----------
model : PLS
A fitted PLS model (must expose ``beta_coefficients_``).
X : array-like of shape (n_samples, n_features)
The data the model's ``predict`` takes, in the units it was fitted
on. A model that centres and scales ``X`` itself (``PLS(scale=True)``,
the default) does so here first; ``PLS(scale=False)`` does not, and is
given ``X`` scaled as its training data was.
response : str or int or None, default=None
Which response (Y column) to project onto. ``None`` is allowed only for
a single-response model; otherwise pass the response label (or its
integer position).
Returns
-------
sklearn.utils.Bunch
With fields ``scores`` (pd.Series, the TP scores per sample),
``loadings`` (pd.Series, the TP loading per feature), ``weights``
(pd.Series, the unit TP weight per feature) and ``response`` (the
resolved response label).
Raises
------
ValueError
If the model is not a fitted PLS, the response selector is invalid, or
the regression vector / TP scores are degenerate (~0).
References
----------
Kvalheim, O. M. and Karstang, T. V. (1989). Interpretation of latent-variable
regression models. *Chemometrics and Intelligent Laboratory Systems*, 7(1-2),
39-51.
Examples
--------
>>> pls = PLS(n_components=3).fit(X_scaled, y_scaled)
>>> tp = pls.target_projection(X_scaled) # bound convenience method
>>> tp.scores.head()
See Also
--------
selectivity_ratio : Per-variable explained/residual ratio on the TP component.
"""
X_df, t_tp, p_tp, w_tp, label = _target_projection_arrays(model, X, response)
return Bunch(
scores=pd.Series(t_tp, index=X_df.index, name="TP score"),
loadings=pd.Series(p_tp, index=X_df.columns, name="TP loading"),
weights=pd.Series(w_tp, index=X_df.columns, name="TP weight"),
response=label,
)
def _selectivity_ratio_one(
model: BaseEstimator, X: DataMatrix, response: str | int | None, conf_level: float
) -> pd.Series:
"""Compute the selectivity ratio per feature for a single response (public-API helper)."""
X_df, t_tp, p_tp, _w_tp, label = _target_projection_arrays(model, X, response)
X_values = X_df.to_numpy(dtype=float)
# Per-feature explained and residual sums of squares on the TP component, over the feature's
# observed cells: they add up to its observed sum of squares, as the loading is fitted there.
ss_explained = (p_tp**2) * ((t_tp**2) @ ~np.isnan(X_values))
ss_residual = ssq(X_values - np.outer(t_tp, p_tp), axis=0)
with np.errstate(invalid="ignore", divide="ignore"):
sr = ss_explained / ss_residual
# Rounding leaves a feature the TP component reproduces exactly a residual of about 1e-32 of its
# sum of squares, not 0. A residual of at most eps times the explained sum of squares (a norm at
# most epsqrt times the feature's) counts as none, so every ratio from 1 / eps up is inf.
sr[sr >= 1.0 / epsqrt**2] = np.inf
sr[np.isnan(sr)] = 0.0 # 0 / 0: no sum of squares, or a ratio that missing values leave undefined
series = pd.Series(sr, index=X_df.columns, name="selectivity_ratio")
# F-based critical value (Rajalahti et al., 2009): SR_j above this is
# "significant" at conf_level, with N-2 and N-3 degrees of freedom. The
# sensory layer prefers a permutation test, so this is advisory metadata.
# With N <= 3 a degree of freedom is below 1, and scipy returns NaN.
# N is the feature's own count of observed rows, over which its ratio is formed (#1144), so the
# value is one per feature, a Series whatever the counts: one type for the code that reads it.
n_observed = (~np.isnan(X_values)).sum(axis=0)
critical = f_dist.ppf(conf_level, dfn=n_observed - 2, dfd=n_observed - 3)
series.attrs["f_critical"] = pd.Series(critical, index=X_df.columns, name="f_critical")
series.attrs["conf_level"] = float(conf_level)
series.attrs["response"] = label
return series
def selectivity_ratio(
model: BaseEstimator,
X: DataMatrix,
response: str | int | None = None,
*,
conf_level: float = 0.95,
) -> pd.Series | pd.DataFrame:
r"""Compute the selectivity ratio of each feature on the target-projected component.
The selectivity ratio (Rajalahti et al., 2009) ranks each feature by how
much of its variance the *predictive* (target-projected) direction explains.
On the TP component (:func:`target_projection`), for feature :math:`j`:
.. math::
\text{SR}_j = \frac{\text{SS}_{\text{explained},j}}
{\text{SS}_{\text{residual},j}}
= \frac{p_{\text{TP},j}^2\, (t_{\text{TP}}^\top t_{\text{TP}})}
{\sum_i (x_{ij} - t_{\text{TP},i}\, p_{\text{TP},j})^2}.
A large SR means the feature is well aligned with the predictive direction.
Unlike VIP, it is a true explained/residual variance ratio and can be
compared against an F distribution. Note that two collinear features carry
near-identical SR: the selectivity ratio ranks predictive relevance, it does
not break ties between mutually collinear features.
A feature the TP component reproduces exactly has no residual, and its SR is
``inf``, so a descending sort ranks it first. The residual counts as none
when its sum of squares is at most machine epsilon (about 2.2e-16) times the
explained one, which is a residual norm of at most about 1.5e-8 times the
feature's: rounding alone leaves about 1e-32 of its sum of squares. Every
finite SR is therefore below ``1 / eps``, about 4.5e15. A feature with no
sum of squares at all, such as a constant column after :class:`MCUVScaler`,
has neither an explained nor a residual part (0 / 0), and its SR is 0.
A missing cell (NaN) is skipped, as in :func:`target_projection`, and both
sums of squares of a feature run over its observed cells, where they add
up to the feature's sum of squares. A feature with nothing observed has
none, and its SR is 0. The F critical value of a feature counts its
observed rows too.
Parameters
----------
model : PLS
A fitted PLS model (must expose ``beta_coefficients_``).
X : array-like of shape (n_samples, n_features)
The data the model's ``predict`` takes, in the units it was fitted
on. A model that centres and scales ``X`` itself (``PLS(scale=True)``,
the default) does so here first; ``PLS(scale=False)`` does not, and is
given ``X`` scaled as its training data was.
response : str or int or None, default=None
Which response to compute SR for. ``None`` returns a feature-by-response
DataFrame when the model has several responses, or a Series for a
single-response model.
conf_level : float, default=0.95
Confidence level for the advisory F-based critical value, attached to
the result's ``.attrs["f_critical"]``: a Series indexed by feature,
with :math:`N - 2` and :math:`N - 3` degrees of freedom, :math:`N`
being the number of rows the feature is observed in. So
``ratio > ratio.attrs["f_critical"]`` compares each feature with its own.
Returns
-------
pd.Series or pd.DataFrame
Selectivity ratios indexed by feature, each at least 0 and ``inf`` for
a feature with no residual. A Series for one response (with
``f_critical`` / ``conf_level`` / ``response`` in ``.attrs``), or a
feature-by-response DataFrame when ``response`` is ``None`` and the
model has several responses.
Raises
------
ValueError
If the model is not a fitted PLS or the response selector is invalid.
References
----------
Rajalahti, T., Arneberg, R., Berven, F. S., Myhr, K.-M., Ulvik, R. J. and
Kvalheim, O. M. (2009). Biomarker discovery in mass spectral profiles by
means of selectivity ratio plot. *Chemometrics and Intelligent Laboratory
Systems*, 95(1), 35-48.
Examples
--------
>>> pls = PLS(n_components=3).fit(X_scaled, y_scaled)
>>> pls.selectivity_ratio(X_scaled).sort_values(ascending=False).head()
See Also
--------
target_projection : The target-projected component the ratio is built on.
vip : Variable Importance in Projection, an alternative importance measure.
"""
_check_probability("conf_level", conf_level)
beta = getattr(model, "beta_coefficients_", None)
if response is None and beta is not None and beta.shape[1] > 1:
return pd.DataFrame({col: _selectivity_ratio_one(model, X, col, conf_level) for col in beta.columns})
return _selectivity_ratio_one(model, X, response, conf_level)
[docs]
def squared_cosine(model: BaseEstimator, n_components: int | None = None) -> pd.DataFrame:
r"""Calculate the squared cosine (cos2): quality of representation of observations.
Works with fitted :class:`PCA` and :class:`PLS` models. The squared cosine
of observation :math:`i` on component :math:`a` is the squared score
divided by that observation's total variation budget:
.. math::
\\cos^2_{ia} = \\frac{t_{ia}^2}
{\\sum_{a=1}^{A} t_{ia}^2 + \\text{SPE}_i^2}
where :math:`t_{ia}` is the score and :math:`\\text{SPE}_i` the residual
(squared prediction error) of the observation. Across all components the
cos2 values plus the residual fraction sum to 1. A value close to 1 means
the observation is well represented on that component. For :class:`PCA`,
whose loadings are orthonormal, the denominator equals the squared distance
of the observation from the origin, matching the classical definition.
cos2 complements the existing diagnostics: Hotelling's T² measures distance
*within* the model plane, SPE measures distance *to* it, and cos2 reports
how much of an observation's total variation a given component captures.
Parameters
----------
model : PCA or PLS
A fitted PCA or PLS model.
n_components : int or None, default=None
Number of components to return. ``None`` returns all fitted components.
Returns
-------
pd.DataFrame
cos2 values of shape (n_samples, n_components), indexed by sample.
Raises
------
ValueError
If the model is not fitted, or if *n_components* is out of range.
Examples
--------
>>> pca = PCA(n_components=3).fit(X_scaled)
>>> pca.squared_cosine() # bound convenience method after fit()
>>> squared_cosine(pca, n_components=2) # or call the function directly
"""
if not hasattr(model, "scores_") or not hasattr(model, "spe_"):
msg = "Model is not fitted. Call fit() before computing the squared cosine."
raise ValueError(msg)
scores = model.scores_
total_components = scores.shape[1]
if n_components is not None and not (1 <= n_components <= total_components):
msg = f"n_components must be between 1 and {total_components}, got {n_components}."
raise ValueError(msg)
score_ss = scores.to_numpy() ** 2
residual_ss = model.spe_.to_numpy()[:, -1] ** 2
total_ss = score_ss.sum(axis=1) + residual_ss
with np.errstate(invalid="ignore"): # 0 / 0 for a row at the origin; total_ss >= score_ss, so never x / 0
cos2 = score_ss[:, :n_components] / total_ss[:, None]
cos2[~np.isfinite(cos2)] = 0.0
return pd.DataFrame(cos2, index=scores.index, columns=scores.columns[:n_components])
[docs]
def observation_contributions(model: BaseEstimator, n_components: int | None = None) -> pd.DataFrame:
r"""Calculate the contribution of each observation to each component.
Works with fitted :class:`PCA` and :class:`PLS` models. The contribution of
observation :math:`i` to component :math:`a` is its squared score divided
by the sum of squared scores of all observations on that component:
.. math::
\\text{contribution}_{ia} = \\frac{t_{ia}^2}{\\sum_{i=1}^{N} t_{ia}^2}
Values lie between 0 and 1 and each column sums to 1, so a contribution
well above the average :math:`1/N` flags an observation that strongly
shapes that component. The exception is a component whose score column
has zero variance (``sum(t_{ia}^2) = 0``): the division cannot be
computed, so the column is returned as zeros rather than NaN, and its
sum is 0 rather than 1.
Note that this is *not* the same diagnostic as the ``score_contributions``
method, despite the similar name. ``score_contributions`` is *per-variable*
and signed: it decomposes one observation's position in score space back
onto the original variables ("which **variables** explain why this
observation sits where it does?"). ``observation_contributions`` is
*per-observation* and non-negative: it reports each observation's share of
a component's total inertia ("which **observations** most strongly shape
this component?"). The two are orthogonal views of the same score matrix
and are not interchangeable.
Parameters
----------
model : PCA or PLS
A fitted PCA or PLS model.
n_components : int or None, default=None
Number of components to return. ``None`` returns all fitted components.
Returns
-------
pd.DataFrame
Contributions of shape (n_samples, n_components), indexed by sample.
Each column sums to 1 except columns for components whose score has
zero variance, which are returned as zeros and sum to 0.
Raises
------
ValueError
If the model is not fitted, or if *n_components* is out of range.
Examples
--------
>>> pca = PCA(n_components=3).fit(X_scaled)
>>> pca.observation_contributions()
>>> observation_contributions(pca, n_components=2)
See Also
--------
PCA.score_contributions : The per-variable counterpart - decomposes one
observation's score-space position back onto the original variables.
"""
if not hasattr(model, "scores_"):
msg = "Model is not fitted. Call fit() before computing observation contributions."
raise ValueError(msg)
scores = model.scores_
total_components = scores.shape[1]
if n_components is not None and not (1 <= n_components <= total_components):
msg = f"n_components must be between 1 and {total_components}, got {n_components}."
raise ValueError(msg)
score_ss = scores.to_numpy() ** 2
column_ss = score_ss.sum(axis=0)
with np.errstate(invalid="ignore"): # 0 / 0 for a score with no spread; column_ss >= score_ss, so never x / 0
contributions = score_ss[:, :n_components] / column_ss[:n_components]
contributions[~np.isfinite(contributions)] = 0.0
return pd.DataFrame(contributions, index=scores.index, columns=scores.columns[:n_components])
def _contribution_inputs(model: BaseEstimator, X: DataMatrix) -> tuple[pd.DataFrame, np.ndarray, np.ndarray]:
"""Return ``(X_rows, R, P)``: ``X`` in the space the model was fitted in, and the model's matrices.
``X`` is the data ``transform`` and ``diagnose`` take, and is read through the
model's own ``_scaled_rows``, as they read it (#1125): checked against the fit,
its columns in the fitted order, centred and scaled when the model does that
itself (``PLS(scale=True)``), and 0 in a column the fit left out (#1098), so
that column contributes 0 (#1114).
``R`` generates the scores (``T = X @ R``) and ``P`` reconstructs the X-block
(``X_hat = T @ P.T``). For PCA both are the loadings; for PLS ``R`` is the
direct (rotated) weights ``direct_weights_`` and ``P`` is ``x_loadings_``.
PLS is recognised by the presence of ``direct_weights_``. OPLS, which has
neither, gives its own through ``_contribution_matrices`` (#1144).
"""
if not hasattr(model, "scores_"):
msg = "Model is not fitted. Call fit() before computing contributions."
raise ValueError(msg)
if hasattr(model, "_contribution_matrices"):
R, P = model._contribution_matrices()
else:
is_pls = hasattr(model, "direct_weights_")
P = np.asarray(model.x_loadings_ if is_pls else model.loadings_)
R = np.asarray(model.direct_weights_ if is_pls else model.loadings_)
# An array's columns are taken by position, labelled as fitted, as they always were here: a frame
# lets _scaled_rows do that without sklearn's warning that the rows carry no feature names. A frame
# passes through pd.DataFrame with its labels and values as they are.
return model._scaled_rows(pd.DataFrame(X)), R, P
def _complete_rows_apart(X_values: np.ndarray, rows_product: Callable[[np.ndarray], np.ndarray]) -> np.ndarray:
"""Evaluate ``rows_product(rows)`` for the complete rows and for the rest, separately.
The last bit of a matrix product can depend on how many rows it has, since BLAS
splits the work by shape (Apple's Accelerate does). Multiplying the complete rows on
their own makes them match, bit for bit, what the same rows give when ``X`` has no
missing cells, so a complete row's contributions do not depend on the other rows.
"""
incomplete = np.isnan(X_values).any(axis=1)
out = np.empty(X_values.shape)
for rows in (~incomplete, incomplete):
out[rows] = rows_product(rows)
return out
def _scores_and_guides(
model: BaseEstimator, X_values: np.ndarray, R: np.ndarray, P: np.ndarray, method: str | None, ridge: float
) -> tuple[np.ndarray, np.ndarray | None]:
"""Return the scores of every row and, when rows have missing cells, the per-row score-generating matrices.
A complete row's score is ``x @ R``. A row with missing cells has its scores
estimated from the observed cells by
:func:`~process_improve.multivariate._projection.project_rows`, the code
``project``, ``transform`` and ``diagnose`` score it with, through the
operator ``M`` of its pattern (``t_hat = M @ x_observed``); its
contributions use ``M`` where a complete row uses ``R``. ``method=None``
takes the estimator the model's ``transform`` and ``diagnose`` use, so the
contributions split the scores they give (#1144). A row whose estimator
cannot be inverted is NaN, with a warning, as there. The second return value
is ``None`` when every row is complete, and otherwise an array of shape
(n_samples, n_features, n_components) holding ``R`` for the complete rows
and, for the others, ``M.T`` on the observed features and zero on the
missing ones.
"""
if not np.isnan(X_values).any():
return X_values @ R, None
estimated = project_rows(
P,
R,
np.asarray(model.explained_variance_),
X_values,
method=getattr(
model,
"_MISSING_CELL_ESTIMATOR",
"tsr", # pragma: no mutate - every estimator of the package sets it, so this default is not read
)
if method is None
else method,
ridge=ridge,
# SCP projects onto the weights, as NIPALS does; for PCA they are the loadings (the operator's default).
x_weights=np.asarray(model.x_weights_) if hasattr(model, "direct_weights_") else None,
x_residuals=getattr(model, "_x_residuals", None),
singular="nan", # pragma: no mutate - project_rows reads any value but "raise" as "nan"
with_guides=True,
)
return estimated.scores, estimated.guides
def _t2_score_sd(score_sd: np.ndarray) -> np.ndarray:
"""Return the score standard deviations Hotelling's T² divides by, ``inf`` for a component it leaves out.
A component with no spread, by :func:`~process_improve.multivariate._common.components_with_spread`,
carries no T²: a divisor of ``inf`` makes it add 0 for any finite score. Every T² of PCA,
PLS, MBPCA and MBPLS goes through this one rule, so a row gets the same T² from ``fit`` as
from ``diagnose`` or ``project``.
Parameters
----------
score_sd : np.ndarray of shape (n_components,)
The training score standard deviations (``scaling_factor_for_scores_``).
Returns
-------
np.ndarray of shape (n_components,)
``score_sd``, with ``inf`` for each component left out.
"""
return np.where(components_with_spread(score_sd), score_sd, np.inf)
def t2_precision(covariance: np.ndarray, score_sd: np.ndarray, *, what: str) -> np.ndarray:
"""Return the matrix ``S^-1`` of the T² ``t' S^-1 t`` for a full score covariance ``S`` (#1124).
A component with no spread, by the rule :func:`_t2_score_sd` applies to the training score
standard deviations ``score_sd``, gets a zero row and column, so it adds 0 to T², and the
inverse is taken over the other components only: their block of ``S``. The whole of ``S``
is singular when a component has no spread, since its scores are all about 0.
Parameters
----------
covariance : np.ndarray of shape (n_components, n_components)
The score covariance ``S``.
score_sd : np.ndarray of shape (n_components,)
The model's training score standard deviations, which decide the components with spread.
what : str
What ``S`` is, for the error :func:`~process_improve._linalg.safe_inverse` raises when the
block of the components with spread is itself singular.
Returns
-------
np.ndarray of shape (n_components, n_components)
"""
keep = np.flatnonzero(components_with_spread(score_sd))
precision = np.zeros_like(covariance, dtype=float)
precision[np.ix_(keep, keep)] = safe_inverse(covariance[np.ix_(keep, keep)], what=what)
return precision
def cumulative_block_t2(scores: np.ndarray) -> np.ndarray:
"""Return the cumulative Hotelling's T² of each row of a multiblock model's block scores (#513).
Each component is scaled by the scores' second moment about zero,
``sum(t**2) / (n - 1)`` over the ``n`` rows that have a score, which is how
``explained_variance_`` scales the super scores, so a block's T² and the super T²
follow one convention. T² measures the distance from the model's centre, the
origin of the score space, so the spread it is scaled by is the spread about that
origin too. The block scores used to be scaled by their variance about their own
mean instead, which differs as soon as missing data leaves them off-centre, and
which a single row without a score turned into NaN, so that every row's T² in the
block went unscaled. A component with no spread in the block is left out by
:func:`_t2_score_sd`, as from every other T² (#1028).
Parameters
----------
scores : np.ndarray of shape (n_samples, n_components)
One block's scores, NaN for a row with nothing observed in the block.
Returns
-------
np.ndarray of shape (n_samples, n_components)
Column ``a`` holds the T² over the first ``a`` components; NaN for a row with
no score. A component with no spread adds nothing.
"""
present = np.isfinite(scores)
second_moment = np.nansum(scores**2, axis=0) / np.maximum(1, present.sum(axis=0) - 1)
return np.cumsum((scores**2) / _t2_score_sd(np.sqrt(second_moment)) ** 2, axis=1)
def t2_contributions(
model: BaseEstimator,
X: DataMatrix,
components: list[int] | None = None,
*,
method: str | None = None,
ridge: float = 0.0,
) -> pd.DataFrame:
r"""Per-variable contributions to Hotelling's :math:`T^2`.
Works with fitted :class:`PCA` and :class:`PLS` models. Decomposes each
observation's :math:`T^2` onto the original variables. The contribution of
variable :math:`k` for observation :math:`i` is
.. math::
c^{T^2}_{ik} = x_{ik} \sum_{a} \frac{t_{ia}}{s_a^2}\, R_{ka},
where :math:`t_{ia}` are the scores, :math:`s_a^2` is the score variance of
component :math:`a` (``scaling_factor_for_scores_`` squared) and :math:`R` is
the score-generating matrix (loadings for PCA, ``direct_weights_`` for PLS,
so that :math:`T = XR`). Summed over the variables this telescopes to
:math:`\sum_a t_{ia}^2 / s_a^2`, i.e. the observation's :math:`T^2`. The
values are signed; a large magnitude flags a variable that drives the
observation away from the model centre. This is the standard MSPC
diagnostic (Westerhuis, Gurden and Smilde, 2000). A component whose score
standard deviation is ~0 contributes nothing, as it adds nothing to the
stored ``hotellings_t2_``.
A row with missing cells (NaN) gets its scores from its observed cells
alone, with the estimator named by ``method``, the missing-data operators
of :meth:`PCA.project` and :meth:`PLS.project`. Its contributions are then
defined at every observed cell and NaN at the missing ones, which add
nothing to the row sums (``sum(axis=1)`` skips NaN). The default estimator
is the one the model's ``transform`` and ``diagnose`` score the row with,
``"tsr"`` for PCA and PLS and ``"scp"`` for the adaptive models, so the
contributions split the scores, T² and SPE they give it (#1144).
``method="scp"`` is the single-component projection the NIPALS fit uses:
on a model fitted with missing data, the training rows then reproduce the
stored ``scores_``, ``hotellings_t2_`` and ``spe_`` to within the NIPALS
convergence tolerance, as the complete rows do. A column the fit left out,
having no observed value, is read at 0, as :meth:`PCA.transform` reads it,
so it contributes 0 whatever ``X`` holds there.
Parameters
----------
model : PCA or PLS
A fitted PCA or PLS model.
X : array-like of shape (n_samples, n_features)
The data the model's ``transform`` and ``diagnose`` take, in the units
it was fitted on. A model that centres and scales ``X`` itself
(``PLS(scale=True)``, the default) does so here first; ``PCA`` and
``PLS(scale=False)`` do not, and are given ``X`` scaled as their
training data was, for example with :class:`MCUVScaler`.
Passing the training data reproduces the model's stored ``hotellings_t2_``
(with ``method="scp"`` when it has missing cells).
Until #1125 the model's own scaling was skipped here: code that
scaled ``X`` by hand for a ``scale=True`` model now scales it twice,
and passes the unscaled data instead.
components : list of int, optional
**1-based** component indices to decompose over, matching the model's
column convention. ``None`` (default) uses all fitted components, so the
row sums equal the cumulative :math:`T^2`.
method : {"tsr", "scp", "pmp"} or None, default=None
Score estimator for the rows with missing cells; ignored when ``X`` is
complete. ``None`` takes the one the model's ``transform`` and
``diagnose`` use. See :meth:`PCA.project`.
ridge : float, default=0.0
Regularisation for the ``"tsr"`` and ``"pmp"`` estimators, as in
:meth:`PCA.project`.
Returns
-------
pd.DataFrame
Signed contributions of shape (n_samples, n_features), NaN at the
missing cells. Each row sums to the observation's :math:`T^2` over the
selected components.
Examples
--------
>>> pca = PCA(n_components=3).fit(X_scaled)
>>> contrib = pca.t2_contributions(X_scaled)
>>> contrib.sum(axis=1) # equals pca.hotellings_t2_.iloc[:, -1]
See Also
--------
spe_contributions : The residual-space counterpart.
PCA.score_contributions : Decomposes a single score-space movement.
"""
X_df, R, P = _contribution_inputs(model, X)
X_values = np.ascontiguousarray(X_df.to_numpy(dtype=float)) # one memory layout, whatever the input's
A = R.shape[1]
if components is None:
idx = np.arange(A)
else:
idx = np.asarray(components, dtype=int) - 1
if idx.size == 0 or idx.min() < 0 or idx.max() >= A:
msg = f"components must be 1-based indices within 1..{A}, got {components}."
raise ValueError(msg)
# A component whose score sd is ~0 carries no T2: fit() leaves it out of
# hotellings_t2_ by this same rule, and it contributes 0 here too, whatever X
# scores on it, rather than inf or NaN.
s2 = _t2_score_sd(np.asarray(model.scaling_factor_for_scores_))[idx] ** 2
scores, guides = _scores_and_guides(model, X_values, R, P, method, ridge)
weights = scores[:, idx] / s2 # (n, len(idx))
factors = _complete_rows_apart(X_values, lambda rows: weights[rows] @ R[:, idx].T)
if guides is not None:
rows = np.flatnonzero(np.isnan(X_values).any(axis=1))
factors[rows] = np.einsum(
"ia,ika->ik", # pragma: no mutate - einsum reads upper-case subscripts as lower-case ones
weights[rows],
guides[rows][:, :, idx],
)
contributions = X_values * factors
return pd.DataFrame(contributions, index=X_df.index, columns=X_df.columns)
def spe_contributions(
model: BaseEstimator,
X: DataMatrix,
*,
method: str | None = None,
ridge: float = 0.0,
) -> pd.DataFrame:
r"""Per-variable squared-prediction-error (SPE / DModX) contributions.
Works with fitted :class:`PCA` and :class:`PLS` models. Returns the signed
residual of each variable after reconstructing the X-block from the full
model:
.. math::
e_{ik} = x_{ik} - \hat{x}_{ik}, \qquad \hat{X} = T P^\top,
where :math:`P` is the reconstruction loadings (``loadings_`` for PCA,
``x_loadings_`` for PLS). The squared residuals sum across variables to the
observation's SPE; equivalently ``(spe_contributions(X) ** 2).sum(axis=1)``
equals the stored ``spe_`` (final column) squared. The signs show whether a
variable sits above or below its reconstruction, which is the standard SPE
contribution plot used to diagnose why an observation has a high residual.
A row with missing cells (NaN) gets its scores from its observed cells
alone, with the estimator named by ``method``, the missing-data operators
of :meth:`PCA.project` and :meth:`PLS.project`. Its contributions are then
defined at every observed cell and NaN at the missing ones, which add
nothing to the row sums (``sum(axis=1)`` skips NaN). The default estimator
is the one the model's ``transform`` and ``diagnose`` score the row with,
``"tsr"`` for PCA and PLS and ``"scp"`` for the adaptive models, so the
contributions split the scores, T² and SPE they give it (#1144).
``method="scp"`` is the single-component projection the NIPALS fit uses:
on a model fitted with missing data, the training rows then reproduce the
stored ``scores_``, ``hotellings_t2_`` and ``spe_`` to within the NIPALS
convergence tolerance, as the complete rows do. A column the fit left out,
having no observed value, is read at 0, as :meth:`PCA.transform` reads it,
so it contributes 0 whatever ``X`` holds there.
Parameters
----------
model : PCA or PLS
A fitted PCA or PLS model.
X : array-like of shape (n_samples, n_features)
The data the model's ``transform`` and ``diagnose`` take, in the units
it was fitted on. A model that centres and scales ``X`` itself
(``PLS(scale=True)``, the default) does so here first; ``PCA`` and
``PLS(scale=False)`` do not, and are given ``X`` scaled as their
training data was, for example with :class:`MCUVScaler`.
Passing the training data reproduces the model's stored ``spe_``
(with ``method="scp"`` when it has missing cells).
Until #1125 the model's own scaling was skipped here: code that
scaled ``X`` by hand for a ``scale=True`` model now scales it twice,
and passes the unscaled data instead.
method : {"tsr", "scp", "pmp"} or None, default=None
Score estimator for the rows with missing cells; ignored when ``X`` is
complete. ``None`` takes the one the model's ``transform`` and
``diagnose`` use. See :meth:`PCA.project`.
ridge : float, default=0.0
Regularisation for the ``"tsr"`` and ``"pmp"`` estimators, as in
:meth:`PCA.project`.
Returns
-------
pd.DataFrame
Signed per-variable residuals of shape (n_samples, n_features), NaN at
the missing cells. The squared row sums equal the observation's SPE
squared, over its observed cells.
Examples
--------
>>> pca = PCA(n_components=2).fit(X_scaled)
>>> resid = pca.spe_contributions(X_scaled)
>>> (resid ** 2).sum(axis=1) # equals pca.spe_.iloc[:, -1] ** 2
See Also
--------
t2_contributions : The :math:`T^2` (score-space) counterpart.
"""
X_df, R, P = _contribution_inputs(model, X)
X_values = np.ascontiguousarray(X_df.to_numpy(dtype=float)) # one memory layout, whatever the input's
scores, _guides = _scores_and_guides(model, X_values, R, P, method, ridge)
residuals = _complete_rows_apart(X_values, lambda rows: X_values[rows] - scores[rows] @ P.T)
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 data X, as transform does, 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)
[docs]
def score_contributions(
model: BaseEstimator,
X: DataMatrix,
component: int = 1,
scaling: str = "none",
*,
method: str | None = None,
ridge: float = 0.0,
**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.
A row with missing cells (NaN) gets its scores from its observed cells
alone, with the estimator named by ``method``, the missing-data operators
of :meth:`PCA.project` and :meth:`PLS.project`. Its contributions are then
defined at every observed cell and NaN at the missing ones, which add
nothing to the row sums (``sum(axis=1)`` skips NaN). The default estimator
is the one the model's ``transform`` and ``diagnose`` score the row with,
``"tsr"`` for PCA and PLS and ``"scp"`` for the adaptive models, so the
contributions split the scores, T² and SPE they give it (#1144).
``method="scp"`` is the single-component projection the NIPALS fit uses:
on a model fitted with missing data, the training rows then reproduce the
stored ``scores_``, ``hotellings_t2_`` and ``spe_`` to within the NIPALS
convergence tolerance, as the complete rows do. A column the fit left out,
having no observed value, is read at 0, as :meth:`PCA.transform` reads it,
so it contributes 0 whatever ``X`` holds there.
Parameters
----------
model : PCA or PLS
A fitted PCA or PLS model.
X : array-like of shape (n_samples, n_features)
The data the model's ``transform`` and ``diagnose`` take, in the units
it was fitted on. A model that centres and scales ``X`` itself
(``PLS(scale=True)``, the default) does so here first; ``PCA`` and
``PLS(scale=False)`` do not, and are given ``X`` scaled as their
training data was, for example with :class:`MCUVScaler`.
Passing the training data reproduces the model's stored ``scores_``
(with ``method="scp"`` when it has missing cells).
Until #1125 the model's own scaling was skipped here: code that
scaled ``X`` by hand for a ``scale=True`` model now scales it twice,
and passes the unscaled data instead.
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.
method : {"tsr", "scp", "pmp"} or None, default=None
Score estimator for the rows with missing cells; ignored when ``X`` is
complete. ``None`` takes the one the model's ``transform`` and
``diagnose`` use. See :meth:`PCA.project`.
ridge : float, default=0.0
Regularisation for the ``"tsr"`` and ``"pmp"`` estimators, as in
:meth:`PCA.project`.
**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), NaN at the
missing cells. 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)
X_values = X_df.to_numpy(dtype=float)
_, R, P = _contribution_inputs(model, X)
_scores, guides = _scores_and_guides(model, X_values, R, P, method, ridge)
contributions = X_values * (r_a if guides is None else guides[:, :, int(component) - 1])
# A missing cell's contribution is NaN, so both scalings skip NaN (as the row sums do). One block is the
# multi-block case with a single block, so the floors are the same in both (#1070).
scaled = _scale_block_contributions({"X": contributions}, scaling)["X"]
return pd.DataFrame(scaled, index=X_df.index, columns=X_df.columns)
[docs]
def group_contributions(
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.
A row with missing cells takes part with the contributions
:func:`score_contributions` gives it, from its observed cells alone. A
missing cell has no part in its row's score, so it adds nothing here,
and the contributions still sum to the weighted combination of the
scores, now the estimated ones for such a row.
Parameters
----------
model : PCA or PLS
A fitted PCA or PLS model.
X : array-like of shape (n_samples, n_features)
The data the model's ``transform`` and ``diagnose`` take, in the units
it was fitted on. A model that centres and scales ``X`` itself
(``PLS(scale=True)``, the default) does so here first; ``PCA`` and
``PLS(scale=False)`` do not, and are given ``X`` scaled as their
training data was, for example with :class:`MCUVScaler`.
Until #1125 the model's own scaling was skipped here: code that
scaled ``X`` by hand for a ``scale=True`` model now scales it twice,
and passes the unscaled data instead.
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.
"""
# Each row's contributions, from its observed cells alone when it has a missing one, as
# score_contributions gives them; a missing cell is NaN there and adds nothing to its row's
# score, so it adds nothing to the combination either, rather than making its column NaN.
per_row = score_contributions(model, X, component=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(per_row):
msg = f"weights must have one entry per row of X ({len(per_row)}), got shape {w.shape}."
raise ValueError(msg)
else:
if group is None:
msg = "Pass group (optionally with reference), or weights."
raise ValueError(msg)
w = _mean_weights(per_row, group, "group")
if reference is not None:
w -= _mean_weights(per_row, reference, "reference")
values = per_row.to_numpy()
return pd.Series(w @ np.where(np.isnan(values), 0.0, values), index=per_row.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))
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)]
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 _complete_through_model(
model: BaseEstimator, X: DataMatrix
) -> tuple[pd.DataFrame, np.ndarray, np.ndarray, np.ndarray]:
"""Align ``X`` to the fit and return ``(X_aligned, x, R, P)`` with every missing cell filled.
A missing cell is replaced by the model's reconstruction of it, from the scores that
the row's observed cells give by single-component projection, as the NIPALS fit scores
such a row, so the cross-products MEDA and oMEDA are built from are of complete data.
"""
X_df, R, P = _contribution_inputs(model, X)
x = X_df.to_numpy(dtype=float)
missing = np.isnan(x)
if missing.any():
scores, _ = _scores_and_guides(
model,
x,
R,
P,
method="scp", # pragma: no mutate - _validate_method lowercases it, so "SCP" is the same
ridge=0.0, # pragma: no mutate - SCP does not use the ridge, so 1.0 is the same
)
x = np.where(missing, scores @ P.T, x)
return X_df, x, R, P
def _meda_index(cross: np.ndarray, modelled: np.ndarray, scale: np.ndarray, *, signed: bool) -> np.ndarray:
"""Return ``(2 cross - modelled) * modelled / scale``, with ``|modelled|`` when signed, and 0 where scale is 0."""
weight = np.abs(modelled) if signed else modelled
# A variable with no sum of squares makes scale 0; a divisor of inf gives its entries 0, without a warning.
return (2.0 * cross - modelled) * weight / np.where(scale > 0, scale, np.inf)
def _seriation_order(matrix: np.ndarray) -> np.ndarray:
"""Order rows and columns so that strongly related variables sit next to each other.
Average-linkage clustering on the distance ``1 - |m|`` (symmetrised), with the
leaves put in the order that maximises the similarity of neighbours (Bar-Joseph,
Gifford and Jaakkola, 2001), so each block of related variables is contiguous.
Two variables have only one order, up to reversal.
"""
n = len(matrix)
if n < 3:
return np.arange(n)
similarity = np.abs((matrix + matrix.T) / 2.0)
distance = np.clip(1.0 - similarity, 0.0, None)
condensed = distance[np.triu_indices(n, 1)] # scipy's condensed form: the upper triangle, row by row
return leaves_list(optimal_leaf_ordering(linkage(condensed, method="average"), condensed))
[docs]
def meda(model: BaseEstimator, X: DataMatrix, *, signed: bool = True, seriate: bool = False) -> pd.DataFrame:
r"""Map how the fitted model relates every pair of variables (MEDA).
MEDA, missing-data methods for exploratory data analysis (Camacho, 2010), asks
for each pair of variables how well one is predicted from the other alone,
through the model. A correlation matrix answers that from all the directions in
the data; MEDA answers it from the model's components only. It therefore shows
the relationships the model represents, and adding a component shows which
relationships that component brings in.
Entry :math:`(j, k)` is the goodness of prediction of variable :math:`k` when
only variable :math:`j` is known and the scores are estimated from it by
regression. With :math:`S = X^\top X` and :math:`S_A = S R P^\top`, the part
of it the model reproduces (``R`` gives the scores, ``P`` the reconstruction:
both the loadings for PCA, the direct weights and X-loadings for PLS),
.. math::
q^2_{jk} = 1 - \frac{\lVert x_k - \hat{x}_{k \mid j}\rVert^2}{\lVert x_k \rVert^2}
= \frac{(2 S_{jk} - S_{A,jk})\, S_{A,jk}}{S_{jj} S_{kk}},
\qquad \hat{x}_{k \mid j} = x_j \frac{S_{A,jk}}{S_{jj}}.
The signed index, the default, replaces the last :math:`S_{A,jk}` by its
absolute value, as the MEDA Toolbox does (Camacho et al., 2015). The size is
unchanged where the model reproduces the relationship, and the sign then says
whether the two variables move together or in opposition. The diagonal is
:math:`r_k (2 - r_k)` for the fraction :math:`r_k` of variable :math:`k`'s sum
of squares that the model reproduces, which lies in :math:`[0, 1]`.
Parameters
----------
model : PCA or PLS
A fitted PCA or PLS model. Every fitted component is used; to see the map
at fewer components, fit a model with fewer.
X : array-like of shape (n_samples, n_features)
The data the model's ``transform`` and ``diagnose`` take, in the units
it was fitted on. A model that centres and scales ``X`` itself
(``PLS(scale=True)``, the default) does so here first; ``PCA`` and
``PLS(scale=False)`` do not, and are given ``X`` scaled as their
training data was, for example with :class:`MCUVScaler`. Usually the
training data itself. A missing cell is filled with the model's
reconstruction of it before the cross-products are formed.
Until #1125 the model's own scaling was skipped here: code that
scaled ``X`` by hand for a ``scale=True`` model now scales it twice,
and passes the unscaled data instead.
signed : bool, default=True
Return the signed index. ``False`` returns the goodness of prediction
:math:`q^2_{jk}` itself.
seriate : bool, default=False
Reorder the rows and columns so that groups of related variables are
contiguous (hierarchical clustering on :math:`1 - |\text{MEDA}|` with
optimal leaf ordering). The labels travel with the values.
Returns
-------
pd.DataFrame
The K x K map, rows and columns labelled by variable. It is symmetric
for PCA, and close to symmetric for PLS.
Examples
--------
>>> pca = PCA(n_components=3).fit(X_scaled)
>>> pca.meda(X_scaled, seriate=True)
>>> pca.meda_plot(X_scaled)
See Also
--------
omeda : Which variables separate one group of observations from another.
References
----------
Camacho, J. (2010). Missing-data theory in the context of exploratory data
analysis. Chemometrics and Intelligent Laboratory Systems 103(1), 8-18.
Camacho, J., Perez-Villegas, A., Rodriguez-Gomez, R. A. and Jimenez-Manas, E.
(2015). Multivariate Exploratory Data Analysis (MEDA) Toolbox for Matlab.
Chemometrics and Intelligent Laboratory Systems 143, 49-57.
"""
X_df, x, R, P = _complete_through_model(model, X)
cross = x.T @ x
modelled = cross @ R @ P.T
sums_of_squares = np.diag(cross)
values = _meda_index(cross, modelled, np.outer(sums_of_squares, sums_of_squares), signed=signed)
names = X_df.columns
if seriate:
order = _seriation_order(values)
values, names = values[np.ix_(order, order)], names[order]
return pd.DataFrame(values, index=names, columns=names)
def _omeda_dummy(
X_df: pd.DataFrame, group: Sequence | None, reference: Sequence | None, weights: Sequence | None
) -> np.ndarray:
"""Build oMEDA's dummy: +1 on ``group`` and -1 on ``reference``, or ``weights``; each side peaks at 1 in size."""
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)
dummy = np.asarray(weights, dtype=float).copy()
if dummy.ndim != 1 or dummy.size != len(X_df):
msg = f"weights must have one entry per row of X ({len(X_df)}), got shape {dummy.shape}."
raise ValueError(msg)
# Each side is scaled so its largest weight is 1 in size (Camacho, 2011), as
# the +1 / -1 dummy of group and reference below already is.
positive, negative = dummy > 0, dummy < 0
if positive.any():
dummy[positive] /= dummy[positive].max()
if negative.any():
dummy[negative] /= -dummy[negative].min()
else:
if group is None:
msg = "Pass group (optionally with reference), or weights."
raise ValueError(msg)
dummy = np.zeros(len(X_df))
in_group = X_df.index.get_indexer_for(_select_rows(X_df, group, "group").index)
dummy[in_group] = 1.0
if reference is not None:
in_reference = X_df.index.get_indexer_for(_select_rows(X_df, reference, "reference").index)
if np.intersect1d(in_group, in_reference).size:
msg = "group and reference share observations; each observation belongs to one side."
raise ValueError(msg)
dummy[in_reference] = -1.0
if not np.any(dummy):
msg = "The comparison selects no observations: every weight is zero."
raise ValueError(msg)
return dummy
[docs]
def omeda(
model: BaseEstimator,
X: DataMatrix,
group: Sequence | None = None,
reference: Sequence | None = None,
*,
weights: Sequence | None = None,
) -> pd.Series:
r"""Find the variables that separate one group of observations from another (oMEDA).
Observation-based MEDA (Camacho, 2011) answers the question a pair of clusters
on a score plot poses: which variables put these observations here and those
there? It compares the groups only along the directions the model trusts, so a
variable that differs between them but that the model does not capture scores
low.
A dummy vector :math:`d` marks the comparison: :math:`+1` on ``group``,
:math:`-1` on ``reference`` and 0 elsewhere, or ``weights`` with each side
rescaled so its largest weight is 1 in size. With :math:`\hat{X} = X R P^\top`
the data reconstructed through the model,
.. math::
\text{oMEDA}_k = \frac{(2\, x_k^\top d - \hat{x}_k^\top d)\,
\lvert \hat{x}_k^\top d \rvert}{\sqrt{d^\top d}}.
This is the MEDA index with the dummy in place of a second variable. A
positive value means the variable is higher in ``group`` than in ``reference``
(or than the model centre, with no reference), a negative value that it is
lower. The size grows with the difference and with how well the model
reproduces it. The dummy sums the observations rather than averaging them, so
a comparison of unequal groups leans towards the larger one, as in the
original method.
Parameters
----------
model : PCA or PLS
A fitted PCA or PLS model; every fitted component is used.
X : array-like of shape (n_samples, n_features)
The data the model's ``transform`` and ``diagnose`` take, in the units
it was fitted on. A model that centres and scales ``X`` itself
(``PLS(scale=True)``, the default) does so here first; ``PCA`` and
``PLS(scale=False)`` do not, and are given ``X`` scaled as their
training data was, for example with :class:`MCUVScaler`. A missing
cell is filled with the model's reconstruction of it.
Until #1125 the model's own scaling was skipped here: code that
scaled ``X`` by hand for a ``scale=True`` model now scales it twice,
and passes the unscaled data instead.
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, or a boolean mask, selecting the observations to compare
against. ``None`` (default) compares the group against the model centre.
weights : sequence, optional
One weight per row of ``X``: a general dummy vector, for comparisons
other than two groups. Mutually exclusive with ``group`` / ``reference``.
Returns
-------
pd.Series
One signed value per variable, named ``"omeda"``.
Examples
--------
>>> pca = PCA(n_components=3).fit(X_scaled)
>>> # The cluster at the top right of the score plot, against the rest:
>>> rest = X_scaled.index.difference(cluster)
>>> pca.omeda(X_scaled, group=cluster, reference=rest).sort_values()
>>> pca.omeda_plot(X_scaled, group=cluster, reference=rest)
See Also
--------
group_contributions : The contributions to one component's score shift.
meda : The map of how the model relates the variables to each other.
References
----------
Camacho, J. (2011). Observation-based missing data methods for exploratory data
analysis to unveil the connection between observations and variables in latent
subspace models. Journal of Chemometrics 25(11), 592-600.
"""
X_df, x, R, P = _complete_through_model(model, X)
dummy = _omeda_dummy(X_df, group, reference, weights)
observed = x.T @ dummy
modelled = (x @ R @ P.T).T @ dummy
values = (2.0 * observed - modelled) * np.abs(modelled) / np.sqrt(dummy @ dummy)
return pd.Series(values, index=X_df.columns, name="omeda")
[docs]
def eigenvalue_summary(model: BaseEstimator) -> pd.DataFrame:
"""Summarize the variance captured by each component as a tidy table.
Works with fitted :class:`PCA` and :class:`PLS` models. Returns one row per
component, collecting ``explained_variance_``, ``r2_per_component_`` and
``r2_cumulative_`` into a single table.
Parameters
----------
model : PCA or PLS
A fitted PCA or PLS model.
Returns
-------
pd.DataFrame
Indexed by component, with columns ``eigenvalue`` (the variance of the
component scores), ``percent_variance`` and ``cumulative_percent``. For
PCA the percentages refer to variance in X; for PLS they refer to the
variance in Y explained by each component.
Raises
------
ValueError
If the model is not fitted.
Examples
--------
>>> pca = PCA(n_components=3).fit(X_scaled)
>>> pca.eigenvalue_summary()
>>> eigenvalue_summary(pca)
"""
if not hasattr(model, "r2_per_component_"):
msg = "Model is not fitted. Call fit() before computing the eigenvalue summary."
raise ValueError(msg)
summary = pd.DataFrame(
{
"eigenvalue": model.explained_variance_,
"percent_variance": model.r2_per_component_.to_numpy() * 100.0,
"cumulative_percent": model.r2_cumulative_.to_numpy() * 100.0,
},
index=model.r2_per_component_.index,
)
summary.index.name = "component"
return summary
[docs]
def project_variables(model: BaseEstimator, supplementary_data: DataMatrix) -> pd.DataFrame:
"""Project supplementary (passive) variables onto a fitted model.
Works with fitted :class:`PCA` and :class:`PLS` models. Supplementary
variables are extra columns that did not take part in fitting the model but
were measured on the *same observations*. Each supplementary variable is
represented by its correlation with each component's scores, the standard
representation for passive quantitative variables. This is the column-wise
counterpart of ``transform``, which projects supplementary *rows* (new
observations).
Parameters
----------
model : PCA or PLS
A fitted PCA or PLS model.
supplementary_data : array-like of shape (n_samples, n_supplementary)
Passive variables measured on the same observations used to fit the
model. Must have the same number of rows as the training data.
Returns
-------
pd.DataFrame
Correlations of shape (n_supplementary, n_components): the coordinate
of each supplementary variable on each component.
Raises
------
ValueError
If the model is not fitted, or if *supplementary_data* does not have
the same number of rows as the training data.
Examples
--------
>>> pca = PCA(n_components=3).fit(X_scaled)
>>> pca.project_variables(passive_columns)
>>> project_variables(pca, passive_columns)
"""
if not hasattr(model, "scores_"):
msg = "Model is not fitted. Call fit() before projecting variables."
raise ValueError(msg)
scores = model.scores_
if not isinstance(supplementary_data, pd.DataFrame):
supplementary_data = pd.DataFrame(supplementary_data)
if supplementary_data.shape[0] != scores.shape[0]:
msg = (
f"Supplementary data must have {scores.shape[0]} rows (the number of "
f"observations used to fit the model), got {supplementary_data.shape[0]}."
)
raise ValueError(msg)
xs = supplementary_data.to_numpy(dtype=float)
t = scores.to_numpy()
xs_centered = xs - xs.mean(axis=0)
t_centered = t - t.mean(axis=0)
xs_norm = np.sqrt((xs_centered**2).sum(axis=0))
t_norm = np.sqrt((t_centered**2).sum(axis=0))
with np.errstate(invalid="ignore", divide="ignore"):
correlations = (xs_centered.T @ t_centered) / np.outer(xs_norm, t_norm)
correlations[~np.isfinite(correlations)] = 0.0
return pd.DataFrame(correlations, index=supplementary_data.columns, columns=scores.columns)
def _column_centred_array(data: DataMatrix, name: str) -> np.ndarray:
"""Coerce *data* to a 2-D float array with every column mean-centred."""
arr = np.asarray(data, dtype=float)
if arr.ndim == 1:
arr = arr[:, np.newaxis]
if arr.ndim != 2:
msg = f"{name} must be 1- or 2-dimensional, got {arr.ndim} dimensions."
raise ValueError(msg)
return np.asarray(center(arr))
def _matrix_correlation(X: DataMatrix, Y: DataMatrix, *, modified: bool = False) -> float:
"""Shared core of the RV and modified-RV (RV2) coefficients."""
x_centred = _column_centred_array(X, "X")
y_centred = _column_centred_array(Y, "Y")
if x_centred.shape[0] != y_centred.shape[0]:
msg = f"X and Y must have the same number of rows; got {x_centred.shape[0]} and {y_centred.shape[0]}."
raise ValueError(msg)
s_x = x_centred @ x_centred.T
s_y = y_centred @ y_centred.T
if modified:
np.fill_diagonal(s_x, 0.0)
np.fill_diagonal(s_y, 0.0)
denominator = np.sqrt(np.sum(s_x * s_x) * np.sum(s_y * s_y))
if denominator == 0.0:
return np.nan
return float(np.sum(s_x * s_y) / denominator)
[docs]
def rv_coefficient(X: DataMatrix, Y: DataMatrix) -> float:
"""Compute the RV coefficient between two data blocks.
The RV coefficient (Robert and Escoufier, 1976) measures how much common
structure two matrices, measured on the *same observations*, share. It is a
multivariate generalisation of the squared Pearson correlation: it compares
the observation-by-observation configuration matrices :math:`XX^T` and
:math:`YY^T` rather than individual variables.
Parameters
----------
X : array-like of shape (n_samples, n_features_x)
First data block.
Y : array-like of shape (n_samples, n_features_y)
Second data block. Must have the same number of rows as *X*; the
number of columns may differ.
Returns
-------
float
The RV coefficient in the range [0, 1]. A value of 1 means the two
blocks describe the same configuration of observations up to a
rotation and an overall scaling; 0 means no shared structure.
``nan`` is returned if either block has no variance.
Notes
-----
Each column is mean-centred internally, since the RV coefficient is defined
on centred data. The blocks are **not** scaled; scale the columns yourself
(for example with :class:`MCUVScaler`) when the variables have different
units.
For high-dimensional data (many more variables than observations) the RV
coefficient is biased upwards and tends towards 1 even for unrelated
blocks. Use :func:`rv2_coefficient` in that regime.
References
----------
Robert, P. and Escoufier, Y. (1976). A unifying tool for linear
multivariate statistical methods: the RV-coefficient. *Journal of the
Royal Statistical Society, Series C*, 25(3), 257-265.
See Also
--------
rv2_coefficient : Modified RV coefficient, unbiased for high-dimensional data.
Examples
--------
>>> rv_coefficient(X, Y)
>>> rv_coefficient(X, X) # 1.0: a block is perfectly correlated with itself
"""
return _matrix_correlation(X, Y)
[docs]
def rv2_coefficient(X: DataMatrix, Y: DataMatrix) -> float:
"""Compute the modified RV coefficient (RV2) between two data blocks.
The modified RV coefficient (Smilde et al., 2009) is a variant of
:func:`rv_coefficient` that removes the diagonals of the configuration
matrices :math:`XX^T` and :math:`YY^T` before comparing them. This removes
the upward bias that makes the ordinary RV coefficient tend towards 1 for
high-dimensional data, so RV2 stays near 0 for genuinely unrelated blocks.
Parameters
----------
X : array-like of shape (n_samples, n_features_x)
First data block.
Y : array-like of shape (n_samples, n_features_y)
Second data block. Must have the same number of rows as *X*; the
number of columns may differ.
Returns
-------
float
The modified RV coefficient, in the range [-1, 1]. A value of 1 means
the two blocks describe the same configuration of observations; values
near 0 mean no shared structure, and small negative values can occur.
``nan`` is returned if either block has no variance.
Notes
-----
Each column is mean-centred internally but the blocks are **not** scaled;
scale the columns yourself (for example with :class:`MCUVScaler`) when the
variables have different units.
References
----------
Smilde, A. K., Kiers, H. A. L., Bijlsma, S., Rubingh, C. M. and van Erk,
M. J. (2009). Matrix correlations for high-dimensional data: the modified
RV-coefficient. *Bioinformatics*, 25(3), 401-405.
See Also
--------
rv_coefficient : The original RV coefficient.
Examples
--------
>>> rv2_coefficient(X, Y)
"""
return _matrix_correlation(X, Y, modified=True)
