# (c) Kevin Dunn, 2010-2026. MIT License. Based on own private work over the years.
"""Multi-block PCA (MBPCA) estimator (ENG-01).
Holds :class:`MBPCA`, the hierarchical / superblock multi-block PCA transformer.
"""
from __future__ import annotations
import dataclasses
import itertools
import logging
import numbers
import time
import typing
import warnings
from collections.abc import Sequence
import numpy as np
import pandas as pd
from sklearn.base import BaseEstimator, TransformerMixin, _fit_context
from sklearn.utils import Bunch
from sklearn.utils.validation import check_is_fitted
from ._base import _HotellingsT2LimitMixin
from ._common import (
DiagnoseResult,
SpecificationWarning,
_nz,
_own_axes,
_rows_with_data,
_scale_block_contributions,
check_at_least_one,
check_x_varies,
components_to_fit,
constant_columns,
epsqrt,
informative_widths,
refuse_empty_block,
resolve_loop_settings,
text_table,
unobserved_columns,
user_stacklevel,
)
from ._diagnostics import _select_rows, _t2_score_sd, cumulative_block_t2
from ._limits import spe_calculation
from ._nipals import quick_regress, row_regression_divisors, ssq
from ._preprocessing import MCUVScaler
try:
import plotly.graph_objects as go
except ImportError: # pragma: no cover - exercised via env-without-plotly
from process_improve._extras import _MissingExtra
go = _MissingExtra("plotly", "plotting") # type: ignore[assignment]
logger = logging.getLogger(__name__)
def _validate_blocks(X: dict[str, pd.DataFrame]) -> None:
"""Reject anything that is not a non-empty dict of equally tall DataFrames."""
if not isinstance(X, dict) or len(X) == 0:
raise TypeError("X must be a non-empty dict[str, pd.DataFrame].")
for name, block in X.items():
if not isinstance(block, pd.DataFrame):
raise TypeError(f"X['{name}'] must be a pandas DataFrame; got {type(block).__name__}.")
refuse_empty_block(block, f'X["{name}"]')
names = list(X)
n_samples = X[names[0]].shape[0]
for name in names:
if X[name].shape[0] != n_samples:
raise ValueError(
f"All X-blocks must have the same row count. Block '{name}' has "
f"{X[name].shape[0]} rows; expected {n_samples}."
)
class _MBPCALoopContext(typing.NamedTuple):
"""Everything the component loop needs that does not change between components."""
algo: str
block_names: list[str]
#: Algorithmic block weighting: X_b / sqrt(K_b), so blocks of unequal width contribute fairly.
sqrt_kb: dict[str, float]
block_has_data: dict[str, np.ndarray]
tol: float
max_iter: int
n_samples: int
@property
def n_blocks(self) -> int:
"""How many X-blocks this fit is over."""
return len(self.block_names)
class _MBPCAComponent(typing.NamedTuple):
"""One converged component, before deflation."""
super_scores: np.ndarray
super_loadings: np.ndarray
block_loadings: dict[str, np.ndarray]
block_scores: dict[str, np.ndarray]
iterations: int
converged: bool
def _seed_super_score(x_def: dict[str, np.ndarray], context: _MBPCALoopContext) -> np.ndarray:
"""Deterministic start (#503): the column, across all blocks, with the largest sum of squares.
Mirrors the single-block PCA / PLS seeding (#195). No RNG is involved, so the fit is
reproducible without a ``random_state``, and the highest-variance column is closest to
the leading component. The sign convention applied after convergence makes the fitted
signs independent of this seed. NaN is replaced by 0 for the missing-data path, in a
copy, so the seed is not a view of the block.
"""
col_ssq = {name: np.nansum(x_def[name] ** 2, axis=0) for name in context.block_names}
start_block = max(context.block_names, key=lambda name: float(np.max(col_ssq[name])))
start_col = int(np.argmax(col_ssq[start_block]))
return np.nan_to_num(x_def[start_block][:, start_col])
def _fit_one_component(x_def: dict[str, np.ndarray], context: _MBPCALoopContext) -> _MBPCAComponent:
"""Iterate the super-score to convergence on the currently deflated blocks.
The first pass always runs. A centred seed changes by about its own size in that pass,
so a tolerance above 1 used to end the loop before it began and leave no loadings.
"""
t_super = _seed_super_score(x_def, context)
t_b_summary = np.zeros((context.n_samples, context.n_blocks))
local_loadings: dict[str, np.ndarray] = {}
local_scores: dict[str, np.ndarray] = {}
# On the dense path each loading is normalised as soon as it is formed, so dividing it by
# the score's sum of squares first would change nothing; nor would dividing a score by
# the sum of squares of a loading of unit length.
for itern in itertools.count(1):
prev = t_super
if context.algo == "nipals":
# Mask-aware NIPALS: each projection is a per-column (or per-row) regression
# that uses only the entries that are not NaN, and divides by the masked sum of
# squares. Reuses the same primitives as single-block PCA NIPALS.
t_super_col = t_super[:, np.newaxis]
for b_idx, name in enumerate(context.block_names):
p_b = quick_regress(x_def[name], t_super_col, per="column").flatten()
p_b = p_b / _nz(float(np.sqrt(ssq(p_b))))
t_b = quick_regress(x_def[name], p_b[:, np.newaxis], per="row").flatten() / context.sqrt_kb[name]
local_loadings[name] = p_b
local_scores[name] = t_b
t_b_summary[:, b_idx] = np.where(
context.block_has_data[name],
t_b,
np.nan, # pragma: no mutate - None reads as NaN in a float array
)
# Masked, as MBPLS's super weight is: a block a row has nothing in has no say in the
# super loading, nor in that row's super score, which comes from the blocks it has.
# A row with nothing observed anywhere scores 0 and has no weight, as in PCA (#1098).
p_s = quick_regress(t_b_summary, t_super_col, per="column").flatten()
p_s = p_s / _nz(float(np.linalg.norm(p_s)))
t_super = quick_regress(t_b_summary, p_s[:, np.newaxis], per="row").flatten()
else:
for b_idx, name in enumerate(context.block_names):
p_b = x_def[name].T @ t_super
p_b = p_b / _nz(float(np.linalg.norm(p_b)))
t_b = x_def[name] @ p_b / context.sqrt_kb[name]
local_loadings[name] = p_b
local_scores[name] = t_b
t_b_summary[:, b_idx] = t_b
p_s = t_b_summary.T @ t_super
p_s = p_s / _nz(float(np.linalg.norm(p_s)))
t_super = t_b_summary @ p_s
# Relative convergence criterion (#504): the change between two successive super-score
# iterations is judged against the size of the current super-score, so the decision is
# invariant to a global rescaling of the data. The denominator is floored via ``_nz``
# so an all-zero super-score cannot divide by zero.
change = np.linalg.norm(prev - t_super) / _nz(float(np.linalg.norm(t_super)))
# Converged means the last pass met the tolerance, even when it was the last pass allowed (#1080).
converged = bool(change <= context.tol)
if converged or itern >= context.max_iter:
break
component = _MBPCAComponent(
super_scores=t_super,
super_loadings=p_s,
block_loadings=local_loadings,
block_scores=local_scores,
iterations=itern,
converged=converged,
)
return _apply_sign_convention(component, context)
def _deflation_loadings(component: _MBPCAComponent, context: _MBPCALoopContext) -> dict[str, np.ndarray]:
"""Scale each block loading by its super loading and block size: the loading that deflates it."""
return {
name: component.block_loadings[name] * component.super_loadings[b_idx] * context.sqrt_kb[name]
for b_idx, name in enumerate(context.block_names)
}
def _apply_sign_convention(component: _MBPCAComponent, context: _MBPCALoopContext) -> _MBPCAComponent:
"""Flip the component so its largest-magnitude deflation loading, over all blocks, is positive.
This is single-block PCA's convention, so the fitted signs do not depend on the seed,
and a one-block model has PCA's signs (#586). The super loading is not flipped, and
cannot decide the sign: its element for block b is proportional to ``t_b' t``, which
flipping the whole component leaves unchanged. (With no missing data that is
``||X_b' t||``, never negative, so the old convention keyed on it never fired.)
"""
loadings = _deflation_loadings(component, context)
stacked = np.nan_to_num(np.concatenate([loadings[name] for name in context.block_names]))
if stacked[int(np.argmax(np.abs(stacked)))] >= 0:
return component
return component._replace(
super_scores=-component.super_scores,
block_loadings={name: -loading for name, loading in component.block_loadings.items()},
block_scores={name: -scores for name, scores in component.block_scores.items()},
)
def _deflate(
x_def: dict[str, np.ndarray], component: _MBPCAComponent, context: _MBPCALoopContext
) -> dict[str, np.ndarray]:
"""Remove this component from every block, using the super-score and the scaled block loading."""
loadings = _deflation_loadings(component, context)
return {name: x_def[name] - np.outer(component.super_scores, loadings[name]) for name in context.block_names}
@dataclasses.dataclass
class _MBPCAArrays:
"""The numpy workspace a fit fills in, one column per component, wrapped in pandas at the end."""
super_scores: np.ndarray
super_loadings: np.ndarray
block_scores: dict[str, np.ndarray]
block_loadings: dict[str, np.ndarray]
block_spe: dict[str, np.ndarray]
r2_x_block_cum: np.ndarray
r2_x_var_cum: dict[str, np.ndarray]
timing: np.ndarray
iterations: np.ndarray
converged: np.ndarray
@classmethod
def allocate(cls, block_widths: dict[str, int], n_samples: int, n_components: int) -> _MBPCAArrays:
"""Allocate the workspace for a fit of the given shape."""
n_blocks = len(block_widths)
per_block = {name: np.zeros((n_samples, n_components)) for name in block_widths}
per_variable = {name: np.zeros((width, n_components)) for name, width in block_widths.items()}
return cls(
super_scores=np.zeros((n_samples, n_components)),
super_loadings=np.zeros((n_blocks, n_components)),
block_scores={name: values.copy() for name, values in per_block.items()},
block_loadings={name: values.copy() for name, values in per_variable.items()},
block_spe={name: values.copy() for name, values in per_block.items()},
r2_x_block_cum=np.zeros((n_blocks, n_components)),
r2_x_var_cum={name: values.copy() for name, values in per_variable.items()},
timing=np.zeros(n_components),
iterations=np.zeros(n_components, dtype=int),
converged=np.zeros(n_components, dtype=bool),
)
def record_component(self, a: int, component: _MBPCAComponent, context: _MBPCALoopContext) -> None:
"""Store one converged component's scores and loadings.
A block a row has nothing observed in gets NaN for its block score, not a zero, which
would place the row at that block's average: MBPLS's rule (#1098).
"""
for b_idx, name in enumerate(context.block_names):
self.block_loadings[name][:, a] = component.block_loadings[name]
self.block_scores[name][:, a] = np.where(
context.block_has_data[name],
component.block_scores[name],
np.nan, # pragma: no mutate - None reads as NaN in a float array
)
self.super_loadings[b_idx, a] = component.super_loadings[b_idx]
self.super_scores[:, a] = component.super_scores
self.iterations[a] = component.iterations
self.converged[a] = component.converged
def record_explained_variation(
self,
a: int,
x_def: dict[str, np.ndarray],
initial: _InitialSumsOfSquares,
context: _MBPCALoopContext,
) -> None:
"""Store the cumulative R2X and SPE left after this component was removed."""
for b_idx, name in enumerate(context.block_names):
ssq_remain_per_var = np.nansum(x_def[name] ** 2, axis=0)
self.r2_x_block_cum[b_idx, a] = 1 - np.sum(ssq_remain_per_var) / initial.per_block[name]
self.r2_x_var_cum[name][:, a] = 1 - ssq_remain_per_var / initial.per_variable[name]
# A row with nothing observed in the block has no SPE there, as in MBPLS (#1028, #1098)
self.block_spe[name][:, a] = np.where(
context.block_has_data[name],
np.sqrt(np.nansum(x_def[name] ** 2, axis=1)),
np.nan, # pragma: no mutate - None reads as NaN in a float array
)
class _InitialSumsOfSquares(typing.NamedTuple):
"""What each block held before any component was removed, the reference every R2X is against.
R2X is undefined for a block or a column with nothing to explain, so its sum of squares is
held as NaN: every R2X divided by it is NaN, without the warning a division by 0 gives.
"""
per_block: dict[str, float]
per_variable: dict[str, np.ndarray]
@classmethod
def of(cls, x_blocks_pp: dict[str, np.ndarray]) -> _InitialSumsOfSquares:
"""Measure the preprocessed blocks before the first deflation."""
per_variable = {name: np.nansum(values**2, axis=0) for name, values in x_blocks_pp.items()}
per_block = {name: float(np.sum(values)) for name, values in per_variable.items()}
return cls(
per_block={name: total if total > 0 else np.nan for name, total in per_block.items()},
per_variable={name: np.where(values > 0, values, np.nan) for name, values in per_variable.items()},
)
[docs]
class MBPCA(_HotellingsT2LimitMixin, TransformerMixin, BaseEstimator):
r"""Multi-block PCA (hierarchical / consensus PCA).
Generic multi-block PCA following the consensus-PCA / hierarchical PCA
formulation of Westerhuis, Kourti & MacGregor (1998). Each X-block is
preprocessed independently (mean-centred and unit-variance scaled),
then divided by ``sqrt(K_b)`` so blocks of unequal width contribute
fairly to the consensus super-score.
The hierarchical NIPALS loop alternates: (i) regress each block on the
super-score to get block loadings and block scores, (ii) collect block
scores into a super-block, (iii) regress the super-block to get a new
super-score / super-loading, repeat to convergence. After convergence,
deflate every block by the super-score and the corresponding block
loading scaled by the super-loading element.
Parameters
----------
n_components : int
Number of super-components (consensus latent variables) to extract.
max_iter : int, default=500
Maximum NIPALS iterations per component in the hierarchical outer loop.
tol : float, default=``epsqrt`` (about 1.5e-8)
Relative convergence tolerance on the super-score change: the norm of
the change between two successive super-score iterations, divided by
the norm of the current super-score vector (#504). The same default,
and the same meaning, as on PCA / PLS / TPLS. It used to default to
``None`` and be resolved to ``epsqrt`` inside ``fit``; the literal is
the default now, which removes a branch and lets ``get_params`` report
the value that will actually be used (#588). The legacy absolute
tolerance ``np.finfo(float).eps ** (9/10)`` (about 8.2e-15) sits below
the floating-point oscillation floor of a relative criterion, so it
would never be reached in practice.
algorithm : str, default="auto"
Algorithm to use for fitting the model.
- ``"auto"``: dense vectorised hierarchical NIPALS when the data is
complete; mask-aware NIPALS (NaN-tolerant) when any block contains
missing values.
- ``"dense"``: dense vectorised hierarchical NIPALS. Raises if any
block contains missing values.
- ``"nipals"``: mask-aware hierarchical NIPALS. Always uses the
NaN-tolerant inner-loop primitives, even when the data is
complete (slower than ``"dense"`` but produces equivalent
results).
missing_data_settings : dict or None, default=None
MBPCA reads no keys from it. The keys ``md_tol`` and ``md_max_iter``
were removed in 2.0 and raise a :class:`ValueError`; set ``tol`` and
``max_iter`` instead.
Attributes (after fitting)
--------------------------
n_components_ : int
The number of components fitted. Everything after ``fit`` reads it, not the
``n_components`` parameter, so ``set_params`` without a refit changes nothing.
block_names_ : list[str]
Ordered list of X-block names (the keys of the input dict).
block_widths_ : dict[str, int]
Number of variables in each X-block.
block_informative_widths_ : dict[str, int]
Number of variables in each block whose values vary (at least 1). The block
weighting divides by its square root, so constant columns do not lower a
block's weight.
n_samples_ : int
Number of rows fitted.
n_features_in_ : int
Total number of X variables summed across blocks.
feature_names_in_ : np.ndarray
Concatenated column names, one per feature, in block order.
preproc_ : dict[str, MCUVScaler]
Per-block preprocessors used to mean-centre and unit-variance
scale each X-block.
super_scores_ : pd.DataFrame, shape (n_samples, n_components)
Super-block (consensus) scores ``T``.
super_loadings_ : pd.DataFrame, shape (n_blocks, n_components)
Super-block loadings ``p_super``; rows indexed by block name.
super_hotellings_t2_ : pd.DataFrame, shape (n_samples, n_components)
Cumulative Hotelling's T^2 on the super-scores per component.
block_scores_ : dict[str, pd.DataFrame]
Per-block scores ``t_b``, each shape ``(n_samples, n_components)``. NaN for a row with
nothing observed in that block, as in :class:`MBPLS`.
block_loadings_ : dict[str, pd.DataFrame]
Per-block loadings ``p_b``, each shape ``(K_b, n_components)``.
block_spe_ : dict[str, pd.DataFrame]
Per-block squared prediction error per sample and component. NaN for a row with
nothing observed in that block.
block_hotellings_t2_ : dict[str, pd.DataFrame]
Per-block cumulative Hotelling's T^2 per sample and component. Each component is
scaled by its block scores' second moment about zero, as ``super_hotellings_t2_`` is
by ``explained_variance_``.
block_vip_ : dict[str, pd.Series]
Per-block variable-importance in projection, indexed by variable
name inside each block. Computed from ``block_loadings_``, each component
weighted by its share of that block's X variance
(``r2_x_per_block_per_component_``), as a PCA VIP is.
r2_x_per_block_cumulative_ : pd.DataFrame, shape (n_blocks, n_components)
Cumulative R^2X per block and component.
r2_x_per_block_per_component_ : pd.DataFrame, shape (n_blocks, n_components)
Incremental R^2X per block and component.
r2_x_per_variable_ : dict[str, pd.DataFrame]
Cumulative R^2X per variable within each block.
explained_variance_ : np.ndarray, shape (n_components,)
Variance of the super-score per component (ddof=1).
scaling_factor_for_super_scores_ : pd.Series
``sqrt(explained_variance_)`` per component.
fitting_info_ : dict
Per-component iteration count and timing.
has_missing_data_ : bool
Whether any X-block had NaN values.
algorithm_ : str
The resolved algorithm actually used for the fit. With
``algorithm="auto"``, this is ``"dense"`` for complete data
and ``"nipals"`` for NaN-containing data.
Notes
-----
The deflation step is :math:`X_b \leftarrow X_b - t_{\rm super}\,
(p_b\,p_s[b]\,\sqrt{K_b})^\top`, derived in Westerhuis et al. 1998.
An earlier implementation of this method had this step marked as broken
by its author; this implementation re-derives it directly from the paper
and is independently validated against the pure-numpy reference oracles
in the test suite.
Missing data
------------
When any block contains NaN entries, the ``"auto"`` algorithm
routes to a mask-aware NIPALS variant. Each per-block projection
in the inner loop is computed as a regression that uses only the
observed entries; the masked sum-of-squares is used as the
denominator so missing values neither bias the loading direction
nor contribute to the score. The mask is preserved across
components automatically because deflation propagates NaN through
subtraction. This is the standard skip-NaN NIPALS update; see
Walczak & Massart (2001) and Arteaga & Ferrer (2002).
A row with nothing observed in a block has no score for that block, and
no say in the super loading; its super score comes from the blocks it
does have. A row with nothing observed in any block is fitted with a
super score of 0 and no weight, as :class:`PCA`, :class:`PLS` and
:class:`TPLS` fit an all-missing row. A column with nothing observed
is left out of the model, with a ``UserWarning`` that names it: its
loadings are 0, and a value a new row has in it adds nothing to that
row's scores or SPE (#1098).
:meth:`transform` and :meth:`diagnose` score a new row as the fit
scores its own: over the cells the row has, so a training row is given
its fitted scores again (#1076). Its block SPE is over those cells too,
and NaN for a block it has nothing observed in. A new row with nothing
observed at all has NaN scores, SPE and T², where the fit places such a
row at 0.
References
----------
Westerhuis, J. A., Kourti, T. & MacGregor, J. F. *Analysis of
multiblock and hierarchical PCA and PLS models.* J. Chemometrics, 12
(1998), 301-321.
Walczak, B. & Massart, D. L. *Dealing with missing data: Part I.*
Chemom. Intell. Lab. Syst., 58 (2001), 15-27.
Arteaga, F. & Ferrer, A. *Dealing with missing data in MSPC: several
methods, different interpretations, some examples.* J. Chemometrics,
16 (2002), 408-418.
"""
_valid_algorithms: typing.ClassVar[list[str]] = ["auto", "dense", "nipals"]
_parameter_constraints: typing.ClassVar = {
"n_components": [numbers.Integral],
"max_iter": [numbers.Integral],
"tol": [numbers.Real],
"algorithm": [str],
"missing_data_settings": [dict, None],
}
def __init__(
self,
n_components: int,
*,
max_iter: int = 500,
tol: float = epsqrt,
algorithm: str = "auto", # pragma: no mutate - _resolve_algorithm lowercases it, so "AUTO" is the same default
missing_data_settings: dict | None = None,
):
super().__init__()
self.n_components = n_components
self.max_iter = max_iter
self.tol = tol
self.algorithm = algorithm
self.missing_data_settings = missing_data_settings
[docs]
@_fit_context(prefer_skip_nested_validation=True)
def fit(self, X: dict[str, pd.DataFrame], y: None = None) -> MBPCA: # noqa: ARG002
"""Fit the multi-block PCA model.
Parameters
----------
X : dict[str, pd.DataFrame]
X-blocks. Keys are block names; values are DataFrames sharing the same row
index (and row count). Each block is preprocessed independently.
y : None
Ignored; accepted so the transformer plugs into a sklearn Pipeline.
Raises
------
NotEnoughVarianceError
If every column of every block is constant or has no observed value, as across the
family (#1144). A block with no variance beside one that varies is fitted with zero
weight: its block scores, loadings and contributions are 0, and a warning names it.
"""
check_at_least_one("n_components", self.n_components)
_validate_blocks(X)
X = {name: _own_axes(block) for name, block in X.items()} # nothing the caller renames later reaches the fit
self._record_data_shape(X)
# A column with nothing observed is left out, with a warning, as across the family (#1098)
self._unobserved_columns = {name: unobserved_columns(X[name], f'X["{name}"]') for name in self.block_names_}
algo = self._resolve_algorithm(X)
# The resolved settings reach the loop. They used to be computed here and dropped,
# so the caller's loop settings were validated and then ignored (#588).
settings = resolve_loop_settings(
tol=self.tol,
max_iter=self.max_iter,
missing_data_settings=self.missing_data_settings,
estimator_name="MBPCA",
)
x_blocks_pp = self._preprocess(X) # refuses a block of fewer than two rows
# An X in which nothing varies tells no row apart: refused, as across the family, before its columns are named
check_x_varies([X[name] for name in self.block_names_], explains_y=False)
for name in self.block_names_:
constant_columns(X[name], f'X["{name}"]', self)
self.block_informative_widths_: dict[str, int] = informative_widths(x_blocks_pp)
context = _MBPCALoopContext(
algo=algo,
block_names=self.block_names_,
# Block weighting: X_b / sqrt(K_b), K_b counting the columns that vary (#913)
sqrt_kb=self._block_scale(),
# A row with nothing observed in a block has no score for that block (#1098), as in MBPLS
block_has_data={name: _rows_with_data(values) for name, values in x_blocks_pp.items()},
tol=float(settings["md_tol"]),
max_iter=int(settings["md_max_iter"]),
n_samples=self.n_samples_,
)
work = _MBPCAArrays.allocate(self.block_widths_, self.n_samples_, self.n_components_)
initial = _InitialSumsOfSquares.of(x_blocks_pp)
x_def = {name: values.copy() for name, values in x_blocks_pp.items()}
for a in range(self.n_components_):
start = time.time()
component = _fit_one_component(x_def, context)
x_def = _deflate(x_def, component, context)
work.record_component(a, component, context)
work.record_explained_variation(a, x_def, initial, context)
work.timing[a] = time.time() - start
self._store_latent_frames(work)
self._report_convergence(work)
self._store_explained_variation(work)
self._store_diagnostics(work)
return self
def _block_scale(self) -> dict[str, float]:
"""Return what each block is divided by: the square root of its count of varying columns (#913).
A model fitted before that count was kept divides by its full width, as it was fitted.
"""
widths = getattr(self, "block_informative_widths_", self.block_widths_)
return {name: float(np.sqrt(widths[name])) for name in self.block_names_}
def _record_data_shape(self, X: dict[str, pd.DataFrame]) -> None:
"""Record the block names, widths, row and column labels, and counts this fit is over."""
self.block_names_: list[str] = list(X.keys())
first = X[self.block_names_[0]]
self.block_widths_: dict[str, int] = {name: int(X[name].shape[1]) for name in self.block_names_}
self._sample_index = first.index
self._block_columns: dict[str, pd.Index] = {name: X[name].columns for name in self.block_names_}
self.n_samples_ = int(first.shape[0])
self.n_features_in_ = int(sum(self.block_widths_.values()))
# feature_names_in_: sklearn convention (#392). Flat concatenation of all blocks'
# column names in block-iteration order. Lets ``Pipeline.get_feature_names_out``
# and SHAP / eli5 / model-card tooling introspect a multiblock fit through the
# same surface as a single-block estimator.
self.feature_names_in_ = np.concatenate([self._block_columns[name].to_numpy() for name in self.block_names_])
# Fitted mirror of the constructor parameter, so shared helpers (the T2 limit
# mixin, spe_limit, the plot pre-checks) read one resolved attribute across
# PCA / PLS / MBPCA / MBPLS (#505).
self.n_components_ = components_to_fit(self.n_components, self.n_samples_, self.n_features_in_)
def _resolve_algorithm(self, X: dict[str, pd.DataFrame]) -> str:
"""Pick between the dense and NIPALS paths, and record what was picked."""
self.has_missing_data_ = any(np.any(X[name].isna().values) for name in self.block_names_)
algo = self.algorithm.lower()
if algo not in self._valid_algorithms:
raise ValueError(
f"Algorithm '{self.algorithm}' is not recognised. Must be one of {self._valid_algorithms}."
)
if algo == "auto":
algo = "nipals" if self.has_missing_data_ else "dense"
if algo == "dense" and self.has_missing_data_:
raise ValueError("Algorithm 'dense' cannot handle missing data. Use 'nipals' or 'auto' instead.")
self.algorithm_ = algo
return algo
def _preprocess(self, X: dict[str, pd.DataFrame]) -> dict[str, np.ndarray]:
"""Mean-centre and unit-variance scale each block independently."""
self.preproc_: dict[str, MCUVScaler] = {name: MCUVScaler().fit(X[name]) for name in self.block_names_}
return {name: self.preproc_[name].transform(X[name]).to_numpy() for name in self.block_names_}
def _component_names(self) -> list[int]:
"""Column labels for every per-component frame: 1-based component numbers."""
return list(range(1, self.n_components_ + 1))
def _store_latent_frames(self, work: _MBPCAArrays) -> None:
"""Wrap the scores and loadings in labelled pandas containers."""
component_names = self._component_names()
self.super_scores_ = pd.DataFrame(work.super_scores, index=self._sample_index, columns=component_names)
self.super_loadings_ = pd.DataFrame(work.super_loadings, index=self.block_names_, columns=component_names)
self.block_scores_ = {
name: pd.DataFrame(work.block_scores[name], index=self._sample_index, columns=component_names)
for name in self.block_names_
}
self.block_loadings_ = {
name: pd.DataFrame(work.block_loadings[name], index=self._block_columns[name], columns=component_names)
for name in self.block_names_
}
# fit() has at least two rows: each block's MCUVScaler refuses fewer
self.explained_variance_ = np.diag(work.super_scores.T @ work.super_scores) / (self.n_samples_ - 1)
self.scaling_factor_for_super_scores_ = pd.Series(
np.sqrt(self.explained_variance_), index=component_names, name="Standard deviation per super-score"
)
def _report_convergence(self, work: _MBPCAArrays) -> None:
"""Record the per-component timing and iteration counts, and warn about any that did not converge."""
converged = work.converged
self.fitting_info_ = {"timing": work.timing, "iterations": work.iterations, "converged": converged}
logger.debug("MBPCA (%s): iterations per component = %s", self.algorithm_, work.iterations.tolist())
if np.all(converged):
return
failed = [int(i + 1) for i, ok in enumerate(converged) if not ok]
warnings.warn(
f"MBPCA NIPALS did not converge within max_iter={self.max_iter} for "
f"component(s) {failed}; results for those components may be unreliable.",
SpecificationWarning,
stacklevel=user_stacklevel(),
)
def _store_explained_variation(self, work: _MBPCAArrays) -> None:
"""Turn the cumulative R2X into the cumulative, per-component and per-variable frames."""
component_names = self._component_names()
r2_x_block_per_a = np.diff(work.r2_x_block_cum, prepend=0.0)
self.r2_x_per_block_cumulative_ = pd.DataFrame(
work.r2_x_block_cum, index=self.block_names_, columns=component_names
)
self.r2_x_per_block_per_component_ = pd.DataFrame(
r2_x_block_per_a, index=self.block_names_, columns=component_names
)
self.r2_x_per_variable_ = {
name: pd.DataFrame(work.r2_x_var_cum[name], index=self._block_columns[name], columns=component_names)
for name in self.block_names_
}
# VIPs (per-block: variance-of-X explanation; the super VIP is the same idea on
# the super-loadings).
self.block_vip_: dict[str, pd.Series] = {}
for b_idx, name in enumerate(self.block_names_):
r2 = r2_x_block_per_a[b_idx, :]
if np.sum(r2) > 0:
p = self.block_loadings_[name].values
vip_b = np.sqrt(self.block_informative_widths_[name] * np.sum(r2 * p**2, axis=1) / np.sum(r2))
else:
vip_b = np.zeros(self.block_widths_[name])
self.block_vip_[name] = pd.Series(vip_b, index=self._block_columns[name], name=f"VIP[{name}]")
def _store_diagnostics(self, work: _MBPCAArrays) -> None:
"""Build the per-block SPE and Hotelling's T2 frames, and the super-block T2."""
component_names = self._component_names()
self.block_spe_ = {
name: pd.DataFrame(work.block_spe[name], index=self._sample_index, columns=component_names)
for name in self.block_names_
}
# Scaled by the second moment about zero, as the super T2 is by explained_variance_
self.block_hotellings_t2_ = {
name: pd.DataFrame(
cumulative_block_t2(self.block_scores_[name].values), index=self._sample_index, columns=component_names
)
for name in self.block_names_
}
# A super score with no spread carries no T², by the rule PCA and PLS use (#1026)
super_score_var = _t2_score_sd(np.sqrt(self.explained_variance_)) ** 2
super_t2 = np.cumsum((work.super_scores**2) / super_score_var, axis=1)
self.super_hotellings_t2_ = pd.DataFrame(super_t2, index=self._sample_index, columns=component_names)
[docs]
def diagnose(self, X: dict[str, pd.DataFrame]) -> Bunch:
"""Project new data; return super_scores, block_scores, block_spe, cumulative_hotellings_t2.
``cumulative_hotellings_t2`` is a DataFrame, n_samples x n_components: column ``a``
holds the T² over the first ``a`` components, as in ``super_hotellings_t2_``; every
``diagnose`` of the family returns it in this shape (#983). ``hotellings_t2``, the
Series of the full model's T², is deprecated since 2.4.0 and is removed in 3.0; it is
``cumulative_hotellings_t2.iloc[:, -1]``, and reading it warns.
A row with missing cells is scored over the cells it has, as the fit
scores a row. A row with nothing observed in a block has NaN for its
block scores and its block SPE there, as in ``block_scores_`` and
``block_spe_``; a row with nothing observed at all has NaN scores and T²
as well.
The rename (since 1.38.4, #395) matches :meth:`PCA.diagnose` and
:meth:`PLS.diagnose`. The deprecated ``predict`` alias was removed in
2.0.
"""
check_is_fitted(self, "super_loadings_")
return self._project(X)
def _scaled_block(self, name: str, block: pd.DataFrame) -> np.ndarray:
"""Centre and scale one block of new rows as the fit did.
A column the fit left out, having no observed value (#1098), is read as missing, as every
cell of it was in the fit: a value there adds nothing to a score or an SPE, and a row has
nothing observed in the block unless it has a value in one of the block's other columns.
"""
# One memory layout whatever pandas hands over, as in MBPLS._scaled_block (#1148)
values = self.preproc_[name].transform(block).to_numpy()
return np.ascontiguousarray(np.where(self._left_out(name), np.nan, values))
def _left_out(self, name: str) -> np.ndarray | bool:
"""Which columns of a block the fit left out (#1098); a model pickled before then left none out."""
return getattr(self, "_unobserved_columns", {}).get(
name,
False, # pragma: no mutate - np.where reads a condition of None as False
)
def _project_component(
self, x_def: dict[str, np.ndarray], a: int
) -> tuple[dict[str, np.ndarray], np.ndarray, dict[str, np.ndarray]]:
"""Score new rows on component ``a`` as the fit scores its own rows, then deflate them by it.
Each score is the fit's regression over the cells a row has, so a missing cell leaves the
row scored from the rest (#1076). A block the row has nothing observed in has no block
score and no say in the super score. A row with nothing observed at all has no super score
either: NaN, where the fit places it at 0. On a complete row each regression divides by a
loading of unit length, so the scores are the plain products.
Returns the block scores, the super score, and the blocks with the component removed.
"""
sqrt_kb = self._block_scale()
p_s = self.super_loadings_.values[:, a]
block_scores: dict[str, np.ndarray] = {}
for name in self.block_names_:
t_b = quick_regress(x_def[name], self.block_loadings_[name].values[:, [a]], per="row").flatten()
block_scores[name] = np.where(
_rows_with_data(x_def[name]),
t_b / sqrt_kb[name],
np.nan, # pragma: no mutate - None reads as NaN in a float array
)
summary = np.column_stack([block_scores[name] for name in self.block_names_])
t_super = np.where(
_rows_with_data(summary),
quick_regress(summary, p_s[:, np.newaxis], per="row").flatten(),
np.nan, # pragma: no mutate - None reads as NaN in a float array
)
# A row with no super score has nothing observed, so the NaN removes nothing from it.
deflated = {
name: x_def[name] - np.outer(t_super, self.block_loadings_[name].values[:, a] * p_s[b_idx] * sqrt_kb[name])
for b_idx, name in enumerate(self.block_names_)
}
return block_scores, t_super, deflated
def _project(self, X: dict[str, pd.DataFrame]) -> Bunch:
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 for prediction: {sorted(missing)}.")
x_def: dict[str, np.ndarray] = {}
sample_index = None
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 block.shape[1] != self.block_widths_[name]:
raise ValueError(f"Block '{name}' must have {self.block_widths_[name]} columns; got {block.shape[1]}.")
x_def[name] = self._scaled_block(name, block)
if sample_index is None:
sample_index = block.index
n_components = self.n_components_
n_new = next(iter(x_def.values())).shape[0]
super_scores = np.zeros((n_new, n_components))
block_scores: dict[str, np.ndarray] = {name: np.zeros((n_new, n_components)) for name in self.block_names_}
for a in range(n_components):
t_b, super_scores[:, a], x_def = self._project_component(x_def, a)
for name in self.block_names_:
block_scores[name][:, a] = t_b[name]
component_names = list(range(1, n_components + 1))
super_scores_df = pd.DataFrame(super_scores, index=sample_index, columns=component_names)
block_scores_df = {
name: pd.DataFrame(block_scores[name], index=sample_index, columns=component_names)
for name in self.block_names_
}
# Over the cells the row has. A block with none has no SPE, NaN, not the 0 of a perfect fit (#1076).
block_spe = {
name: pd.Series(
np.where(
_rows_with_data(x_def[name]),
np.sqrt(np.nansum(x_def[name] ** 2, axis=1)),
np.nan, # pragma: no mutate - None reads as NaN in a float array
),
index=sample_index,
name=f"SPE[{name}]",
)
for name in self.block_names_
}
# A super score with no spread carries no T², by the rule PCA and PLS use (#1026)
super_score_var = _t2_score_sd(np.sqrt(self.explained_variance_)) ** 2
standardised = (super_scores**2) / super_score_var
cumulative = pd.DataFrame(np.cumsum(standardised, axis=1), index=sample_index, columns=component_names)
hotellings_t2 = pd.Series(np.sum(standardised, axis=1), index=sample_index, name="Hotelling's T²")
return DiagnoseResult(
f"multivariate.{type(self).__name__}",
super_scores=super_scores_df,
block_scores=block_scores_df,
block_spe=block_spe,
cumulative_hotellings_t2=cumulative,
hotellings_t2=hotellings_t2,
)
[docs]
def block_spe_limit(self, block: str, conf_level: float = 0.95) -> float:
"""SPE limit for one X-block (Nomikos & MacGregor chi-square approximation)."""
check_is_fitted(self, "block_spe_")
if block not in self.block_spe_:
raise KeyError(f"Unknown block '{block}'. Known blocks: {list(self.block_spe_)}.")
# A row with nothing observed in the block has no SPE there (NaN), and no part in its limit.
return spe_calculation(self.block_spe_[block].iloc[:, -1].dropna().to_numpy(), conf_level=conf_level)
[docs]
def super_spe_limit(self, conf_level: float = 0.95) -> float:
"""SPE limit for the merged super-block.
A block a row has nothing observed in leaves no residual in that row: its NaN SPE
counts as 0 in the sum, as in :meth:`MBPLS.super_spe_limit`.
"""
check_is_fitted(self, "block_spe_")
spe = np.column_stack([self.block_spe_[name].iloc[:, -1].to_numpy() for name in self.block_names_])
return spe_calculation(np.sqrt(np.nansum(spe**2, axis=1)), conf_level=conf_level)
[docs]
def spe_contributions(self, X: dict[str, pd.DataFrame]) -> dict[str, pd.DataFrame]:
"""Per-variable squared residuals for each X-block (SPE contributions).
Reconstruction matches the MBPCA deflation step:
``X_b = T_super @ (P_b * p_super[b] * sqrt(K_b))^T`` summed over
components. Returns squared residuals on the preprocessed scale; sum
across columns equals ``block_spe_[b].iloc[:, -1] ** 2``.
"""
check_is_fitted(self, "block_loadings_")
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)}.")
result = self._project(X)
super_scores = result.super_scores.values # (N, A)
sample_index = next(iter(result.block_scores.values())).index
sqrt_kb = self._block_scale()
out: dict[str, pd.DataFrame] = {}
for b_idx, name in enumerate(self.block_names_):
block = X[name]
if not isinstance(block, pd.DataFrame):
block = pd.DataFrame(block, columns=self._block_columns[name])
x_pp = self._scaled_block(name, block)
# X_b reconstruction summed over components
p_b = self.block_loadings_[name].values # (K_b, A)
p_s = self.super_loadings_.values[b_idx, :] # (A,)
p_eff = p_b * p_s * sqrt_kb[name] # (K_b, A) effective loading
x_hat = super_scores @ p_eff.T
# A column the fit left out is read as missing, NaN, but contributes 0, as in MBPLS, PCA and PLS (#1114)
residuals_sq = np.where(self._left_out(name), 0.0, (x_pp - x_hat) ** 2)
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 = 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._scaled_block(name, block)
for a in range(int(component) - 1):
x_def = self._project_component(x_def, a)[2]
assert sample_index is not None
return x_def, sample_index
[docs]
def score_contributions(
self,
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).
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::
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.
A row with missing cells is scored over the cells it has, as in
:meth:`diagnose`: :math:`p_b^\top p_b` is then taken over the row's
cells in block :math:`b`, and :math:`p_\mathrm{super}^\top
p_\mathrm{super}` over the blocks it has a score in, so its
contributions still add up to its super score (#1144). They are NaN at
the missing cells, which add nothing to the sum.
See :meth:`MBPLS.score_contributions` for the parameter and return
descriptions; the API is identical.
"""
check_is_fitted(self, "block_loadings_")
deflated, sample_index = self._deflated_blocks(X, component)
a = int(component) - 1
p_s = self.super_loadings_.values[:, a]
sqrt_kb = self._block_scale()
# Each score is the fit's masked regression (_project_component), so each term divides by what
# it divides by: p_b'p_b over the row's cells in the block, and p_s'p_s over the blocks the row
# has a score in. On a complete row both are the loadings' unit length (#1144).
has_score = np.column_stack([_rows_with_data(deflated[name]) for name in self.block_names_])
super_divisor = row_regression_divisors(
np.where(
has_score,
0.0, # pragma: no mutate - row_regression_divisors reads only which cells are NaN
np.nan,
),
p_s,
)
raw: dict[str, np.ndarray] = {}
for b_idx, name in enumerate(self.block_names_):
p_b = self.block_loadings_[name].values[:, a]
divisor = row_regression_divisors(deflated[name], p_b) * super_divisor
weight = p_b / sqrt_kb[name] * p_s[b_idx]
# A column the fit left out is read as missing, NaN, but contributes 0, as in MBPLS, PCA and PLS (#1114)
raw[name] = np.where(self._left_out(name), 0.0, deflated[name] * weight / divisor[:, np.newaxis])
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()
}
[docs]
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, contributions in per_block.items():
# A missing cell adds nothing to its row's super score, so it adds nothing to the mean (#1144),
# where skipping its NaN averaged the column over the other rows alone.
frame = contributions.fillna(0.0)
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
[docs]
def super_score_plot(self, pc_horiz: int = 1, pc_vert: int = 2) -> go.Figure:
"""Scatter plot of MBPCA super-scores for two components."""
check_is_fitted(self, "super_scores_")
a_max = self.n_components_
if not (1 <= pc_horiz <= a_max and 1 <= pc_vert <= a_max):
raise ValueError(f"pc_horiz and pc_vert must be in 1..{a_max}.")
x = self.super_scores_[pc_horiz].to_numpy()
y = self.super_scores_[pc_vert].to_numpy()
labels = [str(i) for i in self.super_scores_.index]
fig = go.Figure(
data=[
go.Scatter(
x=x,
y=y,
mode="markers+text",
text=labels,
textposition="top center",
name="Super-scores",
)
]
)
fig.update_layout(
xaxis_title=f"t_super[{pc_horiz}]",
yaxis_title=f"t_super[{pc_vert}]",
title=f"MBPCA super-score plot: PC{pc_horiz} vs PC{pc_vert}",
)
return fig
[docs]
def super_loadings_bar_plot(self, component: int = 1) -> go.Figure:
"""Bar plot of MBPCA super-loadings for a single component."""
check_is_fitted(self, "super_loadings_")
a_max = self.n_components_
if not (1 <= component <= a_max):
raise ValueError(f"component must be in 1..{a_max}.")
loadings = self.super_loadings_[component]
fig = go.Figure(data=[go.Bar(x=list(loadings.index), y=loadings.to_numpy(), name=f"p_super[{component}]")])
fig.update_layout(
xaxis_title="Block",
yaxis_title=f"p_super[{component}]",
title=f"MBPCA super-loadings, component {component}",
)
return fig
[docs]
def display_results(self, show_cumulative: bool = True) -> str:
"""Format a short text summary of per-block R²X, iterations and timing."""
check_is_fitted(self, "super_scores_")
rows: list[str] = [f"MBPCA model: {self.n_components_} component(s), {len(self.block_names_)} X-block(s)"]
src = self.r2_x_per_block_cumulative_ if show_cumulative else self.r2_x_per_block_per_component_
header = ["PC", *(f"R²X[{name}]" for name in self.block_names_)]
body = [
[str(a + 1), *(f"{src.loc[name].iloc[a]:.4f}" for name in self.block_names_)]
for a in range(self.n_components_)
]
rows.extend(text_table(header, body))
rows.append("")
rows.append(f" Iterations per PC: {self.fitting_info_['iterations'].tolist()}")
rows.append(f" Time per PC (ms): {[round(float(t * 1000), 1) for t in self.fitting_info_['timing']]}")
return "\n".join(rows)
