# (c) Kevin Dunn, 2010-2026. MIT License. Based on own private work over the years.
"""T-PLS (PLS for T-shaped data structures) estimator and its data container (ENG-01).
Holds :class:`DataFrameDict`, the container for the partitionable (Z, F) and
static (Y) data blocks, and :class:`TPLS`, the sklearn-compatible T-PLS
regressor. ("TPLS" here means PLS for T-shaped data structures, not "Total PLS"
or "Three-way PLS".)
"""
from __future__ import annotations
import inspect
import logging
import numbers
import time
import typing
import warnings
from collections.abc import Callable, KeysView
from functools import partial
import numpy as np
import pandas as pd
from sklearn.base import BaseEstimator, RegressorMixin, _fit_context
from sklearn.metrics import (
mean_absolute_error,
mean_squared_error,
r2_score,
root_mean_squared_error,
)
from sklearn.utils import Bunch
from sklearn.utils.validation import check_array, check_is_fitted
from ._common import (
DiagnoseResult,
NotEnoughVarianceError,
SpecificationWarning,
_equal_values,
_model_method,
_nz,
_own_axes,
check_at_least_one,
check_variance_left,
check_y_varies,
column_sd,
components_to_fit,
constant_columns,
epsqrt,
largest_ssq_column,
refuse_empty_block,
resolve_loop_settings,
unobserved_columns,
user_stacklevel,
)
from ._diagnostics import _t2_score_sd
from ._limits import JudgedRows, spe_calculation
from ._limits import ellipse_coordinates as _ellipse_coordinates
from ._limits import hotellings_t2_limit as _hotellings_t2_limit
from ._nipals import internal_pls_nipals_fit_one_pc, nan_to_zeros, regress_a_space_on_b_row
from ._optimal_inversion import Probe, finite_or
from ._optimal_inversion import inversion_sensitivity as _inversion_sensitivity
from ._optimal_inversion import optimal_inversion as _optimal_inversion
from .plots import Plot
[docs]
class DataFrameDict(dict):
"""Container for the partitionable (Z, F) and static (Y) data blocks used by TPLS."""
[docs]
def __init__(self, datadict: dict[str, dict[str, pd.DataFrame]]):
"""
Initialize a DataFrameDict to handle partitionable and static dataframes.
datadict: Dictionary with 3 keys, one for each block: Z, F and Y.
Each block is itself a dictionary of dataframes: dict[str, dict[str, pd.DataFrame]]
"""
self.partitionable_blocks: list[str] = ["Z", "F", "Y"]
self.datadict: dict[str, dict[str, pd.DataFrame]] = {}
for block in self.partitionable_blocks:
self.datadict[block] = datadict.get(block, {})
first_group = next(iter(self.datadict["F"].keys()))
self.n_samples = self.datadict["F"][first_group].shape[0]
self.shape = (self.n_samples, len(self.datadict))
# Some basic checks: each dataframe inside each block has the same number of rows
for block in self.partitionable_blocks:
for group, df in self.datadict[block].items():
if not isinstance(df, pd.DataFrame):
raise TypeError(
f"Expected a DataFrame for block {block}, group '{group}'; got {type(df).__name__}."
)
if df.shape[0] != self.n_samples:
raise ValueError(
f"DataFrames in block {block} must have the same number of rows ({self.n_samples}). "
f"Group {group} has {df.shape[0]} rows."
)
[docs]
def keys(self) -> KeysView[str]: # type: ignore[override] # reason: narrows dict.keys to str view
"""Return the keys of the DataFrameDict."""
return self.datadict.keys()
def __setitem__(self, key: str, value: pd.DataFrame | dict) -> None:
"""Set a DataFrame for a specific key in the DataFrameDict."""
if key not in self.partitionable_blocks:
raise KeyError(f"Key {key} is not a valid partitionable block. Valid keys are: {self.partitionable_blocks}")
if not isinstance(value, pd.DataFrame):
raise TypeError(f"Expected a DataFrame for key {key}, got {type(value)}.")
if value.shape[0] != self.n_samples:
raise ValueError(
f"DataFrames in block {key} must have the same number of rows ({self.n_samples}). "
f"Provided DataFrame has {value.shape[0]} rows."
)
# The setter intentionally stores a bare DataFrame here (not the nested
# dict[str, DataFrame] used at construction); preserve that behaviour.
self.datadict[key] = value # type: ignore[assignment] # reason: legacy bare-DataFrame storage
@typing.overload
def __getitem__(self, lookup: str) -> dict[str, pd.DataFrame]:
pass
@typing.overload
def __getitem__(self, lookup: int | list[int] | list[np.int64] | np.ndarray | tuple) -> DataFrameDict:
pass
def __getitem__(
self, lookup: int | list[int] | list[np.int64] | str | np.ndarray | tuple
) -> DataFrameDict | dict[str, pd.DataFrame]:
"""Return a new DataFrameDict with partitioned data."""
if isinstance(lookup, str):
return self.datadict[lookup] # returns the `dict[str, pd.DataFrame]` version of the function
datadict: dict[str, dict[str, pd.DataFrame]] = {}
for block in self.partitionable_blocks:
datadict[block] = {}
for group, df in self.datadict[block].items():
match lookup:
case int() | np.integer():
datadict[block][group] = df.iloc[[lookup]]
case list():
datadict[block][group] = df.iloc[lookup] # type: ignore[index] # reason: list of np.int64
case np.ndarray():
datadict[block][group] = df.iloc[lookup.tolist()]
case tuple():
if lookup[1] == Ellipsis:
datadict[block][group] = df.iloc[[int(item) for item in lookup[0]]]
else:
raise TypeError(f"Invalid tuple structure for lookup: {lookup}")
case _:
raise TypeError(
f"Lookup must be an int, list of ints, or a string. Got {lookup}; {type(lookup)}"
)
return DataFrameDict(datadict)
def __len__(self):
"""Return the number of samples in the DataFrameDict."""
return self.n_samples
@staticmethod
def _value_equal(a: object, b: object) -> bool:
"""Compare two stored block values, which may be DataFrames or nested dicts of them."""
if isinstance(a, pd.DataFrame) and isinstance(b, pd.DataFrame):
return a.equals(b)
if isinstance(a, dict) and isinstance(b, dict):
return a.keys() == b.keys() and all(DataFrameDict._value_equal(a[k], b[k]) for k in a)
# A block is a dict of DataFrames, or a bare DataFrame set through __setitem__: anything else
# here is one of each, which are never equal.
return False
def __eq__(self, other: object) -> bool:
"""Value-based equality over the held data.
Two :class:`DataFrameDict` instances are equal when they share the same
block structure and every contained :class:`pandas.DataFrame` is equal
(via :meth:`pandas.DataFrame.equals`). The inherited ``dict.__eq__``
would instead compare the (always-empty) ``dict`` base - all data lives
in ``self.datadict`` - and so report every instance as equal regardless
of its contents (CodeQL ``py/missing-equals``).
"""
if self is other:
return True
if not isinstance(other, DataFrameDict):
return NotImplemented
return (
self.partitionable_blocks == other.partitionable_blocks
and self.datadict.keys() == other.datadict.keys()
and all(self._value_equal(self.datadict[b], other.datadict[b]) for b in self.datadict)
)
def __ne__(self, other: object) -> bool:
"""Negation of :meth:`__eq__`.
Defined explicitly because the C-level ``dict.__ne__`` would otherwise
bypass the Python ``__eq__`` above and compare the empty ``dict`` bases.
"""
result = self.__eq__(other)
return result if result is NotImplemented else not result
# All data lives in ``self.datadict`` and instances are mutable, so a
# ``DataFrameDict`` is deliberately unhashable (matching ``dict``); make the
# intent explicit now that ``__eq__`` is defined.
__hash__ = None # type: ignore[assignment] # reason: intentionally unhashable, mirrors dict
def __repr__(self):
"""Return a string representation of the DataFrameDict."""
groups_in_block_f = list(self.datadict["F"].keys())
groups_in_block_z = list(self.datadict["Z"].keys())
groups_in_block_y = list(self.datadict["Y"].keys())
output = f"DataFrameDict with {len(self)} samples and {len(self.datadict)} blocks: {list(self.datadict.keys())}"
output += f"\n F groups: {groups_in_block_f}"
output += f"\n Z groups: {groups_in_block_z}"
output += f"\n Y groups: {groups_in_block_y}"
return output
logger = logging.getLogger(__name__)
[docs]
class TPLS(RegressorMixin, BaseEstimator):
"""
TPLS algorithm for T-shaped data structures (we also include standard pre-processing of the data inside this class).
Source: Garcia-Munoz, https://doi.org/10.1016/j.chemolab.2014.02.006, Chem.Intell.Lab.Sys. v133, p 49 to 62, 2014.
We change the notation from the original paper to avoid confusion with a generic "X" matrix, and match symbols
that are more natural for our use.
Notation mapping (paper → this code):
- X^T → D: ``d_matrix`` (external), ``d_mats`` (internal) - Database of properties
- X → D^T: transposed D (not used directly)
- R → F: ``f_mats`` - Formula matrices
- Z → Z: ``z_mats`` - Process conditions
- Y → Y: ``y_mats`` - Quality indicators
Notes on the blocks:
1. Matrices in F, Z and Y must all have the same number of rows.
2. Columns in F must be the same as the **rows** in D. They are paired by name, in any order; a
material that no blend uses has a column of zeros in F.
3. Conditions in Z may be missing (turning it into an L-shaped data structure).
4. A blend with nothing observed in F, Z or Y, or a material with nothing observed in D, is fitted
with a zero score or loading in that block and no weight in that block's regressions, as PCA
and PLS fit an all-missing row.
Parameters
----------
n_components : int
A parameter used to specify the number of components.
d_matrix : dict[str, pd.DataFrame]
A dictionary containing the properties of each group of materials.
Keys are group names; values are DataFrames with properties as columns and materials as rows.
This "D" matrix is provided once at construction and reused for fitting, prediction and
cross-validation.
max_iter : int, optional
The maximum number of iterations for the TPLS algorithm. Default is 500. A component
that has not converged by then is the state of its last pass, and the fit warns.
skip_f_matrix_preprocessing : bool, optional
If True, the F (formula) matrices are used as-is, skipping the internal
centering and scaling of the F block. Default is False. The caller then owns
F's preprocessing: a column with no variance, which scaling would leave out with
a warning, is used as given, without one. A missing cell is still left out, and
a column with no observed value is still named.
tol : float, optional
Relative convergence tolerance for the per-component super-score loop: the
component is converged once the norm of the difference between two successive
super-score vectors, relative to the norm of the previous one, falls below
this. Default is ``epsqrt`` (about 1.5e-8), which is the value that was
hard-coded before this became a parameter, so the default fit is unchanged.
It must lie between machine precision and 10, as on every iterative estimator.
It is listed last so that no existing positional argument changed position.
Notes
-----
The input ``X`` passed to :meth:`fit` and :meth:`diagnose` is a dictionary with 3 keys:
- ``F``: Formula matrices (rows = blends, columns = materials).
``F = {"Group A": df_formulas_a, "Group B": df_formulas_b, ...}``
- ``Z``: Process conditions - one row per blend, one column per condition, in one or more
groups. Each Z group has its own weights and block score, as a block of a multi-block PLS.
- ``Y``: Product quality indicators - one row per blend, one column per indicator.
The ``D`` matrix (database of material properties) is supplied once at
construction via the ``d_matrix`` argument; it is not part of ``X``.
Attributes
----------
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.
max_iter_ : int
The ``max_iter`` the model was fitted with, which :meth:`display_results` reports.
skip_f_matrix_preprocessing_ : bool
Whether the fit used the F block as given; :meth:`diagnose` treats new F blocks the same way.
n_samples : int
Number of samples (rows) in the training data.
n_substances : int
Number of materials (columns) in the F matrix, summed across groups.
n_conditions : int
Number of process-condition columns summed across Z blocks.
n_outputs : int
Number of quality-indicator columns summed across Y blocks.
is_fitted_ : bool
Set to True once ``fit()`` completes.
fitting_statistics : dict
Per-component ``iterations``, ``convergance_tolerance``, ``milliseconds`` and
``converged`` lists. ``converged`` holds whether the component met ``tol`` within
``max_iter`` passes, as ``fitting_info_["converged"]`` does on MBPLS and MBPCA.
preproc_ : dict
Nested per-block per-group preprocessors, indexed as
``preproc_[block][group]`` where ``block`` is ``"D"``, ``"F"``,
``"Y"`` or ``"Z"``. With ``skip_f_matrix_preprocessing`` the ``"F"``
entries are learned but not applied.
sums_of_squares_ : list[dict]
Per-component per-block sums of squares.
r2_frac : list[dict]
Per-component per-block cumulative R^2 fractions.
feature_importance : dict
Per-block per-group variable-importance (populated for ``"D"``
and ``"F"``).
d_mats, f_mats, z_mats, y_mats : dict[str, np.ndarray]
Deflated block matrices (per group for D/F, per Z-block / Y-block
for Z/Y).
not_na_d, not_na_f, not_na_z, not_na_y : dict[str, np.ndarray]
Boolean observed-value masks matching the shapes above. A column excluded for having
no variance is marked absent in every row.
observation_names : pd.Index
Row index shared across the F/Z/Y blocks.
property_names, material_names : dict[str, list[str]]
Column and row names of each D group.
condition_names, quality_names : dict[str, list[str]]
Column names of each Z-block and Y-block respectively.
t_scores_super : pd.DataFrame, shape (n_samples, n_components)
Super-block scores ``T``.
r_loadings_f : dict[str, pd.DataFrame]
F-block weights per group.
w_loadings_z : dict[str, pd.DataFrame]
Z-block weights per Z-block.
w_loadings_super : pd.DataFrame
Super-block weights, rows ``["Z", "F"]`` (or just ``["F"]`` when
there are no process conditions). Each Z group is a block of its own
in the super-level, so several Z groups give one row each, labelled
``Z["<group>"]``, before the ``"F"`` row.
s_loadings_d, v_loadings_d : dict[str, pd.DataFrame]
D-block scores and loadings per group.
p_loadings_f : dict[str, pd.DataFrame]
F-block loadings per group (used for deflation).
p_loadings_z : dict[str, pd.DataFrame]
Z-block loadings per Z-block.
q_loadings_y : dict[str, pd.DataFrame]
Y-block loadings per Y-block.
hat_ : dict[str, pd.DataFrame]
Per-Y-block predictions on the preprocessed (centred / scaled) scale.
hat : dict[str, pd.DataFrame]
Per-Y-block in-sample predictions on the *original* scale.
spe : dict[str, dict[str, pd.DataFrame or pd.Series]]
Nested ``spe[block][group]`` squared prediction error tables.
spe_limit : dict[str, dict[str, Callable]]
Nested ``spe_limit[block][group]`` callables. Each is a
:func:`functools.partial` over :func:`spe_calculation` bound to
the block's SPE array; call it with a confidence level to obtain
the SPE limit. This is a deliberate divergence from the flat
``spe_limit`` method on PCA / PLS.
hotellings_t2 : pd.DataFrame, shape (n_samples, n_components)
Cumulative Hotelling's T^2 per super-component: column ``a`` holds
``sum_{j<=a} (t_j / s_j)^2``, the same form as PCA / PLS and as
:meth:`diagnose` (#502).
scaling_factor_for_scores : pd.Series
Standard deviation of each super-score column, computed with the
unbiased (N - 1) divisor; used by the ellipse / T^2 helpers.
.. note::
TPLS deliberately does **not** follow the sklearn trailing-
underscore convention for its fitted attributes. Names such as
``t_scores_super``, ``spe``, ``hotellings_t2``, and the
``*_loadings_*`` family are written without the trailing ``_``
to keep the chemometrics symbol names readable. Attributes
that are set in ``fit`` for the class's own
bookkeeping (``is_fitted_``, ``preproc_``, ``required_blocks_``)
do carry the underscore.
Example
-------
>>> import numpy as np
>>> import pandas as pd
>>> rng = np.random.default_rng()
>>>
>>> n_props_a, n_props_b = 6, 4 # Two groups of properties: A and B.
>>> n_materials_a, n_materials_b = 12, 8 # Number of materials in each group.
>>> n_formulas = 40 # Number of formulas in matrix F.
>>> n_outputs = 3
>>> n_conditions = 2
>>>
>>> properties = {
>>> "Group A": pd.DataFrame(rng.standard_normal((n_materials_a, n_props_a))),
>>> "Group B": pd.DataFrame(rng.standard_normal((n_materials_b, n_props_b))),
>>> }
>>> formulas = {
>>> "Group A": pd.DataFrame(rng.standard_normal((n_formulas, n_materials_a))),
>>> "Group B": pd.DataFrame(rng.standard_normal((n_formulas, n_materials_b))),
>>> }
>>> process_conditions = {"Conditions": pd.DataFrame(rng.standard_normal((n_formulas, n_conditions)))}
>>> quality_indicators = {"Quality": pd.DataFrame(rng.standard_normal((n_formulas, n_outputs)))}
>>> all_data = {"Z": process_conditions, "F": formulas, "Y": quality_indicators}
>>> estimator = TPLS(n_components=4, d_matrix=properties)
>>> estimator.fit(DataFrameDict(all_data))
"""
# This is a dictionary allowing to define the type of parameters.
# It used to validate parameter within the `_fit_context` decorator.
_parameter_constraints: typing.ClassVar = {
"n_components": [numbers.Integral],
"max_iter": [numbers.Integral],
"tol": [numbers.Real],
"d_matrix": [dict, None],
}
def __init__(
self,
n_components: int,
d_matrix: dict,
max_iter: int = 500,
skip_f_matrix_preprocessing: bool = False,
# Appended rather than slotted in beside ``max_iter`` so that no existing
# positional argument changes position (#588).
tol: float = epsqrt,
):
super().__init__()
# Stored unchecked, as sklearn requires: fit checks them, however they were set (#1071).
self.n_components = n_components
self.d_matrix = d_matrix # This is required input dict containing the properties for each group.
self.max_iter = max_iter
self.tol = tol
self.skip_f_matrix_preprocessing = skip_f_matrix_preprocessing
def __sklearn_is_fitted__(self) -> bool:
"""Whether ``fit`` has completed, for :func:`~sklearn.utils.validation.check_is_fitted`.
A fit that fails part way has set some attributes that end in an underscore, so the
check scikit-learn makes by default, that one exists, would pass on that model.
"""
return getattr(self, "is_fitted_", False)
@property
def plot(self) -> Plot:
"""Score and loading plots of the fitted model: ``model.plot.scores()`` and ``model.plot.loadings()``."""
return Plot(self)
# ENG-05: convenience methods forwarding to the standalone functions. These
# used to be ``functools.partial`` instances bound in ``fit``; defining them
# as real methods keeps ``help`` / ``inspect.signature`` accurate and the
# fitted model picklable.
# Inversion as an optimisation: which lots, at what fractions, and which conditions, under rules.
optimal_inversion = _model_method(_optimal_inversion)
inversion_sensitivity = _model_method(_inversion_sensitivity)
def _inversion_columns(self) -> list[tuple[str, str, list[str]]]:
"""Each input group, in the order optimal inversion reads them: (block, group, its column names)."""
columns = [("F", g, list(self.material_names[g])) for g in self.material_names]
return columns + [("Z", g, list(self.condition_names[g])) for g in self.condition_names]
def _inversion_probe(self, conf_level: float) -> Probe:
"""
Describe this model for optimal inversion: the inputs are each group's lot fractions, then the process
conditions; each F and Z group is a residual block; rows are read through diagnose.
A group the model uses no column of, as a group with no variation, which has zero weight (#1144), is an
input but not a residual block: its SPE is 0 whatever the inputs, so it is no rule, and its limit of 0
would divide the search by 0.
"""
columns = self._inversion_columns()
present = {"F": self.not_na_f, "Z": self.not_na_z}
residual_blocks = [(block, group) for block, group, _ in columns if present[block][group].any()]
names = tuple(name for _, _, group in columns for name in group)
if len(set(names)) != len(names):
# shortcut: inputs are named by their column name alone; qualify them by group if two groups share a name.
msg = "Two TPLS groups share a lot or condition name; optimal inversion needs every input named once."
raise ValueError(msg)
edges = np.cumsum([0] + [len(group) for _, _, group in columns])
# The pragmas below mark the changes no result can show: the rows' names, which diagnose's arrays are read
# without, the zip's strict, as there is one edge more than there are groups, and the float casts of
# frames that hold float64 already.
def diagnose(rows: np.ndarray) -> tuple[np.ndarray, np.ndarray, list[np.ndarray], np.ndarray, np.ndarray]:
index = [f"p{i}" for i in range(len(rows))] # pragma: no mutate - the arrays are read by position
blocks: dict[str, dict[str, pd.DataFrame]] = {}
for (block, group, cols), lo, hi in zip(
columns,
edges[:-1],
edges[1:],
strict=True, # pragma: no mutate - one edge more than there are groups
):
blocks.setdefault(block, {})[group] = pd.DataFrame(
rows[:, lo:hi],
index=index, # pragma: no mutate - the arrays are read by position
columns=cols,
)
out = self.diagnose(DataFrameDict(blocks))
y = np.hstack(
[out.hat[g].to_numpy(dtype=float) for g in self.quality_names] # pragma: no mutate - float64 already
)
residual = [
out.residuals[b][g].to_numpy(dtype=float) # pragma: no mutate - float64 already
for b, g in residual_blocks
]
spe = np.column_stack(
[
out.spe[b][g]
.iloc[:, -1]
.to_numpy(
dtype=float # pragma: no mutate - float64 already
)
for b, g in residual_blocks
]
)
t2 = out.cumulative_hotellings_t2.iloc[:, -1].to_numpy(
dtype=float # pragma: no mutate - float64 already
)
return (
out.t_scores_super.to_numpy(dtype=float), # pragma: no mutate - float64 already
y,
residual,
t2,
spe,
)
def preproc(key: str, block: str, group: str, cols: list[str]) -> np.ndarray:
return (
self.preproc_[block][group][key]
.reindex(cols)
.to_numpy(
dtype=float # pragma: no mutate - float64 already
)
)
return Probe(
names=names, y_names=tuple(n for g in self.quality_names for n in self.quality_names[g]),
groups={g: tuple(self.material_names[g]) for g in self.material_names},
block_names=tuple(f"{b}/{g}" for b, g in residual_blocks),
block_limits=tuple(float(self.spe_limit[b][g](conf_level)) for b, g in residual_blocks),
score_sd=_t2_score_sd(self.scaling_factor_for_scores.to_numpy()),
t2_limit=float(self.hotellings_t2_limit(conf_level, judged_rows="new")), # an input is a new row
y_scale=np.concatenate([finite_or(self.preproc_["Y"][g]["scale"], 1.0) for g in self.quality_names]),
centre=np.concatenate([np.nan_to_num(preproc("center", b, g, c)) for b, g, c in columns]),
spread=np.concatenate([finite_or(preproc("scale", b, g, c), 1.0) for b, g, c in columns]),
diagnose=diagnose,
) # fmt: skip
[docs]
def hotellings_t2_limit(self, conf_level: float = 0.95, *, judged_rows: JudgedRows | None = None) -> float:
"""Hotelling's T2 limit at the given confidence level (see :func:`hotellings_t2_limit`).
``judged_rows`` is ``"fitted"`` to judge the rows the model was fitted to (its
``hotellings_t2``), or ``"new"`` to judge rows it was not fitted to. ``None`` gives the
``"new"`` limit until 3.0 and the ``"fitted"`` one from 3.0, with a :class:`FutureWarning`.
"""
return _hotellings_t2_limit(
conf_level=conf_level,
n_components=self.n_components_,
n_rows=self.hotellings_t2.shape[0],
judged_rows=judged_rows,
)
[docs]
def ellipse_coordinates(
self,
score_horiz: int,
score_vert: int,
conf_level: float = 0.95,
n_points: int = 100,
*,
judged_rows: JudgedRows | None = None,
) -> tuple[np.ndarray, np.ndarray]:
"""Coordinates of the T2 confidence ellipse (see :func:`ellipse_coordinates`)."""
return _ellipse_coordinates(
score_horiz=score_horiz,
score_vert=score_vert,
conf_level=conf_level,
n_points=n_points,
n_components=self.n_components_,
scaling_factor_for_scores=self.scaling_factor_for_scores,
n_rows=self.t_scores_super.shape[0],
judged_rows=judged_rows,
)
[docs]
@_fit_context(prefer_skip_nested_validation=True)
def fit(self, X: DataFrameDict, y: object = None) -> TPLS: # noqa: PLR0915
"""Fit the preprocessing parameters and also the latent variable model from the training data.
Parameters
----------
X : DataFrameDict
The training input samples. A plain ``dict`` raises :class:`TypeError`. Keys that
must be present: ``"F"``, ``"Z"``, and ``"Y"``. See the class docstring for the
shape of each block.
y : object, optional
Must be ``None``. The signature exists only for sklearn API compatibility;
a T-shaped model takes its response from ``X["Y"]``. Anything else raises,
rather than being silently ignored as it was before #565.
Returns
-------
self : object
Returns self.
Raises
------
ValueError
If ``y`` is not ``None``, or if ``n_components``, ``max_iter`` or ``tol`` is out of range, or
if a material of a group of F is not a row of that group of ``d_matrix``, or the reverse.
TypeError
If ``X`` is not a :class:`DataFrameDict`, or ``d_matrix`` is not a dict of DataFrames.
NotEnoughVarianceError
If the first Y block has no variance, or no group has a column of F and a column of D that
vary, as MBPLS raises for an X with no variance; if the first Y block has no variance left
for a component, as when the components before it fit Y exactly, in the words of PLS and
MBPLS; or if a component's scores are not finite because a block has no variation for the
iteration to follow. A group whose F or D has no
variation, beside one whose F and D vary, is fitted with zero weight (#1144): its weights and
loadings are 0, and a warning names it.
Warns
-----
SpecificationWarning
If a component did not converge within ``max_iter`` passes (#1163).
"""
_reject_separate_y(y, "TPLS.fit")
if not isinstance(X, DataFrameDict):
raise TypeError(f"X must be a DataFrameDict; got {type(X).__name__}.")
check_at_least_one("n_components", self.n_components)
# No missing_data_settings, so the resolver's message naming the estimator is never reached.
settings = resolve_loop_settings(tol=self.tol, max_iter=self.max_iter, missing_data_settings=None)
self._input_data_checks(X)
# F in D's order from here on, so the centre and scale learned for F are in the order diagnose reads it.
X = DataFrameDict(
{**X.datadict, "F": {k: _f_in_the_order_of_d(f, self.d_matrix[k], k) for k, f in X["F"].items()}}
)
# A fit that fails part way leaves the model unfitted, not the previous fit with some parts replaced.
self.is_fitted_ = False
self.required_blocks_ = {"D", "F", "Y", "Z"} # "Z" block is optional; an empty one is added if not provided
group_keys = list(self.d_matrix)
# What the fit used, which everything after it reads, so that set_params without a refit changes nothing.
self.max_iter_ = settings["md_max_iter"]
self.skip_f_matrix_preprocessing_ = self.skip_f_matrix_preprocessing
# Storage for pre-processing and the raw matrices
self.fitting_statistics: dict[str, list] = {
"iterations": [],
"convergance_tolerance": [],
"milliseconds": [],
"converged": [],
}
# Per block / per group preprocessing parameters. "center" and "scale" are pd.Series (one entry
# per column); the D-block additionally stores a scalar "block" factor (float). The heterogeneous
# value types make `Any` the honest inner annotation here. See `_preprocess`.
self.preproc_: dict[str, dict[str, dict[str, typing.Any]]] = {key: {} for key in self.required_blocks_}
self.feature_importance: dict[str, dict[str, pd.Series]] = {key: {} for key in self.required_blocks_}
# Float copies: int64 columns are accepted, and the blocks are deflated in place, which
# cannot write a float result back into an integer array.
self.d_mats: dict[str, np.ndarray] = {
key: self.d_matrix[key].to_numpy(dtype=float, copy=True) for key in group_keys
}
self.f_mats: dict[str, np.ndarray] = {key: X["F"][key].to_numpy(dtype=float, copy=True) for key in group_keys}
self.z_mats: dict[str, np.ndarray] = {key: X["Z"][key].to_numpy(dtype=float, copy=True) for key in X["Z"]}
self.y_mats: dict[str, np.ndarray] = {key: X["Y"][key].to_numpy(dtype=float, copy=True) for key in X["Y"]}
# Empty model coefficients
self.n_substances = sum(self.f_mats[key].shape[1] for key in group_keys)
self.n_conditions = sum(self.z_mats[key].shape[1] for key in self.z_mats)
self.n_outputs = sum(self.y_mats[key].shape[1] for key in self.y_mats)
self.n_samples = self.f_mats[group_keys[0]].shape[0]
# Learn the centering and scaling parameters
for key in X["Y"]:
self.preproc_["Y"][key] = {}
self.preproc_["Y"][key]["center"], self.preproc_["Y"][key]["scale"] = (
self._learn_center_and_scaling_parameters(X["Y"][key], f'Y["{key}"]')
)
for key in X["Z"]:
self.preproc_["Z"][key] = {}
self.preproc_["Z"][key]["center"], self.preproc_["Z"][key]["scale"] = (
self._learn_center_and_scaling_parameters(X["Z"][key], f'Z["{key}"]')
)
for key, df_d in self.d_matrix.items():
self.preproc_["D"][key] = {}
self.preproc_["D"][key]["center"], self.preproc_["D"][key]["scale"] = (
self._learn_center_and_scaling_parameters(df_d, f'D["{key}"]')
)
# Block-scale each D-block by 1 / sqrt(P_i * M_i), where P_i = number of lots (rows) and
# M_i = number of properties (columns). After column-wise auto-scaling this makes
# trace(X_i^T X_i) ~= 1 for every block, removing bias toward blocks that simply have more lots
# or more properties. See Garcia-Munoz (2014), section 2.1 (https://doi.org/10.1016/j.chemolab.2014.02.006).
n_lots, n_properties = df_d.shape
self.preproc_["D"][key]["block"] = float(np.sqrt(n_lots * n_properties))
#
# Also do the same for the formula matrix. F used as given is not scaled, so its centre and
# scale are learned only for optimal inversion, which starts and steps from them.
self.preproc_["F"][key] = {}
self.preproc_["F"][key]["center"], self.preproc_["F"][key]["scale"] = (
self._learn_center_and_scaling_parameters(X["F"][key], f'F["{key}"]')
)
# Then implement the preprocessing on the raw data
self._preprocess_data()
# Sum of square values for each column in each block (dicts) per component (elements in the list)
# The first entry is after centering and scale (baseline variance) but before fitting any components.
# The second entry is after fitting one component, and so on.
# You can sum the sums-of-squares values for all columns to get the total variance for each block.
self.sums_of_squares_: list[dict[str, dict[str, np.ndarray]]] = [
{
"D": {key: np.nansum(self.d_mats[key] ** 2, axis=0) for key in group_keys},
"F": {key: np.nansum(self.f_mats[key] ** 2, axis=0) for key in group_keys},
"Z": {key: np.nansum(self.z_mats[key] ** 2, axis=0) for key in X["Z"]},
"Y": {key: np.nansum(self.y_mats[key] ** 2, axis=0) for key in X["Y"]},
}
]
# These are *fractional* R2 values, i.e. always less than or equal to 1.0.
# As a list: entry 0 is zeros; entry 1 is after fitting the first component, and so on.
# The keys are the blocks, and the values are dictionaries with group keys as keys.
# The values are the R2 values for each column in the block.
self.r2_frac: list[dict[str, dict[str, np.ndarray]]] = [
{
"D": {key: np.zeros(self.d_mats[key].shape[1]) for key in group_keys},
"F": {key: np.zeros(self.f_mats[key].shape[1]) for key in group_keys},
"Z": {key: np.zeros(self.z_mats[key].shape[1]) for key in self.z_mats},
"Y": {key: np.zeros(self.y_mats[key].shape[1]) for key in self.y_mats},
}
]
# Create the missing value maps, except we store the opposite, i.e., not missing, since these are more useful.
# We refer to these as `pmaps` in the code (present maps, as opposed to `mmap` or missing maps). They are
# taken from the preprocessed blocks, not the raw ones, so a column excluded for having no variance (NaN
# after scaling) is absent too: it is zeroed below, and counted as present it shrank every score (#513).
self.not_na_d = {key: ~np.isnan(self.d_mats[key]) for key in group_keys}
self.not_na_f = {key: ~np.isnan(self.f_mats[key]) for key in group_keys}
self.not_na_z = {key: ~np.isnan(self.z_mats[key]) for key in X["Z"]}
self.not_na_y = {key: ~np.isnan(self.y_mats[key]) for key in X["Y"]}
# Then set missing data values to zeros (not because we are ignoring the values), but because we will use
# the missing value maps to identify where the missing values are and therefore ignore them. But set to zero,
# so these values have no influence on the calculations.
self.d_mats = {key: nan_to_zeros(self.d_mats[key]) for key in group_keys}
self.f_mats = {key: nan_to_zeros(self.f_mats[key]) for key in group_keys}
self.z_mats = {key: nan_to_zeros(self.z_mats[key]) for key in X["Z"]}
self.y_mats = {key: nan_to_zeros(self.y_mats[key]) for key in X["Y"]}
self._refuse_blocks_with_no_variance()
self._name_blocks_with_no_variance(X)
# Storage for the model objects. Make a copy only of the Numpy values to use in the Estimator.
# A copy, so that a name the caller gives this axis later does not rename the model's frames
self.observation_names = X["F"][group_keys[0]].index.copy()
self.property_names = {key: self.d_matrix[key].columns.to_list() for key in group_keys}
self.material_names = {key: self.d_matrix[key].index.to_list() for key in group_keys}
self.condition_names = {key: X["Z"][key].columns.to_list() for key in X["Z"]}
self.quality_names = {key: X["Y"][key].columns.to_list() for key in X["Y"]}
# Model parameters. Naming convention: x_i_j
# x = block letter (P, W, R, T, etc)
# i = block type: `scores` [for the observations (rows)] or `loadings` [for the variables (columns)]
# j = block name [z, f, d, y, super]
# ----------------
self.t_scores_super: pd.DataFrame = pd.DataFrame(index=self.observation_names)
# The frames start from copies of the blocks' own axes, not the name lists, so each keeps its axis
# name (#1083), and renaming the caller's axis later does not rename the model's.
self.r_loadings_f: dict[str, pd.DataFrame] = {
key: pd.DataFrame(index=self.d_matrix[key].index.copy()) for key in group_keys
}
self.w_loadings_z: dict[str, pd.DataFrame] = {
key: pd.DataFrame(index=X["Z"][key].columns.copy()) for key in self.z_mats
}
self.w_loadings_super = pd.DataFrame(index=self._super_block_names())
# Capture the correlation of the properties in D; for the last component.
self.s_loadings_d: dict[str, pd.DataFrame] = {
key: pd.DataFrame(index=self.d_matrix[key].columns.copy()) for key in group_keys
}
# Captures the deflation of the properties in D; for the last component.
self.v_loadings_d: dict[str, pd.DataFrame] = {
key: pd.DataFrame(index=self.d_matrix[key].columns.copy()) for key in group_keys
}
self.p_loadings_f: dict[str, pd.DataFrame] = {
key: pd.DataFrame(index=self.d_matrix[key].index.copy()) for key in group_keys
}
self.p_loadings_z: dict[str, pd.DataFrame] = {
key: pd.DataFrame(index=X["Z"][key].columns.copy()) for key in self.z_mats
}
self.q_loadings_y: dict[str, pd.DataFrame] = {
key: pd.DataFrame(index=X["Y"][key].columns.copy()) for key in self.y_mats
}
# Model performance
# -----------------
# 1. Prediction matrices (hat matrices for Y-space) in pre-processed space
self.hat_: dict[str, pd.DataFrame] = {
key: pd.DataFrame(index=self.observation_names, columns=X["Y"][key].columns.copy(), dtype=float).fillna(0)
for key in self.y_mats
}
# 2. Squared prediction error (SPE) for each observation, per component, per block
# Most blocks store a (N x 1) DataFrame of SPE values; the D-space stores a Series.
self.spe: dict[str, dict[str, pd.DataFrame | pd.Series]] = {key: {} for key in self.required_blocks_}
self.spe_limit: dict[str, dict[str, Callable]] = {key: {} for key in self.required_blocks_}
# Every component's weights on the materials come through D, which runs out of directions after as
# many components as its rank: later ones were fitted to rounding noise, with scores near 1e15 (#1149).
self.n_components_ = components_to_fit(
self.n_components,
self.n_samples,
max(min(d.shape) for d in self.d_matrix.values()),
columns="the largest rank a group of D can have, the fewer of its materials and its properties",
)
self._fit_iterative_regressions()
self.is_fitted_ = True
return self
[docs]
def diagnose(self, X: DataFrameDict) -> Bunch: # noqa: C901
"""
Model inference on new data.
This will pre-process the new data and apply those subsequently to the latent variable model.
Example
-------
# Training phase:
estimator = TPLS(n_components=2).fit(training_data)
# Testing/inference phase. A plain dict raises TypeError, so wrap it.
# You need at least the F block for a new prediction; "Z" is optional.
new_data = DataFrameDict({"Z": ..., "F": ...})
predictions = estimator.diagnose(new_data)
Parameters
----------
X : DataFrameDict
The input samples.
Returns
-------
sklearn.utils.Bunch
A bunch with the following fields:
- ``hat`` (dict[str, DataFrame]) : predicted Y per Y-group, on
the original (un-scaled) scale, indexed by the observation
names of ``X``.
- ``t_scores_super`` (DataFrame, shape ``(n_new, n_components)``)
: super-scores for the new observations.
- ``spe`` (dict[str, dict[str, DataFrame]]) : per-component
squared prediction error for the ``"Z"`` and ``"F"`` blocks,
keyed first by block name and then by group; each inner
DataFrame has shape ``(n_new, n_components)``.
- ``cumulative_hotellings_t2`` (DataFrame, shape
``(n_new, n_components)``) : cumulative Hotelling's T² per
new observation after each component, as in ``hotellings_t2``,
the fitted attribute. Every ``diagnose`` of the family returns it
in this shape (#983).
- ``hotellings_t2`` : deprecated since 2.4.0, to be removed in 3.0;
the same DataFrame as ``cumulative_hotellings_t2``. Reading it warns.
- ``residuals`` (dict[str, dict[str, DataFrame]]) : what is left of
each ``"Z"`` and ``"F"`` group after all the components, in the
preprocessed (centred and scaled) units, keyed like ``spe``. The
last column of ``spe`` is the length of each row; a missing cell,
and a column excluded for having no variance, leave 0.
"""
check_is_fitted(self) # Check if fit had been called
if not isinstance(X, DataFrameDict):
raise TypeError(f"X must be a DataFrameDict; got {type(X).__name__}.")
# DataFrameDict holds every group of every block at one row count (its constructor
# and __setitem__ both enforce it), so only the column names need checking here.
for key, df_f in X["F"].items():
if set(df_f.columns) != set(self.material_names[key]):
raise ValueError(
f"Columns in block F, group [{key}] must match training data column names for each material."
)
for key, df_z in X["Z"].items():
if set(df_z.columns) != set(self.condition_names[key]):
raise ValueError(f"Column names in block Z, group [{key}] must match training data column names.")
# In the training column order, since the scores below read the values by position: F used as
# given was scored in the order it came. As floats, so a nullable column's NA reads as missing.
x_f = {key: frame[self.material_names[key]].to_numpy(dtype=float) for key, frame in X["F"].items()}
x_z = {key: frame[self.condition_names[key]].to_numpy(dtype=float) for key, frame in X["Z"].items()}
for key, values in x_f.items():
if not self.skip_f_matrix_preprocessing_:
x_f[key] = (values - self.preproc_["F"][key]["center"].to_numpy()) / self.preproc_["F"][key][
"scale"
].to_numpy()
for key, values in x_z.items():
x_z[key] = (values - self.preproc_["Z"][key]["center"].to_numpy()) / self.preproc_["Z"][key][
"scale"
].to_numpy()
# From the preprocessed blocks, as in fit(), so a column excluded for having no variance is absent.
# A missing cell is then 0, so it adds nothing to a score's numerator, as in fit().
not_na_f = {key: ~np.isnan(values) for key, values in x_f.items()}
not_na_z = {key: ~np.isnan(values) for key, values in x_z.items()}
x_f = {key: np.where(not_na_f[key], values, 0.0) for key, values in x_f.items()}
x_z = {key: np.where(not_na_z[key], values, 0.0) for key, values in x_z.items()}
names_observations = X["F"][next(iter(X["F"]))].index
n_rows, n_components = len(names_observations), self.n_components_
scores = np.empty((n_rows, n_components))
spe = {"Z": {key: np.empty((n_rows, n_components)) for key in x_z}, "F": {key: scores.copy() for key in x_f}}
weights_super = self.w_loadings_super.to_numpy()
for pc_a in range(n_components):
# The F score: each row of all the F groups at once regressed on r, over its observed cells. A
# row with nothing observed in F has no F score: NaN (0 / 0), as PCA and PLS transform an
# all-missing row.
r = {key: self.r_loadings_f[key].to_numpy()[:, pc_a] for key in x_f}
with np.errstate(invalid="ignore"):
score_f = sum(x_f[key] @ r[key] for key in x_f) / sum(not_na_f[key] @ r[key] ** 2 for key in x_f)
# One Z score per Z group likewise, from that group's own weights, as fit() scores each group on its
# own (#981), in the training order of the rows of w_loadings_super; none for a model without Z. A
# row with no value in a group scores 0 there, as in fit() (#513). Z is optional here: a group left
# out of X scores 0.
score_z = {key: np.zeros(n_rows) for key in self.w_loadings_z}
for key, values in x_z.items():
w = self.w_loadings_z[key].to_numpy()[:, pc_a]
rows = not_na_z[key].any(axis=1)
score_z[key][rows] = values[rows] @ w / (not_na_z[key][rows] @ w**2)
# The super score: the block scores, [Z groups..., F] or [F], weighted by the super weights.
scores[:, pc_a] = np.column_stack([*score_z.values(), score_f]) @ weights_super[:, pc_a]
# Deflate each block, to score the next component and to compute SPE. Only the observed cells,
# as in fit(): a missing cell has no residual (#1028).
for block, values, observed, loadings in (
("F", x_f, not_na_f, self.p_loadings_f),
("Z", x_z, not_na_z, self.p_loadings_z),
):
for key in values:
explained = np.outer(scores[:, pc_a], loadings[key].to_numpy()[:, pc_a])
values[key] = values[key] - explained * observed[key]
spe[block][key][:, pc_a] = np.sqrt(np.sum(values[key] ** 2, axis=1))
# Predictions are returned in un-scaled form, so they are in the same units as the training data.
components = list(range(1, n_components + 1))
hat = {
key: pd.DataFrame(
scores @ self.q_loadings_y[key].to_numpy().T * self.preproc_["Y"][key]["scale"].to_numpy()
+ self.preproc_["Y"][key]["center"].to_numpy(),
index=names_observations,
columns=self.q_loadings_y[key].index,
)
for key in self.y_mats
}
# The T2 values, with the same cumulative per-component form (and the same rank-deficiency
# guard) as ``fit`` stores in ``self.hotellings_t2``, so the two agree for the same rows (#502).
t2 = pd.DataFrame(self._cumulative_hotellings_t2(scores), index=names_observations, columns=components)
return DiagnoseResult(
f"multivariate.{type(self).__name__}",
hat=hat,
t_scores_super=pd.DataFrame(scores, index=names_observations, columns=components),
spe={
block: {
key: pd.DataFrame(values, index=X[block][key].index, columns=components)
for key, values in spe[block].items()
}
for block in ("Z", "F")
},
cumulative_hotellings_t2=t2,
hotellings_t2=t2,
# The blocks as deflated by the last component: the vectors whose lengths are SPE, for optimal inversion.
residuals={
"F": {
key: pd.DataFrame(values, index=X["F"][key].index, columns=self.p_loadings_f[key].index)
for key, values in x_f.items()
},
"Z": {
key: pd.DataFrame(values, index=X["Z"][key].index, columns=self.p_loadings_z[key].index)
for key, values in x_z.items()
},
},
)
[docs]
def display_results(self, show_cumulative_stats: bool = True) -> str:
"""Display the results of the model fitting."""
if not self.__sklearn_is_fitted__():
raise RuntimeError("The model is not fitted yet. Please call `fit` first.")
t2_limits = [self.hotellings_t2_limit(level, judged_rows="new") for level in (0.95, 0.99)]
output = f"Hotelling's T2 limit for new rows [95% limit]: {t2_limits[0]:.4g}\n"
output += f" [99% limit]: {t2_limits[1]:.4g}\n"
# output += f"SPE limits: {self.spe_limit['Y'](self.spe['Y'])}\n"
sep = "------ ---------- ---------- ---------- ---------- -------------\n"
output += sep
if show_cumulative_stats:
header = "LV # sum(R2: D) sum(R2: Z) sum(R2: F) sum(R2: Y)| ms [iter]"
else:
header = "LV # R2: D R2: Z R2: F R2: Y| ms [iter]"
output += header + "\n" + sep
# Each block's R2 per component: the mean, over its groups, of the mean R2 of a group's columns.
r2 = {
block: np.array([np.mean([np.nanmean(v) for v in entry[block].values()]) for entry in self.r2_frac[1:]])
for block in ["D", "F", "Y"] + (["Z"] if self.n_conditions > 0 else [])
}
if show_cumulative_stats:
r2 = {block: np.cumsum(values) for block, values in r2.items()}
for a, iterations in enumerate(self.fitting_statistics["iterations"]):
time_ms = self.fitting_statistics["milliseconds"][a]
r2_z = f"{r2['Z'][a] * 100:>10.1f}" if "Z" in r2 else f"{'-':>10}"
time_iter = f"{time_ms:>5.0f} [{iterations:>4d}]"
line = (
f"LV {a + 1:<2} {r2['D'][a] * 100:>10.1f} {r2_z} {r2['F'][a] * 100:>10.1f} {r2['Y'][a] * 100:>10.1f}"
)
line += f"|{time_iter:>13}"
if iterations >= self.max_iter_:
line += "** (max iter reached)"
output += line + "\n"
output += sep
ms_per_iter = round(
sum(self.fitting_statistics["milliseconds"]) / sum(self.fitting_statistics["iterations"]), 2
)
output += f"Timing: {ms_per_iter} ms/iter; {sum(self.fitting_statistics['iterations'])} iterations required\n"
output += f"Total time: {sum(self.fitting_statistics['milliseconds']) / 1000:.2f} seconds\n"
output += f"Average tolerance: {np.mean(self.fitting_statistics['convergance_tolerance']):.4g}\n"
output += "Settings\n---------\n"
output += f"n_components: {self.n_components_}\n"
output += f"max_iter: {self.max_iter_}\n"
output += f"skip_f_matrix_preprocessing: {self.skip_f_matrix_preprocessing_}\n"
return output
[docs]
def score(self, X: DataFrameDict, y: object = None, sample_weight: np.ndarray | None = None) -> float:
"""Return the :func:`~sklearn.metrics.r2_score` of the Y block on test data.
See RegressorMixin.score for the general contract.
Parameters
----------
X : DataFrameDict
Test samples. The nested ``"Y"`` block supplies the actual response
values; ``X["Z"]`` and ``X["F"]`` drive the prediction.
y : object, optional
Must be ``None``. The Y-data comes from ``X["Y"]``, not from a separate
argument (the sklearn ``RegressorMixin.score`` signature is preserved
only for API compatibility). Anything else raises, rather than being
silently ignored as it was before #565.
sample_weight : np.ndarray or None
Optional per-sample weight forwarded to
:func:`~sklearn.metrics.r2_score`.
Returns
-------
score : float
:func:`~sklearn.metrics.r2_score` of ``X["Y"]``, the one Y block a TPLS
model fits, against ``self.diagnose(X).hat``: the mean of the columns'
:math:`R^2`.
Raises
------
ValueError
If ``y`` is not ``None``, or if ``X["Y"]`` holds no blocks.
Notes
-----
Only sklearn's **default** scoring reaches this method. A named scorer *string*
never does::
cross_val_score(TPLS(...), X=DataFrameDict(blocks), cv=5) # works
cross_val_score(TPLS(...), X=DataFrameDict(blocks), cv=5,
scoring="r2") # all NaN
sklearn builds ``scoring="r2"`` into a ``_Scorer`` whose ``__call__`` requires a
``y_true`` argument. TPLS carries its response inside ``X["Y"]``, so
``cross_val_score`` receives no ``y`` to hand the scorer, and the call fails
inside sklearn before this method is reached: instrumentation shows this method
called 3 times out of 3 folds under default scoring and 0 times under
``scoring="r2"``. Because TPLS never gets control it cannot turn that into a clear
error, and sklearn's ``error_score`` (default ``np.nan``) records every fold as
``NaN`` behind a ``UserWarning``.
Use :func:`make_tpls_scorer` for any metric other than the default. sklearn passes
a *callable* ``scoring=`` through untouched and calls it as
``scorer(estimator, X_test)``, which a callable with an optional ``y`` accepts::
cross_val_score(TPLS(...), X=DataFrameDict(blocks), cv=5,
scoring=make_tpls_scorer("r2")) # works
cross_val_score(TPLS(...), X=DataFrameDict(blocks), cv=5,
scoring=make_tpls_scorer("neg_mean_squared_error"))
``make_tpls_scorer("r2")`` reproduces this method's value fold for fold, so it is
a drop-in replacement for the broken string form.
If a string scorer is used anyway, pass ``error_score="raise"`` to see the
underlying ``TypeError`` instead of a silent ``NaN``::
cross_val_score(TPLS(...), X=DataFrameDict(blocks), cv=5,
scoring="r2", error_score="raise")
# TypeError: _Scorer._score() missing 1 required positional argument: 'y_true'
See Also
--------
make_tpls_scorer : Build a ``scoring=`` callable that TPLS can honour.
"""
_reject_separate_y(y, "TPLS.score")
y_pred = self.diagnose(X).hat
y_actual = X["Y"]
if not y_actual:
msg = "y_actual is empty: X['Y'] must contain at least one block to compute a score."
raise ValueError(msg)
return float(np.mean([r2_score(y_actual[key], y_pred[key], sample_weight=sample_weight) for key in y_actual]))
[docs]
def help(self) -> str:
"""Help for the TPLS Estimator.
Data organization
-----------------
Quick tips
----------
Build model: tpls = TPLS(n_components=2, d_matrix=d_matrix).fit(X)
Get model's predictions: tpls.hat <-- the hat-matrix, i.e., the predictions
Predict on new data: tpls.diagnose(X_new)
See model summary: tpls.display_results()
This help page: tpls.help()
Statistical values
------------------
.t_scores_super Super scores for the entire model [pd.DataFrame]
.hotellings_t2 Hotelling's T2 values for each observation, per component [pd.DataFrame]
.spe Squared prediction error for each block [dict of pd.DataFrames]
.hotellings_t2_limit() Returns the Hotelling's T2 limit for the model [float]
.spe_limit[block][group]() Return the SPE limit for a group in a block; e.g. .spe_limit["Y"][group]() [float]
.vip() Variable importance (VIP) for the D- and F-blocks [dict]
"""
# This method's docstring, dedented. Removing every run of eight spaces, as this did, also
# removed the spaces that line up the columns, and ran ".spe" into its description.
return inspect.cleandoc(self.help.__doc__ or "")
def _input_data_checks(self, X: DataFrameDict) -> None:
"""Check the incoming data, and the ``d_matrix`` it is fitted with.
``fit``, the only caller, has already checked that ``X`` is a DataFrameDict, and a
DataFrameDict always holds exactly the Z, F and Y blocks, so neither is checked here.
``d_matrix`` is checked here rather than in ``__init__``, so that a value set through
``set_params`` is checked too (#1071).
"""
if not isinstance(self.d_matrix, dict):
raise TypeError(f"d_matrix must be a dict of DataFrames; got {type(self.d_matrix).__name__}.")
if not all(isinstance(df, pd.DataFrame) for df in self.d_matrix.values()):
raise TypeError("d_matrix must contain pandas DataFrames as values.")
# Group names are labels, used exactly as given. Passing them through ``str()`` in some
# places and not in others meant a non-string group name could never be fitted.
if set(X["F"]) != set(self.d_matrix):
raise ValueError("The keys in F must match the keys in D.")
# The inner PLS model regresses the first Y block alone on the super scores: a second Y block
# was deflated with the first one's loadings, without a word.
if len(X["Y"]) != 1:
msg = f"TPLS fits one Y block; got {list(X['Y'])}."
raise ValueError(msg)
# A block of blends needs two of them, as across the family: one would be taken for no variance (#1149)
for key in X["Y"]:
self._validate_df(X["Y"][key], f'Y["{key}"]', min_rows=2)
for key in X["Z"]:
self._validate_df(X["Z"][key], f'Z["{key}"]', min_rows=2)
for key in self.d_matrix:
self._validate_df(self.d_matrix[key], f'd_matrix["{key}"]')
self._validate_df(X["F"][key], f'F["{key}"]', min_rows=2)
def _learn_center_and_scaling_parameters(self, y: pd.DataFrame, label: str) -> tuple[pd.Series, pd.Series]:
"""
Learn the centering and scaling parameters for one block.
Parameters
----------
y : pd.DataFrame
The block to learn from.
label : str
Where this block sits, e.g. ``'Z["Conditions"]'``, used only in the
warnings for excluded columns.
Returns
-------
centering : pd.Series
The centering parameters.
scaling : pd.Series
The scaling parameters.
Warns
-----
UserWarning
If a column has no observed value. Its scale is set to ``NaN``, which makes the whole
column ``NaN`` and then zero, so the column is excluded from the fit (#982). A column
with no variance, whose observed values are all equal or which has one, is excluded the
same way, and is named once the fit has not refused the data
(:meth:`_name_blocks_with_no_variance`). It used to be excluded without a word (#513).
"""
# pandas 2 hands a reduction the caller's own column labels, so a name given to them after
# the fit renamed the centre and scale too; pandas 3 does not, and both keep a copy now.
y = _own_axes(y)
centering = y.mean()
# ddof=1, as pandas' std, and in any units: pandas' std was inf near 1e160 and 0 near 1e-160 (#1070)
scaling = pd.Series(column_sd(y), index=y.columns)
# An all-NaN column has a NaN mean, so it ends up excluded exactly as a constant one does;
# it is named separately, in the warning the whole family gives (#1098), because "no
# variance" would send the caller to the wrong cause.
unobserved = unobserved_columns(y, label)
# No variance: its observed values are all equal, in any units, as MCUVScaler reads a
# column. An sd below a floor (it was sqrt(eps)) also took a column in small units for a
# constant one. A block of one row, such as a group of D with one material, has no column
# that varies (#1144). An all-NaN column has a NaN maximum and minimum, which are not
# equal, so it is not counted here.
return centering, scaling.mask(_equal_values(y) | unobserved) # a column with no information has a NaN scale
def _name_blocks_with_no_variance(self, X: DataFrameDict) -> None:
"""Warn about each block that has columns with no variance, in the family's words (#1114).
Called once the fit has not refused the data, so data with no variance at all gets that
error alone, as in MBPLS. A group of D, F or Z with no variance gets zero weight (#1144),
and is named as a group. F used as given keeps a constant column in the fit, so the column
is not named (#1084); a group of it with no variation at all gets zero weight, and is.
"""
for key, block in X["Y"].items():
constant_columns(block, f'Y["{key}"]', self)
for key, block in X["Z"].items():
constant_columns(block, f'Z["{key}"]', self)
for key, block in X["F"].items():
constant_columns(self.d_matrix[key], f'D["{key}"]', self)
if not self.skip_f_matrix_preprocessing_ or self.preproc_["F"][key]["scale"].isna().all():
constant_columns(block, f'F["{key}"]', self)
def _validate_df(self, df: pd.DataFrame, label: str, min_rows: int = 1) -> None:
"""
Refuse a block that is not numbers TPLS can use: no columns, a column that is not float64 or
int64, an infinite value, or fewer than ``min_rows`` rows (`check_array` from scikit-learn, whose
other defaults hold for a DataFrame). A missing value is allowed.
Parameters
----------
df : {pd.DataFrame}
label : str
Where the block sits, as ``'Z["group"]'``, named in the message for a block with no columns.
min_rows : int, default=1
The fewest rows the block may have.
"""
refuse_empty_block(df, label)
# Ensure all columns are dtype "float64" or "int64". (``np.dtypes.IntDType`` is the
# C ``int``, which is int32, so it must not stand in for int64 here.)
if not all(good_cols := [isinstance(col, (np.dtypes.Float64DType, np.dtypes.Int64DType)) for col in df.dtypes]):
bad_columns = df.columns[[not item for item in good_cols]].to_list()
raise ValueError(
f"All columns in the DataFrame must be of type float64 or int64. Bad columns: {bad_columns}"
)
check_array(df, ensure_all_finite="allow-nan", ensure_min_samples=min_rows, estimator=self)
def _has_converged(self, starting_vector: np.ndarray, revised_vector: np.ndarray) -> bool:
"""Whether the scores converged: the norm of their change is less than ``tol`` times the norm before it.
The loop that calls this also stops after ``max_iter`` passes.
"""
# Floor ``||starting_vector||`` -- when both vectors collapse to
# zero (degenerate block / fully-deflated component) the bare
# division produced NaN, the convergence test became permanently
# False, and the loop ran to max_iter. SEC-21 (#270) sub-item 1.
delta_gap = float(
np.linalg.norm(starting_vector - revised_vector) / _nz(float(np.linalg.norm(starting_vector)))
)
return delta_gap < self.tol
@staticmethod
def _regress_f_group_on_u(f_group: np.ndarray, u: np.ndarray, present: np.ndarray) -> np.ndarray:
"""Step 2 for one F group: ``h = F'u / u'u``, or NaN where ``F'u`` is only rounding error (#1144).
A group with no variation over the blends, used as given, meets a centred ``u`` with an
``F'u`` of about 1e-17 times ``||F|| ||u||``: exactly 0 with some numpy builds, not with
others. Regressing D on that remainder gave F weights near 1e15 and took F out of the model
without a word, or refused the fit, by build. Such a group now gets zero weight in place of
this (:meth:`_zero_weight`). Judged against ``||F|| ||u||``, a group that varies, but whose
``F'u`` is as small, gives no direction on every build, and the fit refuses it.
"""
h = regress_a_space_on_b_row(f_group.T, u.T, present.T)
# The two sides are equal only where F is all 0, and then h is 0 and D cannot be regressed on it: the fit
# refuses the component either way, so < and <= give one result.
if (
np.linalg.norm(h) * np.linalg.norm(u)
< epsqrt * np.linalg.norm(f_group) # pragma: no mutate - <= differs only where F is all 0
):
return np.full_like(h, np.nan) # pragma: no mutate - a fill of None is NaN in a float array
return h
def _refuse_blocks_with_no_variance(self) -> None:
"""Refuse nothing to explain, or nothing to explain it by, before the iteration, as MBPLS does (#1144).
Scaling leaves a constant column absent, and an absent cell is 0 here, so a block with
no variance is all zeros, in any units. The iteration explains the first Y block. A group
of F carries something only with its D, as its blend scores weigh F by the properties: a
group whose F, or whose D, has no column that varies gets zero weight, as a block with no
variance does in MBPLS. Its F is judged by the scale learned from it, so that F used as
given, which is not scaled, is judged as scaled F is. With no group left, the blends have
no F score, which the model needs whether or not Z varies.
"""
y_key, y_block = next(iter(self.y_mats.items()))
check_y_varies(y_block, label=f'Y["{y_key}"]', explained_by="F, D and Z")
self._groups_without_variation = {
key
for key in self.d_mats
if self.preproc_["F"][key]["scale"].isna().all() or self.preproc_["D"][key]["scale"].isna().all()
}
if self._groups_without_variation == set(self.d_mats):
raise NotEnoughVarianceError(
"There is no variance in F and D: no group has a column of F and a column of D that vary, so the "
"blends have no F score to explain Y by. A TPLS model needs one, even when Z varies."
)
def _zero_weight(self, per_group: dict[str, np.ndarray]) -> dict[str, np.ndarray]:
"""Give each group with no variation a zero vector in place of its regression: zero weight (#1144).
Regressed as the others are, such a group gave 0 / 0 (its F or D all absent, so all 0) or
rounding noise (F used as given, whose ``F'u`` is about 1e-17 times its size), and either
one stopped the fit, or took the other groups out of the model, by numpy build.
"""
return {
key: np.zeros_like(vector) if key in self._groups_without_variation else vector
for key, vector in per_group.items()
}
def _refuse_scores_that_are_not_finite(self, component: int, t_super_i: np.ndarray) -> None:
"""Raise ``NotEnoughVarianceError`` when a component's scores are not finite.
A NaN never leaves the iteration once it is in it, and the convergence test can never
pass on it, so the fit would run to ``max_iter`` and return NaN everywhere (#709). A
row or a material with nothing observed no longer causes one: it gets a zero score or
loading and no weight (#513), and nor does a group whose F or D has no variation: it gets
zero weight (#1144). What is left is a direction the iteration cannot take, because the
vector it regresses on is zero: a group of F that does not covary with Y, or a block that
the earlier components have used up. It was a ``ValueError``; it is the family's error
for data without the variance to fit, as the checks before the iteration are (#1144).
"""
if np.isfinite(t_super_i).all():
return
msg = (
f"TPLS cannot fit component {component}: its scores are not finite, because a block has "
"no variation for the iteration to follow. Check that each group of F covaries with the "
"first Y block, or fit fewer components."
)
raise NotEnoughVarianceError(msg)
def _super_block_names(self) -> list[str]:
"""Row labels of ``w_loadings_super``: one per Z group, as the fit scores each on its own, then F."""
if self.n_conditions == 0:
return ["F"]
# A single Z group keeps the plain "Z" label it always had.
return ["Z", "F"] if len(self.z_mats) == 1 else [*(f'Z["{key}"]' for key in self.z_mats), "F"]
def _store_model_coefficients(
self,
pc_a_column: int, # one-based index for the component
t_super_i: np.ndarray,
r_i: dict[str, np.ndarray],
w_i_z: dict[str, np.ndarray],
w_super_i: np.ndarray,
s_i: dict[str, np.ndarray],
) -> None:
"""Store the model coefficients for later use."""
self.t_scores_super = self.t_scores_super.join(
pd.DataFrame(t_super_i, index=self.observation_names, columns=[pc_a_column])
)
# These are loadings really, not scores, for each group in the F block.
self.r_loadings_f = {
key: self.r_loadings_f[key].join(
pd.DataFrame(r_i[key], index=self.material_names[key], columns=[pc_a_column])
)
for key in r_i
}
# These are the loadings for the Z space
self.w_loadings_z = {
key: self.w_loadings_z[key].join(
pd.DataFrame(w_i_z[key], index=self.condition_names[key], columns=[pc_a_column])
)
for key in w_i_z
}
self.w_loadings_super = self.w_loadings_super.join(
pd.DataFrame(w_super_i, index=self._super_block_names(), columns=[pc_a_column])
)
self.s_loadings_d = {
key: self.s_loadings_d[key].join(
pd.DataFrame(s_i[key], index=self.property_names[key], columns=[pc_a_column])
)
for key in s_i
}
def _calculate_and_store_deflation_matrices(
self,
pc_a: int,
t_super_i: np.ndarray,
q_super_i: np.ndarray,
r_i: dict[str, np.ndarray],
) -> None:
"""
Calculate and store the deflation matrices for the TPLS model.
Deflate the matrices stored in the instance object.
Returns the prediction matrices in a dictionary.
"""
# Step 13: Deflate the Z matrix with a loadings vector, pz_b (_b is for block)
pz_b = {
key: regress_a_space_on_b_row(df_z.T, t_super_i.T, self.not_na_z[key].T)
for key, df_z in self.z_mats.items()
}
for key in self.z_mats:
self.z_mats[key] -= (t_super_i @ pz_b[key].T) * self.not_na_z[key]
self.p_loadings_z = {
key: self.p_loadings_z[key].join(pd.DataFrame(pz_b[key], index=self.condition_names[key], columns=[pc_a]))
for key in pz_b
}
# Step 13. p_i = F_i' t_i / t_i't_i. Regress the columns of F_i on t_i; store slope coeff in vectors p_i.
# Note: the "t" vector is the t_i vector from the inner PLS model, marked as "Tt" in figure 4 of the paper.
# It is the score column from the super score matrix regression onto Y.
# A group with no variation takes no part in the model: its F and D are not deflated (#1144)
pf_i = self._zero_weight(
{
key: regress_a_space_on_b_row(df_f.T, t_super_i.T, self.not_na_f[key].T)
for key, df_f in self.f_mats.items()
}
)
self.p_loadings_f = {
key: self.p_loadings_f[key].join(pd.DataFrame(pf_i[key], index=self.material_names[key], columns=[pc_a]))
for key in pf_i
}
# Step 13: v_i = D_i' r_i / r_i'r_i. Regress the rows of D_i (properties) on r_i; store slopes in v_i.
v_i = self._zero_weight(
{
key: regress_a_space_on_b_row(df_d.T, r_i[key].T, self.not_na_d[key].T)
for key, df_d in self.d_mats.items()
}
)
self.v_loadings_d = {
key: self.v_loadings_d[key].join(pd.DataFrame(v_i[key], index=self.property_names[key], columns=[pc_a]))
for key in v_i
}
# Step 14. Do the actual deflation.
for key in self.d_mats:
# Step to deflate F matrix
self.f_mats[key] -= (t_super_i @ pf_i[key].T) * self.not_na_f[key]
# Two sets of matrices to deflate: properties D and formulas F.
self.d_mats[key] -= (r_i[key] @ v_i[key].T) * self.not_na_d[key]
# Deflate the Y-space as well
self.q_loadings_y = {
key: self.q_loadings_y[key].join(pd.DataFrame(q_super_i, index=self.quality_names[key], columns=[pc_a]))
for key in self.y_mats
}
for key in self.y_mats:
self.hat_[key] += t_super_i @ q_super_i.T
self.y_mats[key] -= (t_super_i @ q_super_i.T) * self.not_na_y[key]
def _update_performance_statistics(self) -> None:
"""Calculate and store the performance statistics of the model, such as SSQ, R2, etc."""
# Calculate the sums of squares for each block, per column.
# Note: the `ddof=0` is used to calculate the population variance, which is proportional to the SSQ.
calc_ssq = {
"D": {key: np.nansum(self.d_mats[key] ** 2, axis=0) for key in self.d_mats},
"F": {key: np.nansum(self.f_mats[key] ** 2, axis=0) for key in self.f_mats},
"Z": {key: np.nansum(self.z_mats[key] ** 2, axis=0) for key in self.z_mats},
"Y": {key: np.nansum(self.y_mats[key] ** 2, axis=0) for key in self.y_mats},
}
self.sums_of_squares_.append(calc_ssq)
# Calculate the incremental (not cumulative!) R2 values for each block, per column:
# Cumulative R2 values can be found by summation. The R2 values are **always** fractional (between 0 and 1).
ssq_prior_pc = self.sums_of_squares_[-2]
ssq_start_0 = self.sums_of_squares_[0]
# A column left out for having no variance has a sum of squares of 0, and an R2 of NaN.
with np.errstate(invalid="ignore"):
calc_r2 = {
"D": {
key: (ssq_prior_pc["D"][key] - calc_ssq["D"][key]) / ssq_start_0["D"][key] for key in self.d_mats
},
"F": {
key: (ssq_prior_pc["F"][key] - calc_ssq["F"][key]) / ssq_start_0["F"][key] for key in self.f_mats
},
"Z": {
key: (ssq_prior_pc["Z"][key] - calc_ssq["Z"][key]) / ssq_start_0["Z"][key] for key in self.z_mats
},
"Y": {
key: (ssq_prior_pc["Y"][key] - calc_ssq["Y"][key]) / ssq_start_0["Y"][key] for key in self.y_mats
},
}
self.r2_frac.append(calc_r2)
# VIP for each block, for given number of components we currently have. VIP are cumulative.
# For the D-block and F-block: you want to be able to use the VIPs to compare across the entire block, without
# regards to the banding in groups that might have to be done. So use the total R2 for that block, and do not
# work per group.
# Default ("vip") feature importance: standard VIP computed on the raw block loadings. This matches
# the variable importance reported by Garcia-Munoz (2014) (e.g. Table 4). The opt-in "deflated"
# variant, exposed through `vip(method="deflated")`, instead uses the direct (rotated) weights
# S(V^T S)^-1 and P(P^T P)^-1; see `_feature_importance`.
loadings_s = np.concatenate(list(self.s_loadings_d.values()))
loadings_f = np.concatenate(list(self.p_loadings_f.values()))
vip_split_d, vip_split_f = self._split_block_vip(loadings_s, loadings_f)
self.feature_importance["D"] = vip_split_d
self.feature_importance["F"] = vip_split_f
def _split_block_vip(
self, loadings_d: np.ndarray, loadings_f: np.ndarray
) -> tuple[dict[str, pd.Series], dict[str, pd.Series]]:
"""Compute VIP for the D- and F-blocks from (n_features x A) loading/weight matrices.
The per-component R2 weighting is taken from `self.r2_frac` (the first entry is the pre-fit
baseline and is skipped). The flat VIP vector for each block is split back into the original
per-group :class:`pandas.Series`, indexed by feature name. Shared by the default VIP (raw block
loadings) and the opt-in deflated direct-weights variant.
"""
r2_d_a = [_mean_over_groups(r2_frac["D"]) for r2_frac in self.r2_frac]
r2_f_a = [_mean_over_groups(r2_frac["F"]) for r2_frac in self.r2_frac]
vip_d = self._calculate_vip(loadings_d, np.array(r2_d_a[1:]), self._informative_width("D"))
vip_f = self._calculate_vip(loadings_f, np.array(r2_f_a[1:]), self._informative_width("F"))
# Split each flat VIP vector back into the original groups it was concatenated from.
# For example: with two groups of 17 and 3 columns, `vip_d` has 20 values; the first 17 belong to
# the first group and the last 3 to the second.
vip_split_d = {}
start = 0
for key in self.property_names:
end = start + len(self.property_names[key])
vip_split_d[key] = pd.Series(vip_d[start:end], index=self.s_loadings_d[key].index)
start = end
vip_split_f = {}
start = 0
for key in self.material_names:
end = start + len(self.material_names[key])
vip_split_f[key] = pd.Series(vip_f[start:end], index=self.r_loadings_f[key].index)
start = end
return vip_split_d, vip_split_f
def _feature_importance(self, method: str) -> dict[str, dict[str, pd.Series]]:
"""Return the D-/F-block feature importance for the requested `method`.
Parameters
----------
method : {"vip", "deflated"}
``"vip"`` returns the standard VIP stored at fit time (computed from the raw block loadings,
matching Garcia-Munoz, 2014). ``"deflated"`` computes VIP on the direct (rotated) weights that
map the *original* variables onto the scores while accounting for the deflation across
components: ``S(V^T S)^-1`` for the D-block and ``P(P^T P)^-1`` for the F-block, where ``S`` and
``V`` are the D-block correlation- and deflation-loadings and ``P`` is the F-block deflation
loading (equation 7 in the paper).
"""
if method == "vip":
return self.feature_importance
s_mat = np.concatenate(list(self.s_loadings_d.values())) # (K_D x A)
v_mat = np.concatenate(list(self.v_loadings_d.values())) # (K_D x A)
p_mat = np.concatenate(list(self.p_loadings_f.values())) # (K_F x A)
# Direct (rotated) weights; `pinv` keeps this well-defined if the A x A core is near-singular.
w_d = s_mat @ np.linalg.pinv(v_mat.T @ s_mat)
w_f = p_mat @ np.linalg.pinv(p_mat.T @ p_mat)
split_d, split_f = self._split_block_vip(w_d, w_f)
return {"D": split_d, "F": split_f, "Z": {}, "Y": {}}
[docs]
def vip(
self, block: str | None = None, method: str = "vip"
) -> dict[str, dict[str, pd.Series]] | dict[str, pd.Series]:
"""Return Variable Importance in Projection (VIP) scores for TPLS blocks.
VIP scores are computed during fitting for the D-block (material properties) and
F-block (formulation variables) and stored in :attr:`feature_importance`. Each
component is weighted by its share of that block's own variance (its per-column
R2, averaged over the columns of each group and then over the groups with a column
that varies), not of Y as a PLS VIP is. A group with no variation has zero weight, and a
VIP of 0.
Parameters
----------
block : str or None, default=None
Which block to return. Must be ``"D"`` or ``"F"``, or ``None`` to
return all blocks.
method : {"vip", "deflated"}, default="vip"
How importance is measured.
- ``"vip"`` (default): standard VIP on the raw block loadings, as reported by
Garcia-Munoz (2014). These are the values stored in :attr:`feature_importance`.
- ``"deflated"``: VIP computed on the direct (rotated) weights that map the *original*
variables onto the scores while accounting for the deflation across components,
``S(V^T S)^-1`` for the D-block and ``P(P^T P)^-1`` for the F-block. This is an
alternative importance measure; it does not change :attr:`feature_importance`.
Returns
-------
dict
If *block* is ``None``: a four-key dict keyed by block, of the shape
``{"D": {group: pd.Series, ...}, "F": {group: pd.Series, ...},
"Z": {}, "Y": {}}``. VIP is defined only for the ``D``- and
``F``-blocks, so the ``"Z"`` and ``"Y"`` entries are always empty
dicts; they are kept for shape-stability across every ``vip()``
caller. If *block* is ``"D"`` or ``"F"``: the inner dict
``{group: pd.Series, ...}`` for that block, where each
``pd.Series`` is indexed by feature names.
Raises
------
ValueError
If the model is not fitted, *block* is not ``"D"``, ``"F"``, or ``None``, or *method*
is not ``"vip"`` or ``"deflated"``.
Examples
--------
>>> tpls = TPLS(...).fit(data)
>>> tpls.vip() # all blocks, standard VIP
>>> tpls.vip("D") # D-block only → {group_name: pd.Series, ...}
>>> tpls.vip("D", method="deflated") # D-block deflated direct-weights importance
"""
check_is_fitted(self)
if method not in ("vip", "deflated"):
msg = f"method must be 'vip' or 'deflated'; got {method!r}."
raise ValueError(msg)
if block is not None and block not in ("D", "F"):
msg = f"block must be 'D', 'F', or None; got {block!r}."
raise ValueError(msg)
importance = self._feature_importance(method)
if block is None:
return importance
return importance[block]
@property
def d_block_scaling_(self) -> dict[str, float]:
"""Block-scaling factor applied to each D-block (read-only).
After column-wise centring and auto-scaling, every D-block ``X_i`` is additionally divided by
``sqrt(P_i * M_i)`` (``P_i`` = number of lots/rows, ``M_i`` = number of properties/columns) so that
``trace(X_i^T X_i) ~= 1``, removing bias toward blocks with more lots or properties
(Garcia-Munoz, 2014, section 2.1).
Returns
-------
dict[str, float]
Mapping of D-block group name to its scalar block-scaling factor.
Raises
------
AttributeError
If the model has not been fitted yet.
"""
check_is_fitted(self)
return {key: self.preproc_["D"][key]["block"] for key in self.preproc_["D"]}
def _informative_width(self, block: str) -> int | None:
"""Return the count of columns of D or F that vary, in the groups with weight: a VIP's width (#913).
A column with no variance, or one of a group with zero weight (#1144), has a loading of 0, as in
MBPCA and MBPLS, so counted in the width it raised the VIP of every other column, and the fit's VIP
changed although nothing it used had. ``None`` for a model fitted before groups had zero weight:
it keeps the full width it was fitted with.
"""
if not hasattr(self, "_groups_without_variation"):
return None
return sum(
int(self.preproc_[block][key]["scale"].notna().sum())
for key in self.d_mats
if key not in self._groups_without_variation
)
@staticmethod
def _calculate_vip(loadings: np.ndarray, r2_vector: np.ndarray, width: int | None) -> np.ndarray:
"""Calculate the VIP values for the current component.
The `loadings` has as many rows as there are feature varaibles, and A columns, where A = number of components.
The `r2_vector` is a vector of fractional R^2 values for the current component, with `A` entries.
The `r2_vector` values should be between 0 and 1; the fraction of variance explained by the component for that
given `loadings` matrix.
The VIP values are calculated as follows:
VIP = sqrt(n * sum((r2_vector * (loadings ** 2)) / sum(r2_vector)))
where n is ``width``, the number of columns that vary (:meth:`_informative_width`), or, when it is
``None``, the number of rows in the loadings matrix.
"""
n = loadings.shape[0] if width is None else width
return np.sqrt(n * np.sum(r2_vector * (loadings**2), axis=1) / np.sum(r2_vector))
def _cumulative_hotellings_t2(self, t_scores_super: np.ndarray) -> np.ndarray:
"""Cumulative Hotelling's T2 per component from super-scores.
Matches the PCA / PLS convention (#502): column ``a`` holds
``sum_{j<=a} (t_j / s_j)^2`` where ``s_j`` is the training score
standard deviation (unbiased, N - 1 divisor) stored in
``scaling_factor_for_scores``. A component whose score standard
deviation is at most epsqrt times the largest (a rank-deficient fit)
carries no T2 information and adds 0, by the rule every other T2
follows (:func:`~process_improve.multivariate._diagnostics._t2_score_sd`).
Parameters
----------
t_scores_super : np.ndarray of shape (n_obs, n_components)
Super-scores of the observations to evaluate.
Returns
-------
np.ndarray of shape (n_obs, n_components)
Cumulative Hotelling's T2 after each component.
"""
score_sd = _t2_score_sd(self.scaling_factor_for_scores.to_numpy())
return np.cumsum(t_scores_super**2 / score_sd**2, axis=1)
def _calculate_model_statistics_and_limits(self) -> None:
"""Calculate and store the model limits.
Limits calculated:
1. Hotelling's T2 limits
2. Squared prediction error limits
Other calculations:
1. The model's Y-space predictions are scaled back to the original space.
"""
# Hotelling's T2, cumulative per component, in the same form as PCA / PLS
# and as ``diagnose`` (#502): column ``a`` holds ``sum_{j<=a} (t_j / s_j)^2``.
# The score standard deviations use the unbiased (N - 1) divisor, which is the
# estimator the F-distribution limit in ``hotellings_t2_limit`` is derived for;
# the previous N divisor inflated T2 by N / (N - 1). N is at least 2: one blend's
# Y column centres to 0, so the first regression on it is 0 / 0 and fit refuses it.
scores_np = self.t_scores_super.to_numpy()
n_obs = scores_np.shape[0]
self.scaling_factor_for_scores = pd.Series(
np.sqrt(np.diag(scores_np.T @ scores_np) / (n_obs - 1)),
index=[a + 1 for a in range(self.n_components_)],
name="Standard deviation per score",
)
self.hotellings_t2 = pd.DataFrame(
self._cumulative_hotellings_t2(scores_np),
index=self.observation_names,
columns=list(range(1, self.n_components_ + 1)),
)
# Squared prediction error limits. This is a measure of the prediction error = difference between the actual
# and predicted values. Since the matrices are deflated by the predictive part of the model already, the
# data in these matrices is already the prediction error. Calculate the **squared** portion, and store it.
column_name = [f"SPE with A={self.n_components_}"]
def row_spe(residuals: np.ndarray, present: np.ndarray) -> np.ndarray:
# A row with nothing observed in the block has no SPE there, as in MBPCA and MBPLS (#1114). A block
# the model uses no column of, as a group with no variation, has none left in any row: 0, as MBPLS
# gives a block with no variance, and as diagnose gives it (#1144).
return np.where(present.any(axis=1) | ~present.any(), np.sqrt(np.sum(np.square(residuals), axis=1)), np.nan)
self.spe["Y"] = {
key: pd.DataFrame(
row_spe(self.y_mats[key], self.not_na_y[key]), index=self.observation_names, columns=column_name
)
for key in self.y_mats
}
self.spe["Z"] = {
key: pd.DataFrame(
row_spe(self.z_mats[key], self.not_na_z[key]), index=self.observation_names, columns=column_name
)
for key in self.z_mats
}
# SPE for the D-space, one per material: the rows of D.
self.spe["D"] = {
key: pd.Series(row_spe(self.d_mats[key], self.not_na_d[key]), index=self.r_loadings_f[key].index)
for key in self.d_mats
}
self.spe["F"] = {
key: pd.DataFrame(
row_spe(self.f_mats[key], self.not_na_f[key]), index=self.observation_names, columns=column_name
)
for key in self.f_mats
}
# Each limit over the rows that have an SPE
for block in ("Y", "Z", "D", "F"):
self.spe_limit[block] = {
key: partial(spe_calculation, spe.dropna().to_numpy()) for key, spe in self.spe[block].items()
}
# The Y-space predictions, scaled back from the preprocessed space to the original one.
self.hat: dict[str, pd.DataFrame] = {
key: self.hat_[key] * self.preproc_["Y"][key]["scale"].to_numpy()
+ self.preproc_["Y"][key]["center"].to_numpy()
for key in self.y_mats
}
def _preprocess_data(self) -> None:
"""Pre-process the training data."""
for key in self.f_mats:
if not self.skip_f_matrix_preprocessing_:
self.f_mats[key] = (
self.f_mats[key] - self.preproc_["F"][key]["center"].to_numpy()[None, :]
) / self.preproc_["F"][key]["scale"].to_numpy()[None, :]
self.d_mats[key] = (
(self.d_mats[key] - self.preproc_["D"][key]["center"].to_numpy()[None, :])
/ self.preproc_["D"][key]["scale"].to_numpy()[None, :]
/ self.preproc_["D"][key]["block"] # scalar block-scaling factor
)
for key in self.z_mats:
self.z_mats[key] = (
self.z_mats[key] - self.preproc_["Z"][key]["center"].to_numpy()[None, :]
) / self.preproc_["Z"][key]["scale"].to_numpy()[None, :]
for key in self.y_mats:
self.y_mats[key] = (
self.y_mats[key] - self.preproc_["Y"][key]["center"].to_numpy()[None, :]
) / self.preproc_["Y"][key]["scale"].to_numpy()[None, :]
def _fit_iterative_regressions(self) -> None:
"""Fit the model via iterative regressions and store the model coefficients in the class instance."""
# Formula matrix: assemble all not-na maps from blocks in F: make a single matrix.
pmap_f = np.concatenate(list(self.not_na_f.values()), axis=1)
# No variance left in the Y block the iteration explains is judged against its size before the first
# component, as PLS and MBPLS judge it: a Y fitted exactly is refused in their words (#1128, #1149).
y_size = float(np.sqrt(np.sum(next(iter(self.y_mats.values())) ** 2)))
# Follow the steps in the paper on page 54
for pc_a in range(self.n_components_):
check_variance_left("Y", next(iter(self.y_mats.values())), pc_a, y_size)
milliseconds_start = time.time()
# Step 1: Start from the Y column with the largest sum of squares, as PLS does (#195). The first
# column can be constant, all zero after preprocessing, and then the iteration has nothing to follow (#954).
u_super_i = largest_ssq_column(next(iter(self.y_mats.values())))
# The first pass always runs. A test before it, against a previous u primed as u + 1, ended the
# loop before it began once tol was above about 0.7, and the fit failed for want of a score.
for n_iter in range(1, self.max_iter_ + 1): # noqa: B007 - the statistics after the loop read it
u_prior = u_super_i # rebound below, never changed in place
# Step 2. h_i = F_i' u / u'u. Regress the columns of F on u_i, and store the slope coeff in vectors h_i
h_i = self._zero_weight(
{
key: self._regress_f_group_on_u(df_f, u_super_i, self.not_na_f[key])
for key, df_f in self.f_mats.items()
}
)
# A NaN here makes every later step and the scores NaN, through numpy's warnings: refuse it now.
self._refuse_scores_that_are_not_finite(pc_a + 1, np.concatenate(list(h_i.values())))
# Step 3. s_i = D_i' h_i / h_i'h_i. Regress the rows of D_i on h_i, and store slope coeff in vectors s_i
s_i = self._zero_weight(
{
key: regress_a_space_on_b_row(df_d.T, h_i[key].T, self.not_na_d[key].T)
for key, df_d in self.d_mats.items()
}
)
# Step 4: combine the entries in s_i to form a joint `s` and normalize it to unit length.
joint_s_normalized = np.linalg.norm(np.concatenate(list(s_i.values())))
s_i = {key: s / joint_s_normalized for key, s in s_i.items()}
# Step 5: r_i = D_i s_i / s_i's_i. Regress columns of D_i on s_i, and store slope coefficients in r_i.
r_i = self._zero_weight(
{
key: regress_a_space_on_b_row(df_d, s_i[key].T, self.not_na_d[key])
for key, df_d in self.d_mats.items()
}
)
# Step 6: Combine the entries in r_i to form a joint r (which is the name of the method in the paper).
# Horizontally concatenate all matrices in F_i to form a joint F matrix.
# Regress rows of the joint F matrix onto the joint r vector. Store coeff in block scores, t_f
joint_r = np.concatenate(list(r_i.values()))
joint_f = np.concatenate(list(self.f_mats.values()), axis=1)
t_f = regress_a_space_on_b_row(joint_f, joint_r.T, pmap_f)
# If there is a Condition matrix (non-empty Z block)
if self.n_conditions > 0:
# Step 7: w_i = Z_i' u / u'u. Regress the columns of Z on u_i, and store the slope coefficients
# in vectors w_i.
w_i_z = {
key: regress_a_space_on_b_row(df_z.T, u_super_i.T, self.not_na_z[key].T)
for key, df_z in self.z_mats.items()
}
# Step 8: Normalize joint w to unit length. See MB-PLS by Westerhuis et al. 1998. This is normal.
# Floor each per-block norm so a degenerate Z block doesn't yield NaN. SEC-21 (#270) sub-item 1.
w_i_z = {key: w / _nz(float(np.linalg.norm(w))) for key, w in w_i_z.items()}
# Step 9: regress rows of Z on w_i, and store slope coefficients in t_z. There is an error in the
# paper here, but in figure 4 it is clear what should be happening.
t_zb = {
key: regress_a_space_on_b_row(df_z, w_i_z[key].T, self.not_na_z[key])
for key, df_z in self.z_mats.items()
}
t_z = np.hstack(list(t_zb.values()))
else:
# Step 7: No Z block. Take an empty matrix across to the the superblock.
w_i_z = {}
t_z = np.zeros((t_f.shape[0], 0)) # empty matrix: in other words, no Z block
# Step 10: Combine t_z and t_f to form a joint t matrix.
t_combined = np.concatenate([t_z, t_f], axis=1)
# Step 11: Build an inner PLS model: using the t_combined as the X matrix, and the Y (quality space)
# as the Y matrix.
# Every block score is present. A mask of float ones would weigh the cells as these booleans do.
complete = np.ones(t_combined.shape, dtype=bool) # pragma: no mutate - float ones weigh the same
inner_pls = internal_pls_nipals_fit_one_pc(
x_space=t_combined,
y_space=np.array(next(iter(self.y_mats.values()))),
x_present_map=complete,
y_present_map=np.array(next(iter(self.not_na_y.values()))),
)
u_super_i = inner_pls["u_i"] # steps 2 and 7 regress F and Z on it, and it is the convergence check
t_super_i = inner_pls["t_i"]
q_super_i = inner_pls["q_i"]
w_super_i = inner_pls["w_i"]
self._refuse_scores_that_are_not_finite(pc_a + 1, t_super_i)
# Read after the loop: the last pass allowed can converge too (#1080).
converged = self._has_converged(starting_vector=u_prior, revised_vector=u_super_i)
if converged:
break
logger.debug(
"TPLS: super-component %d %s in %d iterations",
pc_a + 1,
"converged" if converged else "did not converge",
n_iter,
)
# After convergance. Step 12: Now store information.
# =================
# Floor ``||u_prior||`` -- see SEC-21 (#270) sub-item 1.
delta_gap = float(np.linalg.norm(u_prior - u_super_i) / _nz(float(np.linalg.norm(u_prior))))
self.fitting_statistics["iterations"].append(n_iter)
self.fitting_statistics["convergance_tolerance"].append(delta_gap)
self.fitting_statistics["milliseconds"].append((time.time() - milliseconds_start) * 1000)
self.fitting_statistics["converged"].append(converged)
# Store model coefficients
self._store_model_coefficients(
pc_a + 1, t_super_i=t_super_i, r_i=r_i, w_i_z=w_i_z, w_super_i=w_super_i, s_i=s_i
)
# Calculate and store the deflation vectors. See equation 7 on page 55.
self._calculate_and_store_deflation_matrices(pc_a + 1, t_super_i=t_super_i, q_super_i=q_super_i, r_i=r_i)
# Update performance statistics for this component
self._update_performance_statistics()
# Step 15: Calculate the final model limits (after all components have been fitted).
self._calculate_model_statistics_and_limits()
self._warn_about_components_not_converged()
def _warn_about_components_not_converged(self) -> None:
"""Warn about each component whose loop stopped at ``max_iter``, in the words of MBPLS and MBPCA.
The loop is not a power iteration, and it can alternate between two states for ever: the
component is then the state of the last pass, which differs from one platform to another (#1163).
"""
failed = [a + 1 for a, converged in enumerate(self.fitting_statistics["converged"]) if not converged]
if failed:
warnings.warn(
f"TPLS 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(),
)
# Named regression metrics understood by :func:`make_tpls_scorer`, each paired with the
# sign that turns it into a "higher is better" score. The signs and names mirror
# sklearn's own ``neg_*`` convention, so a TPLS scorer is read the same way as the string
# scorer it stands in for.
_TPLS_METRICS: dict[str, tuple[Callable[..., float], int]] = {
"r2": (r2_score, 1),
"neg_mean_absolute_error": (mean_absolute_error, -1),
"neg_mean_squared_error": (mean_squared_error, -1),
"neg_root_mean_squared_error": (root_mean_squared_error, -1),
}
def _mean_over_groups(r2_per_column: dict[str, np.ndarray]) -> float:
"""Return a block's R² for one component: the mean over its groups of each group's mean over its columns.
A group with no column that varies has no R² (0 / 0, NaN) and is left out. It made the
block's R², and so the VIP of every group of the block, NaN (#1169). Its own VIP is 0, as
in MBPCA and MBPLS: its loadings are 0, as it has zero weight (#1144).
"""
return float(np.mean([np.nanmean(r2) for r2 in r2_per_column.values() if not np.isnan(r2).all()]))
def _f_in_the_order_of_d(f_group: pd.DataFrame, d_group: pd.DataFrame, group: object) -> pd.DataFrame:
"""Return F's columns in the order of D's rows, each material paired with its properties by name (#1197).
They were paired by position, so F with the materials of a group in another order than D fitted
each material with another one's properties, without a word. A material that only one of them
has, or that D names twice, has no pair, and is refused by name.
Parameters
----------
f_group : pd.DataFrame
One group of F: blends as rows, materials as columns.
d_group : pd.DataFrame
The same group of D: materials as rows, properties as columns.
group : object
The group's name, used in the messages.
Returns
-------
pd.DataFrame
``f_group`` with its columns in the order of ``d_group``'s rows.
Raises
------
ValueError
If D names a material twice, or a material of F or of D is not in the other.
"""
materials = d_group.index
if materials.has_duplicates:
twice = materials[materials.duplicated()].unique().to_list()
msg = f'd_matrix["{group}"] names a material more than once: {twice}.'
raise ValueError(msg)
if unknown := f_group.columns.difference(materials, sort=False).to_list():
msg = (
f'F["{group}"] has materials that d_matrix["{group}"] does not have: {unknown}. '
"Each column of F is paired with the row of d_matrix of the same name."
)
raise ValueError(msg)
if unused := materials.difference(f_group.columns, sort=False).to_list():
msg = (
f'd_matrix["{group}"] has materials that F["{group}"] has no column for: {unused}. '
"Give a material that no blend uses a column of zeros in F."
)
raise ValueError(msg)
return f_group[materials]
def _reject_separate_y(y: object, where: str) -> None:
"""Raise if a caller supplies a conventional ``y`` to a T-shaped model.
TPLS reads its response out of ``X["Y"]``, so a separate ``y`` is never the target
it scores against. Silently ignoring it (the behaviour before #565) let callers
believe a response had been supplied when it had not.
Parameters
----------
y : object
The value passed as ``y``. Anything other than ``None`` is an error.
where : str
Qualified name of the caller, used to open the error message.
Raises
------
ValueError
If ``y`` is not ``None``.
"""
if y is None:
return
msg = (
f"{where} does not take a separate `y`: a T-shaped model carries its response "
'inside X["Y"]. Put the response in the "Y" block of the DataFrameDict and leave '
"`y` as None. For cross-validation with a named metric, pass "
"scoring=make_tpls_scorer('r2') rather than scoring='r2'."
)
raise ValueError(msg)
[docs]
def make_tpls_scorer(
metric: str | Callable[..., float] = "r2",
*,
greater_is_better: bool = True,
**metric_kwargs: object,
) -> Callable[..., float]:
"""Build a ``scoring=`` callable that works with :class:`TPLS` (#565).
A *named* scorer string cannot be used with TPLS. sklearn turns ``scoring="r2"`` into
a ``_Scorer`` whose ``__call__`` takes ``y_true`` as a required positional argument,
but a T-shaped model carries its response inside ``X["Y"]`` rather than in a separate
``y``, so :func:`~sklearn.model_selection.cross_val_score` has no ``y`` to hand over
and calls the scorer as ``scorer(estimator, X_test)``. The call then fails on the
missing argument *before TPLS is reached*, and sklearn's ``error_score`` records every
fold as ``NaN`` behind a ``UserWarning``.
sklearn passes a *callable* ``scoring=`` through untouched and invokes it the same
way, so a callable whose ``y`` is optional receives that two-argument call cleanly and
can read the response out of ``X["Y"]`` itself. That is what this factory returns.
Parameters
----------
metric : str or Callable[..., float], optional
Either one of ``"r2"``, ``"neg_mean_absolute_error"``, ``"neg_mean_squared_error"``
or ``"neg_root_mean_squared_error"``, whose sign convention matches sklearn's, or
a callable ``metric(y_true, y_pred, **kwargs)``
returning a float. Default is ``"r2"``, which reproduces :meth:`TPLS.score`.
greater_is_better : bool, optional
Consulted only when ``metric`` is a callable: ``False`` flips the sign so the
result still obeys sklearn's "higher is better" contract. The named metrics carry
their own sign and ignore this. Default is True.
**metric_kwargs : object
Extra keyword arguments forwarded to ``metric`` on every call.
Returns
-------
scorer : Callable[..., float]
A callable with signature ``scorer(estimator, X, y=None, sample_weight=None)``,
usable as the ``scoring=`` argument of
:func:`~sklearn.model_selection.cross_val_score`,
:func:`~sklearn.model_selection.cross_validate` and the ``*SearchCV`` classes.
Several Y blocks are combined the way :meth:`TPLS.score` combines them: an
unweighted mean over the blocks in ``X["Y"]``, not a pooled multiblock statistic.
Raises
------
ValueError
If ``metric`` is a string that is not a known metric name. The returned scorer
raises :class:`ValueError` in turn if it is called with a non-``None`` ``y``, or
with an ``X`` whose ``"Y"`` block is empty.
See Also
--------
TPLS.score : The default scoring path, equivalent to ``make_tpls_scorer("r2")``.
Examples
--------
>>> from sklearn.model_selection import cross_val_score
>>> scorer = make_tpls_scorer("neg_root_mean_squared_error")
>>> model = TPLS(n_components=2, d_matrix=d_matrix) # doctest: +SKIP
>>> cross_val_score(model, X=blocks, cv=5, scoring=scorer) # doctest: +SKIP
array([-1.07, -0.98, -1.11, -1.02, -1.05])
"""
if callable(metric):
score_func: Callable[..., float] = metric
sign = 1 if greater_is_better else -1
name = getattr(metric, "__name__", "metric")
else:
try:
score_func, sign = _TPLS_METRICS[metric]
except KeyError:
msg = (
f"Unknown metric {metric!r} for make_tpls_scorer. Choose one of "
f"{sorted(_TPLS_METRICS)}, or pass a callable metric(y_true, y_pred)."
)
raise ValueError(msg) from None
name = metric
def tpls_scorer(
estimator: TPLS,
X: DataFrameDict,
y: object = None,
sample_weight: np.ndarray | None = None,
) -> float:
"""Score a fitted :class:`TPLS` model against the response inside ``X["Y"]``."""
_reject_separate_y(y, "A TPLS scorer")
check_is_fitted(estimator)
y_actual = X["Y"]
if not y_actual:
msg = 'X["Y"] must contain at least one block to score a TPLS model.'
raise ValueError(msg)
y_pred = estimator.diagnose(X).hat
per_block = [
score_func(y_actual[key], y_pred[key], sample_weight=sample_weight, **metric_kwargs) for key in y_actual
]
return sign * float(np.mean(per_block))
tpls_scorer.__name__ = f"tpls_scorer({name})"
tpls_scorer.__qualname__ = tpls_scorer.__name__
return tpls_scorer
