Multivariate Analysis#

projection_matrix(observed, *, method='tsr', ridge=0.0)[source]#

Build the fixed linear operator mapping observed columns to score estimates.

For a fixed missingness pattern, every estimator in project() is a fixed linear map t_hat = M @ z_observed. This method exposes that matrix so callers that reuse one pattern many times (an online monitor at time sample k, or an optimiser treating candidate columns as observed) can precompute it once.

Parameters:

observed (array-like) – Either a boolean mask of length n_features_in_ (True = observed), or a list of feature labels to treat as observed.
method ({"tsr", "scp", "pmp"}, default="tsr")
ridge (float, default=0.0)

Returns:

result – With keys matrix (DataFrame, n_components x n_observed, columns labelled by the observed features), condition_number (float) and method.

Return type:

score(X, y=None)[source]#

Negative mean squared reconstruction error (higher is better).

Follows the sklearn convention where higher scores indicate better model fit. This makes PCA compatible with cross_val_score, GridSearchCV, and other sklearn model-selection utilities.

Parameters:

X (array-like of shape (n_samples, n_features)) – Test data to score.
y (ignored)

Returns:

score – Negative mean squared reconstruction error.

Return type:

Examples

>>> from sklearn.model_selection import cross_val_score
>>> scores = cross_val_score(PCA(n_components=2), X_scaled, cv=5)
>>> print(f"Mean CV score: {scores.mean():.4f}")

classmethod minka_mle(X)[source]#

Minka (2000) automatic-dimensionality estimate for PCA.

Closed-form Bayesian model selection on the PPCA evidence (Minka, T. P. 2000. Automatic Choice of Dimensionality for PCA. NIPS 13, pp. 598-604). Operates only on the covariance eigenvalues of X and is therefore very cheap; in the simulations Minka reports it beats cross-validation. Use it alongside the ekf-CV recommendation from select_n_components() as a fast cross-check.

Internally X is mean-centred before estimation (a PPCA assumption); it is not unit-variance scaled, because dividing each column by its standard deviation compresses the noise eigenvalues to near-zero values the MLE misreads as additional latent signal. If your columns are on wildly different scales, pass the analysis-scale X produced by your own preprocessing (e.g. SNV for spectral data) and accept the centring this method applies.

Parameters:: X (array-like of shape (n_samples, n_features)) – Data matrix.
Returns:: n_components – The MLE estimate of the effective dimensionality.
Return type:: int

References

Minka, T. P. (2000). Automatic Choice of Dimensionality for PCA. Advances in Neural Information Processing Systems, 13, 598-604.

See also

parallel_analysis: Horn (1965) eigenvalue-vs-null retention.
select_n_components: ekf cross-validation; pass return_consensus=True to report all three side by side.

classmethod parallel_analysis(X, *, n_simulations=200, quantile=0.95, surrogate='normal', scale=True, random_state=None)[source]#

Horn (1965) parallel analysis component-count estimate.

Generates n_simulations random matrices of the same shape as X, computes their eigenvalues, and retains every observed component whose eigenvalue exceeds the quantile of the null distribution at the same rank. Widely regarded in psychometrics as the best simple retention rule for PCA.

Parameters:

X (array-like of shape (n_samples, n_features)) – Data matrix.
n_simulations (int, default 200) – Number of random matrices drawn to build the null eigenvalue distribution.
surrogate ({"normal", "permutation"}, default "normal") –
How the null matrices are built.
- "normal": independent standard-normal entries, which is Horn’s original proposal. Fast, and exactly right when the columns really are Gaussian.
- "permutation": each column of the real data is permuted independently (Buja and Eyuboglu, 1992). This breaks the correlation between columns, which is what parallel analysis is testing for, while leaving every column’s own distribution untouched. Prefer it on process data, where a tag may be skewed, heavy-tailed, bounded at zero or quantised by its instrument: a Gaussian null then answers a question about Gaussian data rather than about this block.
quantile (float, default 0.95) – Quantile of the null eigenvalues used as the retention threshold. Horn’s original proposal was the mean (0.5); the more conservative 95th-percentile threshold is the modern recommendation.
scale (bool, default True) – Mean-centre and unit-variance scale X before estimation. Unlike minka_mle(), which is mean-centred but never unit-variance scaled, parallel analysis defaults to autoscaling here so wildly different column scales do not dominate the null comparison.
random_state (int, optional) – Seed for the null-matrix simulations.

Returns:

result – With keys:

n_components - number of components retained (can be 0 on pure noise).
observed_eigenvalues - eigenvalues of X after centring/scaling (np.ndarray of length min(n, p)).
null_threshold - per-rank quantile of the null eigenvalue distribution (same length as observed_eigenvalues). Plot it against observed_eigenvalues for the scree-versus-null picture the method is usually read from.
surrogate - which null was used, echoed back so a stored result says how it was produced.

Return type:

References

Horn, J. L. (1965). A rationale and test for the number of factors in factor analysis. Psychometrika, 30(2), 179-185.

Buja, A., & Eyuboglu, N. (1992). Remarks on parallel analysis. Multivariate Behavioral Research, 27(4), 509-540.

See also

minka_mle: closed-form PPCA evidence rule.
select_n_components: ekf cross-validation; pass return_consensus=True to report all three side by side.

classmethod select_n_components(X, *, max_components=None, cv=5, cv_scheme='ekf', n_repeats=1, selection_rule='min', min_q2_increase=0.01, scale_inside_folds=True, n_iter=50, tol=1e-06, random_state=None, return_consensus=False, threshold=None, **pca_kwargs)[source]#

Select the number of PCA components via cross-validation.

Evaluates every component count 1, 2, ..., max_components and recommends one via the configured selection_rule. The default cv_scheme="ekf" is the element-wise k-fold algorithm of Bro, Kjeldahl, Smilde & Kiers (2008, Anal. Bioanal. Chem. 390:1241-1251), which holds out individual cells of X and predicts them via EM-style imputation from a model that never sees their true values. This restores the prediction-independence requirement the legacy row-wise scheme violates, fixing the trivial-fit pathology where PRESS shrinks monotonically with components.

Parameters:

X (array-like of shape (n_samples, n_features)) – Training data. With the default scale_inside_folds=True pass the raw, unscaled X; mean-centring and unit-variance scaling are fit on each fold’s in-fold cells. Pre-scale it yourself only with scale_inside_folds=False.
max_components (int, optional) – Maximum number of components to evaluate. Default is min(n_samples - 1, n_features).
cv (int or sklearn CV splitter, default 5) – For cv_scheme="ekf": the integer number of element-folds (splitter objects are ignored). For cv_scheme="ek": the number of row groups, and of column groups. For cv_scheme="row_wise": either an integer K (fed to KFold) or any sklearn splitter. Ignored by "sacv" and "gcv", which hold nothing out.
cv_scheme ({"ekf", "ckf", "ek", "sacv", "gcv", "row_wise"}, default "ekf") –
How a held-out value is produced. Every one of these except "row_wise" keeps the prediction independent of the value being predicted; they differ in what they hold out and what they cost.
- "ekf", the default: element-wise k-fold with EM imputation. Scattered cells are held out and each is imputed from a model that never saw it. Fits its centring and scaling inside every fold, and is the only scheme here that takes a block with missing cells. Cost: n_folds * n_repeats * max_components decompositions.
- "ckf": column-wise k-fold. Groups of columns are held out, and their values are predicted from scores computed on the retained columns only, so no held-out value appears in the score that predicts it. Information still reaches the model through the loadings, which are fitted on the whole block, and that is the respect in which "ekf" is stricter. Offered because a great deal of published chemometrics uses it, so a number that has to line up with a paper may need it, and because on a block with far more columns than samples it costs one decomposition rather than n_folds * n_repeats * max_components.
- "ek": the two-model scheme of Eastment and Krzanowski (1982). An element is predicted by a score from a model without its column and a loading from a model without its row. This is what Simca-P reports and what pcaMethods::Q2 computes by default, so use it when a number has to line up with either. Cost: 2 * n_folds decompositions.
- "sacv": leave-one-cell-out, approximated by inflating each residual by the leverage of the cell that produced it, after Josse and Husson (2012). The cheap version of holding cells out one at a time. Cost: one decomposition.
- "gcv": the same idea with a single averaged leverage instead of one per cell. Blunter, and it tends to keep more components. Cost: one decomposition. This and "sacv" are the defaults in FactoMineR and missMDA.
- "row_wise": deprecated, removed in 2.0. See the warning admonition below.
"ek", "sacv" and "gcv" factorise the matrix directly, so they raise on a block with missing cells and they ignore scale_inside_folds, n_repeats, n_iter and tol. They also report no per-fold spread, so selection_rule="1se" has nothing to work with; "min" takes the global optimum, where FactoMineR instead stops at the first local worsening, which can return a smaller count on the same data.
n_repeats (int, default 1) – Repeat the ekf pass with a fresh random fold permutation this many times. Each repeat covers every cell exactly once; n_repeats > 1 narrows the per-component PRESS standard error (helpful when the 1-SE rule sits on a borderline) at roughly linear extra runtime. Ignored under cv_scheme="row_wise".
selection_rule ({"min", "1se", "q2_increment"}, default "min") – How the recommended component count is chosen. "min" is the GlobalMin criterion Bro 2008 pairs with ekf - the component count with the lowest pooled PRESS. "1se" is the one- standard-error rule (needs per_fold_press, available under both schemes). "q2_increment" is the Wold’s-R-style cumulative-$Q^2$ threshold from PR #371; min_q2_increase sets the threshold.
min_q2_increase (float, default 0.01) – Threshold used only when selection_rule="q2_increment".
scale_inside_folds (bool, default True) –
With the default, mean-centring and unit-variance scaling constants are fit on each fold’s in-fold cells and applied to the whole matrix before EM, removing the centring/scaling leakage of the prior implementation. Set to False to reproduce the previous behaviour (column mean recomputed each EM iteration from the imputed matrix, no scaling); this is useful only when X is already pre-scaled, and a SpecificationWarning is emitted because scaling constants fit on the full matrix leak into every element-fold. Ignored under cv_scheme="row_wise".

Pass the raw, unscaled X under the default. In-fold re-standardisation overwrites whatever scaling the caller applied, so two deliberately different strategies (autoscale versus Pareto, say) become the same model and report the same PRESS: a comparison between them shows no difference for reasons that have nothing to do with the data. A SpecificationWarning is emitted when X arrives already centred and unit-variance scaled, which is the detectable half of that case; a block scaled some other way cannot be recognised, so the rule is the caller’s to keep. Same contract as PLS.select_n_components().
n_iter (int and float, default 50 and 1e-6) – EM iteration cap and convergence tolerance for the ekf imputation step. Ignored under cv_scheme="row_wise".
tol (int and float, default 50 and 1e-6) – EM iteration cap and convergence tolerance for the ekf imputation step. Ignored under cv_scheme="row_wise".
random_state (int, optional) – Seed for the ekf element-fold permutation.
return_consensus (bool, default False) – When True, also cross-check the CV recommendation against two cheap alternative selectors: Minka’s PPCA MLE (minka_mle()) and Horn’s parallel analysis (parallel_analysis()). The result Bunch then gains the minka_n_components, parallel_analysis_n_components, consensus, and consensus_counts keys (see Returns).
threshold (float, optional) – Deprecated. The original Wold PRESS-ratio cutoff. Passing it emits a DeprecationWarning; the value is ignored. Use selection_rule="q2_increment" (and tune min_q2_increase) for a comparable parsimony preference.
**pca_kwargs – Additional keyword arguments passed to the PCA() constructor under cv_scheme="row_wise" (e.g. algorithm="nipals"). Ignored under cv_scheme="ekf" because ekf runs its own SVD loop.

Returns:

result – With keys:

n_components - recommended number of components (int).
press - pooled PRESS per component count (pd.Series, indexed 1..A_max). Under cv_scheme="ekf" with scale_inside_folds=True this is measured in the space each fold was fitted in, so every variable weighs the same; see press_input_units for the other scale.
press_input_units - the same curve in the units of the matrix that was passed in, for comparing prediction error against instrument error (pd.Series, indexed 1..A_max).
per_fold_press - per-fold PRESS contributions (pd.DataFrame, A_max rows x n_folds * n_repeats columns under ekf; a single fold_1 column under row-wise).
se_press - standard error of the per-fold PRESS curve (pd.Series, indexed 1..A_max). Drives the 1-SE rule.
q2_se - the same standard error rescaled onto the Q2 scale (se_press / null_model_ss, pd.Series indexed 1..A_max), i.e. the half-width of a +/-1 SE band around q2.
press_ratio - PRESS_a / PRESS_{a-1} for inspection (pd.Series, indexed 2..A_max).
q2 - cross-validated $R^2_X$ per component count (pd.Series, indexed 1..A_max). Computed as 1 - press / null_model_ss, where the null model predicts each held-out cell by its in-fold column mean, measured the same way press is. Directly comparable to r2_cumulative_ and to PLS’s r2y_validated.
q2_per_variable - that same quantity split by variable (pd.DataFrame, A_max rows x K columns), which is what shows whether one column is carrying the pooled figure. All NaN under cv_scheme="row_wise", which has no per-cell error to split.
cv_scores - alias of per_fold_press under ekf, or per-fold negative MSE from cross_val_score under row-wise (preserved for back-compat).
cv_scheme - the scheme used ("ekf" or "row_wise").
selection_rule - the rule used to pick n_components.

When return_consensus=True, the Bunch additionally carries:

minka_n_components - the Minka PPCA MLE estimate (int).
parallel_analysis_n_components - Horn’s parallel-analysis estimate (int).
consensus - "agree" if the three integer estimates (CV recommendation, Minka, parallel analysis) span at most 1, otherwise "disagree".
consensus_counts - the tuple (recommended, minka_n, parallel_analysis_n).

Return type:

References

Bro, R., Kjeldahl, K., Smilde, A. K., & Kiers, H. A. L. (2008). Cross-validation of component models: a critical look at current methods. Anal. Bioanal. Chem., 390(5), 1241-1251.

Camacho, J., & Ferrer, A. (2012). Cross-validation in PCA models with the element-wise k-fold (ekf) algorithm: theoretical aspects. J. Chemometrics, 26(7), 361-373.

Eastment, H. T., & Krzanowski, W. J. (1982). Cross-validatory choice of the number of components from a principal component analysis. Technometrics, 24(1), 73-77.

Josse, J., & Husson, F. (2012). Selecting the number of components in principal component analysis using cross-validation approximations. Computational Statistics & Data Analysis, 56(6), 1869-1879.

Warning

cv_scheme="row_wise" is deprecated since 1.84 and will be removed in 2.0. It emits a DeprecationWarning and a SpecificationWarning. It suffers from the trivial-fit problem: holding out whole rows and projecting them back via transform() lets the held-out row’s own values reach its prediction, so PRESS shrinks monotonically with the component count and reaches zero once the components equal the variables. It measures compression, not prediction, and cannot select a component count. Use "ekf", "ek", "sacv" or "gcv".

detect_outliers(conf_level=0.95)[source]#

Detect outlier observations using SPE and Hotelling’s T² diagnostics.

Combines two approaches:

Statistical limits - observations exceeding the SPE or T² limit at conf_level are flagged.
Robust ESD test - the generalized ESD test identifies observations that are unusual relative to the rest of the data, even if they fall below the statistical limit. The mean/std variant is used here; the underlying detect_outliers_esd also offers an opt-in robust median/MAD variant.

An observation can be flagged for one or both reasons.

Parameters:

conf_level (float, default 0.95) – Confidence level in [0.8, 0.999]. Controls both the statistical limits and the ESD test’s significance level (alpha = 1 - conf_level).

Returns:

outliers – Sorted from most severe to least. Each dict contains:

observation - index label of the observation
outlier_types - list of "spe" and/or "hotellings_t2"
spe - SPE value for this observation
hotellings_t2 - T² value for this observation
spe_limit - SPE limit at the given confidence level
hotellings_t2_limit - T² limit at the given confidence level
severity - max(spe/spe_limit, t2/t2_limit), rounded to 4 decimals. A ratio whose denominator is 0 (perfect-fit SPE limit) or non-finite (T2 limit when A == N) is treated as 0 and does not contribute to the ranking.

Return type:

list of dict

Examples

>>> pca = PCA(n_components=3).fit(X_scaled)
>>> outliers = pca.detect_outliers(conf_level=0.95)
>>> for o in outliers:
...     print(f"{o['observation']}: {o['outlier_types']} (severity={o['severity']})")

PLS#

class process_improve.multivariate.methods.PLS(n_components, *, scale=True, max_iter=1000, tol=np.float64(1.4901161193847656e-08), copy=True, warn_on_uncentred=True, missing_data_settings=None)[source]#

Bases: _LatentVariableModel, RegressorMixin, TransformerMixin, BaseEstimator

Projection to Latent Structures (PLS) regression with diagnostics.

Implements PLS via the NIPALS algorithm with production diagnostics: SPE, Hotelling’s T², score contributions, and outlier detection. The API mirrors PCA so that model.scores_, model.spe_, and model.detect_outliers() work identically for both model types.

Parameters:

n_components (int) – Number of latent components to extract. None asks for as many components as the data supports: it is resolved at fit time to min(n_samples, n_features), which is also the ceiling an explicit request is clamped to (with a SpecificationWarning). The resolved count is available on the fitted attribute n_components_.
scale (bool, default=True) –
Mean-center and unit-variance-scale both the X and Y blocks internally before fitting (ddof=1, done with MCUVScaler). This mirrors sklearn.cross_decomposition.PLSRegression, whose scale=True default also scales X and Y; the parameter exists so PLS is a drop-in for the sklearn estimator. Predictions, predictions_ and beta_coefficients_ are returned on the original (un-scaled) data scale. When you scale externally (e.g. with MCUVScaler), set scale=False to avoid the (idempotent) double scaling. Note: the cross-validation helpers (select_n_components()) always re-fit an MCUVScaler inside each training fold regardless of this flag.

scale=False fits no intercept, so both blocks must already be centred. A response left on its natural scale is the trap: predictions come out offset by the response mean, and R² / Q² go large and negative on data that does contain a relationship. fit emits an UncentredDataWarning when either block’s column means are large relative to their spread; it does not centre for you, because scale=False means “touch nothing”. Set warn_on_uncentred=False when that fit is deliberate.
max_iter (int, default=1000) – Maximum number of iterations for the NIPALS algorithm.
tol (float, default=sqrt(machine epsilon)) – Relative convergence tolerance for the NIPALS algorithm: the change between two successive score-vector iterations, relative to the norm of the current score vector (see terminate_check()).
copy (bool, default=True) – Whether to copy X and Y before fitting.
warn_on_uncentred (bool, default=True) –
Emit the UncentredDataWarning described under scale when scale=False and a block arrives un-centred. Set it to False for a fit that is un-centred on purpose (a demonstration of the offset, or a test that some other centring check fires), where the diagnostic is the expected outcome rather than a problem.

This is the narrowest of the three opt-outs, and the one to reach for first. It silences the check for this model only, so an unrelated SpecificationWarning raised elsewhere in the same block, or by this same fit call, still arrives. Filtering UncentredDataWarning as a category is next narrowest; suppressing all of SpecificationWarning is the blunt instrument, and hides clamped component counts and NIPALS non-convergence along with it.

Has no effect when scale=True: the model centres both blocks itself, so the condition cannot arise. Like every constructor parameter it is stored verbatim and survives clone(), so a deliberately un-centred fit stays quiet inside a Pipeline or a grid search.
missing_data_settings (dict or None, default=None) –
Settings for the NIPALS fit when the data has missing cells. Keys:
- md_method: "nipals" (the default, and the only one implemented), or "tsr" / "pmp", which are recognised and raise NotImplementedError. Any other value is refused. This is a different and smaller set than the method= accepted by project() and the contribution helpers, which do implement "tsr", "scp" and "pmp".
- md_tol and md_max_iter: the NIPALS convergence tolerance and iteration cap. They default to this model’s tol and max_iter, so set those instead unless you need the fit and the missing-data path to differ.
fitting) (Attributes (after)
--------------------------
n_components – The resolved number of components actually fitted (the constructor parameter clamped to min(n_samples, n_features); the parameter itself is left as the user set it, including None).
scores (pd.DataFrame of shape (n_samples, n_components)) – X-block score matrix (T). This is the primary score matrix; equivalent to x_scores in older versions.
y_scores (pd.DataFrame of shape (n_samples, n_components)) – Y-block score matrix (U).
x_loadings (pd.DataFrame of shape (n_features, n_components)) – X-block loading matrix (P).
y_loadings (pd.DataFrame of shape (n_targets, n_components)) – Y-block loading matrix (C).
x_weights (pd.DataFrame of shape (n_features, n_components)) – X-block weight matrix (W).
y_weights (pd.DataFrame of shape (n_targets, n_components)) – Y-block weight matrix.
direct_weights (pd.DataFrame of shape (n_features, n_components)) – Direct (W*) weights: W (P'W)^{-1}. Used for direct projection T = X @ W*.
beta_coefficients (pd.DataFrame of shape (n_features, n_targets)) – Regression coefficients linking X directly to Y.
predictions (pd.DataFrame of shape (n_samples, n_targets)) – Y predictions from the training data.
spe (pd.DataFrame of shape (n_samples, n_components)) – Per-row SPE diagnostic; stored as the square root of the row sum-of-squared X-residuals (so it is on the residual scale, not the squared scale). One column per component, not one value per row: column a is the SPE of the model truncated at a components, and the last column is the value at the full fitted model. Reach for a single number per observation with model.spe_.iloc[:, -1], not with np.asarray(model.spe_).ravel(): ravel happens to give the right answer at one component and silently gives n_samples * n_components values above it.
hotellings_t2 (pd.DataFrame of shape (n_samples, n_components)) – Cumulative Hotelling’s T² statistic. Per-component, exactly as spe_ above: column a uses the first a components and the last column is the value at the full fitted model.
r2_per_component (pd.Series of length n_components) – Fractional R² (on Y) explained by each component.
r2_cumulative (pd.Series of length n_components) – Cumulative R² (on Y) after each component.
r2_per_variable (pd.DataFrame of shape (n_features, n_components)) – Per-variable cumulative R² for X after each component.
r2y_per_variable (pd.DataFrame of shape (n_targets, n_components)) – Per-variable R² for Y after each component.
rmse (pd.DataFrame of shape (n_targets, n_components)) – Root mean squared error of Y predictions per component, on the original (un-scaled) Y scale, consistent with predictions_ and prediction_interval.
explained_variance (np.ndarray of shape (n_components,)) – Variance explained by each component in X.
scaling_factor_for_scores (pd.Series of length n_components) – Standard deviation per score (sqrt of explained variance).
has_missing_data (bool) – Whether the training data contained missing values.
fitting_info (dict) – Timing and iteration info from the fitting algorithm.

See also

PCA: Principal Component Analysis.
MCUVScaler: Mean-center unit-variance scaler.

References

Abdi, “Partial least squares regression and projection on latent structure regression (PLS Regression)”, 2010, DOI: 10.1002/wics.51

Examples

>>> import pandas as pd
>>> from process_improve.multivariate.methods import PLS, MCUVScaler
>>> X = pd.DataFrame({"A": [1, 2, 3, 4], "B": [4, 3, 2, 1]})
>>> Y = pd.DataFrame({"y": [2.1, 3.9, 6.2, 7.8]})
>>> pls = PLS(n_components=1)
>>> pls = pls.fit(MCUVScaler().fit_transform(X), MCUVScaler().fit_transform(Y))
>>> pls.scores_.shape
(4, 1)

get_feature_names_out(input_features=None)[source]#

Return the output column names of transform().

PLS’s transform returns the X scores (T matrix), labelled ["T1", "T2", ..., "T{n_components}"]. The input_features argument is accepted (Pipeline introspection passes it through) but unused: the output column count is the fitted n_components, not the input feature count.

Used by set_output() (sklearn 1.2+) to label the DataFrame view of the scores when set_output(transform="pandas") is on, and by Pipeline introspection.

Return type:: ndarray

predictions_vs_observed_plot(*, y_observed, variable=None, settings=None, fig=None)#

Generate an observed-vs-predicted (parity) plot for a fitted PLS model.

Plots the calibration predictions against the observed Y values, with a y = x reference line and an RMSE annotation. Points lying close to the reference line indicate good predictions.

Parameters:

model (PLS object) – A fitted PLS model generated by this library.
y_observed (array-like of shape (n_samples, n_targets)) – The observed Y values, on the same scale as the data used to fit the model (for example the scaled Y from MCUVScaler).
variable (str, optional) – Which Y-variable to plot. Defaults to the first Y-variable.

settings (dict) –

Default settings:

{
    "title": "Observed vs predicted ...",  # str: overall plot title
    "marker_color": None,           # str|None: data-marker colour; None uses the theme
    "reference_color": "#9CA3AF",   # str: colour of the y = x line
    "html_image_height": 500,       # int: image height in pixels
    "html_aspect_ratio_w_over_h": 1.0,   # float: width as ratio of height
    "template": "pi_journal",         # str: registered Plotly theme name
}

fig (go.Figure, optional) – An existing figure to draw onto. A new figure is created if omitted.

Return type:

Figure

Examples

>>> pls.predictions_vs_observed_plot(y_observed=Y_scaled)
>>> pls.predictions_vs_observed_plot(y_observed=Y_scaled, variable="quality")

coefficient_plot(variable=None, settings=None, fig=None)#

Generate a bar plot of the PLS regression coefficients.

Shows beta_coefficients_ for one Y-variable: one bar per X-variable, mapping the (preprocessed) X onto the predicted Y. Tall bars mark the X-variables that most strongly drive the prediction.

Parameters:

model (PLS object) – A fitted PLS model generated by this library.
variable (str, optional) – Which Y-variable’s coefficients to plot. Defaults to the first one.

settings (dict) –

Default settings:

{
    "title": "Regression coefficients ...",  # str: overall plot title
    "bar_color": None,              # str|None: bar colour; None uses the theme
    "html_image_height": 500,       # int: image height in pixels
    "html_aspect_ratio_w_over_h": 16/9,  # float: width as ratio of height
    "template": "pi_journal",         # str: registered Plotly theme name
}

fig (go.Figure, optional) – An existing figure to draw onto. A new figure is created if omitted.

Return type:

Figure

Examples

>>> pls.coefficient_plot()
>>> pls.coefficient_plot(variable="quality")

target_projection(X, response=None)#

Target-projected (TP) component of a fitted PLS model for one response.

Target projection (Kvalheim and Karstang, 1989) rotates the PLS solution so that a single latent component carries all of the predictive information for one response. The component points along the regression vector $b$ (the column of beta_coefficients_ for that response):

\[w_{\text{TP}} = \frac{b}{\lVert b \rVert}, \qquad t_{\text{TP}} = X\, w_{\text{TP}}, \qquad p_{\text{TP}} = \frac{X^\top t_{\text{TP}}}{t_{\text{TP}}^\top t_{\text{TP}}}.\]

The TP component is the basis for the selectivity ratio (selectivity_ratio()).

Parameters:

model (PLS) – A fitted PLS model (must expose beta_coefficients_).
X (array-like of shape (n_samples, n_features)) – Preprocessed data, scaled the same way as the training data (for example with MCUVScaler).
response (str or int or None, default=None) – Which response (Y column) to project onto. None is allowed only for a single-response model; otherwise pass the response label (or its integer position).

Returns:

With fields scores (pd.Series, the TP scores per sample), loadings (pd.Series, the TP loading per feature), weights (pd.Series, the unit TP weight per feature) and response (the resolved response label).

Return type:

Raises:

ValueError – If the model is not a fitted PLS, the response selector is invalid, or the regression vector / TP scores are degenerate (~0).

References

Kvalheim, O. M. and Karstang, T. V. (1989). Interpretation of latent-variable regression models. Chemometrics and Intelligent Laboratory Systems, 7(1-2), 39-51.

Examples

>>> pls = PLS(n_components=3).fit(X_scaled, y_scaled)
>>> tp = pls.target_projection(X_scaled)        # bound convenience method
>>> tp.scores.head()

See also

selectivity_ratio: Per-variable explained/residual ratio on the TP component.

selectivity_ratio(X, response=None, *, conf_level=0.95)#

Compute the selectivity ratio of each feature on the target-projected component.

The selectivity ratio (Rajalahti et al., 2009) ranks each feature by how much of its variance the predictive (target-projected) direction explains. On the TP component (target_projection()), for feature $j$:

\[\text{SR}_j = \frac{\text{SS}_{\text{explained},j}} {\text{SS}_{\text{residual},j}} = \frac{p_{\text{TP},j}^2\, (t_{\text{TP}}^\top t_{\text{TP}})} {\sum_i (x_{ij} - t_{\text{TP},i}\, p_{\text{TP},j})^2}.\]

A large SR means the feature is well aligned with the predictive direction. Unlike VIP, it is a true explained/residual variance ratio and can be compared against an F distribution. Note that two collinear features carry near-identical SR: the selectivity ratio ranks predictive relevance, it does not break ties between mutually collinear features.

Parameters:

model (PLS) – A fitted PLS model (must expose beta_coefficients_).
X (array-like of shape (n_samples, n_features)) – Preprocessed data, scaled the same way as the training data (for example with MCUVScaler).
response (str or int or None, default=None) – Which response to compute SR for. None returns a feature-by-response DataFrame when the model has several responses, or a Series for a single-response model.
conf_level (float, default=0.95) – Confidence level for the advisory F-based critical value, attached to the result’s .attrs["f_critical"].

Returns:

Selectivity ratios indexed by feature. A Series for one response (with f_critical / conf_level / response in .attrs), or a feature-by-response DataFrame when response is None and the model has several responses.

Return type:

pd.Series or pd.DataFrame

Raises:

ValueError – If the model is not a fitted PLS or the response selector is invalid.

References

Rajalahti, T., Arneberg, R., Berven, F. S., Myhr, K.-M., Ulvik, R. J. and Kvalheim, O. M. (2009). Biomarker discovery in mass spectral profiles by means of selectivity ratio plot. Chemometrics and Intelligent Laboratory Systems, 95(1), 35-48.

Examples

>>> pls = PLS(n_components=3).fit(X_scaled, y_scaled)
>>> pls.selectivity_ratio(X_scaled).sort_values(ascending=False).head()

See also

target_projection: The target-projected component the ratio is built on.
vip: Variable Importance in Projection, an alternative importance measure.

y_scores_: ndarray | DataFrame#

y_weights_: ndarray | DataFrame#

y_loadings_: ndarray | DataFrame#

fitting_info_: dict[str, ndarray | int | float]#

scores_#

Expose a private ndarray as a lazily-built, cached pandas.DataFrame (ENG-18).

Declared as a class attribute, e.g.:

scores_   = _LazyFrame("_scores",   index="_sample_index",  columns="_component_names")
loadings_ = _LazyFrame("_loadings", index="_feature_names", columns="_component_names")

On an unfitted model the backing ndarray is absent, so getattr raises AttributeError - the same “not fitted” signal as before this change, so hasattr / check_is_fitted behave identically.

spe_#

Expose a private ndarray as a lazily-built, cached pandas.DataFrame (ENG-18).

Declared as a class attribute, e.g.:

scores_   = _LazyFrame("_scores",   index="_sample_index",  columns="_component_names")
loadings_ = _LazyFrame("_loadings", index="_feature_names", columns="_component_names")

On an unfitted model the backing ndarray is absent, so getattr raises AttributeError - the same “not fitted” signal as before this change, so hasattr / check_is_fitted behave identically.

x_loadings_#

Expose a private ndarray as a lazily-built, cached pandas.DataFrame (ENG-18).

Declared as a class attribute, e.g.:

scores_   = _LazyFrame("_scores",   index="_sample_index",  columns="_component_names")
loadings_ = _LazyFrame("_loadings", index="_feature_names", columns="_component_names")

On an unfitted model the backing ndarray is absent, so getattr raises AttributeError - the same “not fitted” signal as before this change, so hasattr / check_is_fitted behave identically.

x_weights_#

Expose a private ndarray as a lazily-built, cached pandas.DataFrame (ENG-18).

Declared as a class attribute, e.g.:

scores_   = _LazyFrame("_scores",   index="_sample_index",  columns="_component_names")
loadings_ = _LazyFrame("_loadings", index="_feature_names", columns="_component_names")

On an unfitted model the backing ndarray is absent, so getattr raises AttributeError - the same “not fitted” signal as before this change, so hasattr / check_is_fitted behave identically.

fit(X, Y, sample_weight=None)[source]#

Fit a projection to latent structures (PLS) model to the data.

Parameters:

X (array-like, shape (n_samples, n_features)) – Training data, where n_samples is the number of samples (rows) and n_features is the number of features (columns).
Y (array-like, shape (n_samples, n_targets)) – Training data, where n_samples is the number of samples (rows) and n_targets is the number of target outputs (columns).
sample_weight (array-like of shape (n_samples,), optional) – Non-negative row weights for a weighted PLS fit (#394). NIPALS is run on sqrt(w)-rescaled X and Y, which is equivalent to weighting the cross-products X' W u and Y' W t. Loadings, weights and beta are computed correctly; scores are returned on the original sample scale. Zero weights effectively exclude the corresponding rows (sample_weight=[1,1,0,0,1] reproduces the unweighted fit on rows [0,1,4]). Forwarded to score() / r2_score for any caller that also threads it through.

Returns:

Model object.

Return type:

PLS

References

Abdi, “Partial least squares regression and projection on latent structure regression (PLS Regression)”, 2010, DOI: 10.1002/wics.51

transform(X, Y=None)[source]#

Project X (and optionally Y) into the latent space.

Parameters:

X (array-like of shape (n_samples, n_features)) – Data to transform.
Y (array-like of shape (n_samples, n_targets), optional) – Ignored. Present for API compatibility with sklearn pipelines.

Returns:

X_scores – Projected X data (scores).

Return type:

pd.DataFrame of shape (n_samples, n_components)

fit_transform(X, Y=None)[source]#

Fit the model and return X scores.

Parameters:

X (array-like of shape (n_samples, n_features))
Y (array-like of shape (n_samples, n_targets)) – Required despite the None default: the default is present only for signature-compatibility with the sklearn fit_transform protocol, and PLS cannot be fitted without responses. Omitting Y raises ValueError, under python -O as well.

Returns:

X_scores

Return type:

pd.DataFrame of shape (n_samples, n_components)

score(X, Y, sample_weight=None)[source]#

Return the R² score for the prediction.

Parameters:

X (array-like of shape (n_samples, n_features))
Y (array-like of shape (n_samples, n_targets)) – True target values.
sample_weight (array-like of shape (n_samples,), optional)

Returns:

score – R² of self.predict(X) w.r.t. Y.

Return type:

predict(X)[source]#

Predict Y for new observations.

Returns just the predicted y_hat so the call satisfies the scikit-learn RegressorMixin contract (and therefore composes inside Pipeline, cross_val_score(), and GridSearchCV). For the rich diagnostic view (scores, Hotelling’s T², SPE, plus y_hat), see diagnose().

Parameters:: X (array-like of shape (n_samples, n_features))
Returns:: y_hat – Predicted target values, indexed by X’s rows and labelled with the target column names captured during fit.
Return type:: pd.DataFrame of shape (n_samples, n_targets)

See also

diagnose: richer per-prediction diagnostics.

Examples

>>> y_pred = pls.predict(X_new)
>>> diag = pls.diagnose(X_new)   # for scores / T² / SPE

diagnose(X)[source]#

Project new data and compute predictions plus diagnostics.

This is the rich view that predict() used to return before 1.35.0: alongside y_hat it reports the X scores, cumulative Hotelling’s T², and SPE for every row of X so the user can flag out-of-model observations and read their predicted Y from one call.

Parameters:: X (array-like of shape (n_samples, n_features))
Returns:: result – With keys scores, hotellings_t2, spe, y_hat.
Return type:: sklearn.utils.Bunch

See also

predict: sklearn-compatible call returning just y_hat.

Examples

>>> result = pls.diagnose(scaler_x.transform(X_new))
>>> result.y_hat           # Predicted Y values
>>> result.spe             # SPE for each new observation
>>> result.hotellings_t2   # T² for each new observation

project(X, *, method='tsr', ridge=0.0)[source]#

Estimate scores, prediction and diagnostics for rows with missing values.

Whereas transform() and diagnose() propagate NaN into the scores, this method estimates the scores of partially-observed rows from the observed columns only, using the missing-data estimators of Arteaga and Ferrer (2002): trimmed score regression ("tsr", the default and statistically the strongest), single-component projection ("scp", the score step of NIPALS itself: project onto the observed part of each weight vector, deflate with the loadings), or projection to the model plane ("pmp"). Rows with no missing values take the standard complete-data path, so their scores are bitwise identical to transform(). As the observed part grows to the whole row, TSR and SCP tend to the model’s own scores; PMP, the least-squares fit of the observed columns onto the loadings, does not for a PLS model, whose scores come from the weights rather than the loadings, so prefer the other two here.

This is the “batch so far” primitive for predicting the final quality of a running batch: the future part of the unfolded row is missing by construction, and the prediction at each decision point is this projection followed by the ordinary Y regression (Garcia-Munoz, Kourti and MacGregor, 2004; Flores-Cerrillo and MacGregor, 2004).

Parameters:

X (array-like of shape (n_samples, n_features)) – New observations on the original (unscaled) X units when the model was fitted with scale=True, or in the model’s scaled space otherwise, exactly as transform() expects. NaN marks a missing entry; rows that are entirely NaN are rejected.
method ({"tsr", "scp", "pmp"}, default="tsr") – The score estimator; see process_improve.multivariate._projection.
ridge (float, default=0.0) – Non-negative regularisation added to the matrix inverted by the "tsr" and "pmp" estimators. Raise it above zero when condition_number reports near-singularity (typically very early in a batch, when few columns are observed).

Returns:

result – With keys scores (DataFrame), y_hat (DataFrame, on the original Y units), hotellings_t2 (Series; total over all components), spe (Series; square root of the residual sum of squares over the observed columns only), condition_number (Series; 1.0 when nothing is missing) and n_observed (Series). SPE and T2 of a partially-observed row must be compared against limits built from the same missingness pattern, not the full-observation limits.

Return type:

projection_matrix(observed, *, method='tsr', ridge=0.0)[source]#

Build the fixed linear operator mapping observed columns to score estimates.

For a fixed missingness pattern, every estimator in project() is a fixed linear map t_hat = M @ z_observed on the model’s scaled X space. This method exposes that matrix so callers that reuse one pattern many times (an online monitor at time sample k, or a mid-course optimiser treating the candidate future columns as observed) can precompute it once. Note the matrix acts on scaled values: when the model was fitted with scale=True, apply the internal centring and scaling first (as project() does).

Parameters:

observed (array-like) – Either a boolean mask of length n_features_in_ (True = observed), or a list of feature labels to treat as observed.
method ({"tsr", "scp", "pmp"}, default="tsr")
ridge (float, default=0.0)

Returns:

result – With keys matrix (DataFrame, n_components x n_observed, columns labelled by the observed features), condition_number (float) and method.

Return type:

invert(y_desired, *, null_space_coordinates=None)[source]#

Invert the PLS model: find inputs that yield a desired response.

PLS is normally used in the forward direction (predict()): given inputs X, predict the response Y. Model inversion runs the model backwards: fix the response you want (y_desired) and solve for an input vector that the model predicts will achieve it. This is the basis of latent-variable product and process design (Jaeckle and MacGregor, 2000).

Because a PLS model usually retains more components A than the rank r of the response, the target pins down only r of the A score directions and the inversion is underdetermined: a whole (A - r)-dimensional family of input vectors yields the same prediction. That family is the null space. This method returns the minimum-norm (direct-inversion) solution together with an orthonormal basis for the null space, so callers can move along it to satisfy secondary criteria (cost, safety, operability) without changing the predicted response.

For a single response (r = 1), García-Carrión et al. (2025) proved that this null space is the same linear space as the orthogonal space isolated by an O-PLS model with the same total number of components.

Parameters:

y_desired (float, array-like, pandas Series/DataFrame, or dict) – The desired response, on the original (un-scaled) Y scale. A scalar is accepted for a single-target model; otherwise supply one value per target. A Series/DataFrame/dict is aligned to the fitted target names; a plain array must follow the fitted target order.
null_space_coordinates (np.ndarray, optional) – Coordinates along the null-space basis, of length A - r (the null-space dimension). When given, the returned solution is tau_direct_inversion + null_space_basis @ null_space_coordinates reconstructed into the input space. All such solutions yield the same predicted response. When omitted, the minimum-norm (direct-inversion) solution is returned.

Returns:

result – With keys:

x_newpd.Series of shape (n_features,): The estimated input vector, on the original (un-scaled) X scale.
scorespd.Series of length A: The score vector (tau) of the solution.
y_hatpd.Series of length n_targets: The model’s prediction at x_new; equals y_desired up to numerical error, a check that the inversion is consistent.
null_space_basispd.DataFrame of shape (A, A - r): Orthonormal basis of the null space, in score coordinates. Empty (zero columns) when A == r and the solution is unique.
null_space_dimensionint: A - r, the number of free directions.
hotellings_t2float: Hotelling’s T² of the solution, to flag extrapolation beyond the calibration data. Compare against hotellings_t2_limit().

Return type:

See also

predict: the forward direction, X -> Y.
hotellings_t2_limit: confidence limit to judge hotellings_t2.

References

C. M. Jaeckle and J. F. MacGregor, “Industrial applications of product design through the inversion of latent variable models”, Chemometrics and Intelligent Laboratory Systems, 50 (2000): 199-210, DOI: 10.1016/S0169-7439(99)00058-1.

S. García-Carrión et al., “On the equivalence between null space and orthogonal space in latent variable regression modeling”, Journal of Chemometrics, 39 (2025): e70057, DOI: 10.1002/cem.70057.

Examples

>>> result = pls.invert(y_desired=25.0)
>>> result.x_new              # input vector giving the target response
>>> result.null_space_basis   # directions that leave the response fixed
>>> pls.predict(result.x_new.to_frame().T)   # ~= 25.0

classmethod select_n_components(X, Y, *, max_components=None, cv=5, n_repeats=None, random_state=None, selection_rule='1se', scale_inside_folds=True, min_q2_increase=0.01, n_permutations=999, alpha=0.01, stability_threshold=0.6, **pls_kwargs)[source]#

Select the number of PLS components via cross-validation.

Fits PLS models on cross-validation training folds and evaluates the out-of-fold prediction error for every component count 1, 2, ..., max_components. Reports per-fold and pooled RMSECV plus the validated cumulative R² curves, and recommends a component count from one of three rules (see selection_rule below).

The defaults are the research-backed combination: the one-standard-error rule on top of repeated, shuffled K-fold CV, with MCUVScaler re-fit inside every training fold so test data never leaks into the centring/scaling estimates.

Unlike the calibration statistics stored on a fitted model (rmse_, r2_cumulative_), the metrics returned here estimate performance on unseen data and are therefore suitable for choosing n_components.

Parameters:

X (array-like of shape (n_samples, n_features)) – Training X. With the default scale_inside_folds=True the raw, unscaled X may be passed; scaling is fit inside every training fold.
Y (array-like of shape (n_samples, n_targets)) – Training Y. Same treatment as X under scale_inside_folds.
max_components (int, optional) – Maximum number of components to evaluate. Default is the largest value supported by every cross-validation training fold, min(min_fold_size, n_features).
cv (int or sklearn CV splitter, default 5) – If an integer, used as the n_splits of a shuffled KFold (or RepeatedKFold when n_repeats > 1). Any sklearn splitter object (e.g. KFold(10, shuffle=True) or LeaveOneOut()) is also accepted and is used as-is (n_repeats is then ignored).
n_repeats (int, optional) – Number of times the K-fold split is repeated with a fresh shuffle, used only when cv is an integer. The signature default is None, which is resolved to 10 inside the function (giving a cv * 10 per-fold sample for the 1-SE rule); pass 1 to disable repeats. Repeated K-fold’s standard errors are slightly optimistic because test folds overlap across repeats; that is fine for the 1-SE selection rule but should not be reported as an unbiased generalisation variance.
random_state (int, optional) – Seed forwarded to KFold / RepeatedKFold for reproducible shuffling. Ignored when cv is a pre-built splitter.
selection_rule ({"1se", "min", "q2_increment", "randomization"}, default "1se") – How the recommended component count is chosen. See SelectionRule for the rule semantics. "1se" is the default; "min" is the argmin RMSECV (the pre-1.28 default, prone to running to the maximum component count); "q2_increment" is the Wold’s-R-style cumulative-Q² threshold; "randomization" is Van der Voet’s (1994) permutation test (uses n_permutations and alpha) that picks the smallest model whose predictive ability is statistically indistinguishable from the reference (argmin RMSECV) one.
scale_inside_folds (bool, default True) –
When True (the default), fit a fresh MCUVScaler on each training fold’s X and Y, apply it to the held-out rows, fit PLS in scaled space, then inverse-transform the predictions so RMSECV is reported on the original Y scale. This removes the centring / scaling leakage of the prior default. Set to False to keep the pre-1.28 behaviour, in which case X and Y should already be scaled; a SpecificationWarning is emitted.

Pass the raw, unscaled blocks under the default. In-fold re-standardisation overwrites whatever scaling the caller applied, so two deliberately different strategies (autoscale versus Pareto, say) become the same model and report RMSECV identical to several decimal places: a comparison between them shows no difference for reasons that have nothing to do with the data. A SpecificationWarning is emitted when X arrives already centred and unit-variance scaled, which is the detectable half of that case; a block scaled some other way cannot be recognised, so the rule is the caller’s to keep.
min_q2_increase (float, default 0.01) – Threshold used only when selection_rule="q2_increment": the smallest increase in cumulative validated $Q^2_Y$ that justifies keeping an extra component.
n_permutations (int, default 999) – Used only when selection_rule="randomization": number of sign-flip permutations driving the Van der Voet test.
alpha (float, default 0.01) – Used only when selection_rule="randomization": significance level. The smallest component count whose Van der Voet p-value exceeds alpha is recommended. R’s pls::selectNcomp uses the same default; smaller values pick more parsimonious models.
stability_threshold (float, default 0.6) – For the per-repeat stability-selection diagnostic ("1se" / "min" rules with n_repeats > 1 only): the recommendation is judged selection_is_stable=True iff the modal vote share in selection_distribution is at least this fraction. Meinshausen & Bühlmann (2010, JRSS-B) suggest 0.6-0.9 for their variable-selection analogue; we default to the permissive end.
**pls_kwargs – Additional keyword arguments passed to the PLS() constructor (e.g. missing_data_settings).

Returns:

result – With keys:

n_components - recommended number of components (int).
rmsecv - pooled RMSECV per component count (pd.DataFrame, indexed 1..A; columns are the Y-variable names plus "total").
per_fold_rmsecv - per-fold total RMSECV (pd.DataFrame, indexed 1..A; one column per fold across all repeats). Drives the 1-SE rule.
se_rmsecv - standard error of the per-fold RMSECV per component count (pd.Series, indexed 1..A).
q2_se - standard error on the Q2 scale (the per-fold total PRESS standard error divided by the total Y sum-of-squares), i.e. the half-width of a +/-1 SE band around r2y_validated["total"] (pd.Series, indexed 1..A).
r2y_validated - validated cumulative $R^2_Y$ (pd.DataFrame, indexed 1..A; one column per Y-variable, then "total" and "scaled_total"). "total" pools the targets on the original Y scale, so a wide-ranging target dominates it; "scaled_total" weights every target equally, which is the pooling a fitted model’s r2_y_cumulative_ uses, so those two are the columns to compare fitted against held-out.
r2x_validated - validated cumulative $R^2_X$ (pd.DataFrame, indexed 1..A; columns are the X-variable names plus "total").
press - pooled Y prediction error sum of squares per component count (pd.Series, indexed 1..A).
cv_predictions - out-of-fold predictions of Y at the recommended component count, on the original Y scale (pd.DataFrame). For repeated K-fold, the first repeat’s held-out predictions are reported so each row appears exactly once.
selection_rule - the rule used to pick n_components.
randomization_pvalues - per-component Van der Voet right-tail p-values when selection_rule="randomization"; None otherwise.
selection_distribution - per-repeat vote share over candidate component counts (pd.Series indexed 1..A). Populated only for selection_rule in {"1se", "min"} and n_repeats > 1; None otherwise. A concentrated distribution signals a confident recommendation; a flat or multi-modal one flags it for review.
selection_mode - the most-voted component count, or None when selection_distribution is None.
selection_is_stable - True iff the modal vote share meets stability_threshold; None when no distribution was computed.

Return type:

Notes

The pooled RMSECV in rmsecv["total"] is the square root of the total PRESS over all fold-test rows divided by (N_eff * M) where N_eff = N * n_repeats under repeated CV; the per_fold_rmsecv column for fold f is the square root of fold-f’s sum-of-squared residuals over its own test rows.

References

Breiman, Friedman, Olshen & Stone (1984), CART, sec.3.4.3 (1-SE rule). Hastie, Tibshirani & Friedman, ESL, sec.7.10. Kohavi (1995, IJCAI) recommends 10-fold stratified CV for model selection.

Examples

>>> from sklearn.model_selection import KFold
>>> # Default: 1-SE on 10 x 5-fold repeated CV with in-fold scaling.
>>> result = PLS.select_n_components(X, Y, max_components=6, random_state=0)
>>> result.n_components, result.selection_rule
>>> # Opt-in to the older argmin-RMSECV rule:
>>> PLS.select_n_components(X, Y, max_components=6, selection_rule="min")
>>> # Caller-supplied splitter (n_repeats is ignored here):
>>> PLS.select_n_components(X, Y, cv=KFold(10, shuffle=True, random_state=0))

classmethod nested_cv(X, Y, *, max_components=None, outer_cv=5, inner_cv=5, n_inner_repeats=10, selection_rule='1se', scale_inside_folds=True, min_q2_increase=0.01, n_permutations=999, alpha=0.01, random_state=None, **pls_kwargs)[source]#

Nested cross-validation for an honest PLS performance estimate.

Outer loop splits the data into outer-train / outer-test; the inner loop runs select_n_components() on the outer-train (with the configured selection_rule over inner_cv * n_inner_repeats folds) to pick the component count; a final PLS is fit on the outer-train at that count and used to predict the outer-test. The accumulated out-of-fold predictions give RMSEP that is not optimism-biased by the selection decision - the headline number to report when a clean test set is not available.

Parameters:

X (array-like) – Training data. Treated as in select_n_components() (raw if scale_inside_folds=True, pre-scaled otherwise).
Y (array-like) – Training data. Treated as in select_n_components() (raw if scale_inside_folds=True, pre-scaled otherwise).
max_components (int, optional) – Forwarded to the inner select_n_components().
outer_cv (int or sklearn splitter, default 5) – Number of outer folds (or a custom splitter).
inner_cv (int, default 5) – Number of inner folds passed to select_n_components().
n_inner_repeats (int, default 10) – Number of inner-CV repeats per outer fold; the inner random_state is offset by the outer-fold index so each outer fold sees a fresh inner shuffle.
selection_rule (str, default "1se") – Selection rule applied inside the inner loop. See SelectionRule.
scale_inside_folds (bool, default True) – Mirrors select_n_components(). Also applied to the final outer-train fit, with the test-fold predictions inverse- transformed to the original Y scale before RMSEP accumulates.
min_q2_increase (float) – Forwarded to the inner select_n_components() per rule.
n_permutations (int) – Forwarded to the inner select_n_components() per rule.
alpha (float) – Forwarded to the inner select_n_components() per rule.
random_state (int, optional) – Seed for the outer-fold shuffle and the inner CV. The inner seed is offset per outer fold so each outer split sees a fresh shuffled inner CV.
**pls_kwargs – Forwarded to PLS for both the inner CV and the final outer-train fits.

Returns:

result – With keys:

rmsep - honest held-out RMSEP per Y column plus a "total" entry (pd.Series).
q2y - validated $Q^2_Y$ per Y column plus "total" (pd.Series).
cv_predictions - out-of-fold predictions of Y at the per-outer-fold selected component counts (pd.DataFrame on the original Y scale).
selected_components_per_fold - list of inner recommendations, one per outer fold.
selected_components_distribution - vote share over candidate counts (pd.Series).

Return type:

Notes

Runtime is roughly outer_cv * inner_cv * n_inner_repeats * max_components PLS fits. With the defaults that is 5 * 5 * 10 * max_components fits per call; for max_components=10 and a moderate dataset that completes in seconds. Drop n_inner_repeats if you need to bring it down further.

Examples

>>> from process_improve.multivariate import PLS
>>> result = PLS.nested_cv(X, Y, max_components=8, random_state=0)
>>> result.rmsep["total"]

detect_outliers(conf_level=0.95)[source]#

Detect outlier observations using SPE and Hotelling’s T² diagnostics.

Same approach as PCA.detect_outliers: combines statistical limits with the robust generalized ESD test.

Parameters:: conf_level (float, default 0.95) – Confidence level in [0.8, 0.999].
Returns:: outliers – Sorted from most severe to least. Each dict contains observation, outlier_types, spe, hotellings_t2, spe_limit, hotellings_t2_limit, severity.
Return type:: list of dict

Examples

>>> pls = PLS(n_components=3).fit(X_scaled, Y_scaled)
>>> outliers = pls.detect_outliers(conf_level=0.95)
>>> for o in outliers:
...     print(f"{o['observation']}: {o['outlier_types']}")

cross_validate(X, Y, *, cv='loo', n_bootstrap=0, conf_level=0.95, random_state=None, show_progress=True, sample_weight=None)[source]#

Cross-validate the PLS model and compute error bars for beta coefficients.

Refits the model on data subsets (jackknife, K-fold, or bootstrap), collects beta_coefficients_ from each refit, and computes confidence intervals. Also returns cross-validated predictions and prediction-error metrics (RMSE, Q²).

Parameters:

X (array-like of shape (n_samples, n_features)) – Predictor matrix (same data used for fit).
Y (array-like of shape (n_samples, n_targets)) – Response matrix (same data used for fit).
cv (int or "loo", default "loo") –
Cross-validation strategy:
- "loo" - leave-one-out (jackknife). Produces N resamples.
- int - number of folds for K-fold CV.
n_bootstrap (int, default 0) – If > 0, use bootstrap resampling instead of CV folds. The value specifies the number of bootstrap rounds. Overrides the cv parameter when set.
conf_level (float, default 0.95) – Confidence level for the beta-coefficient intervals, in (0, 1).
random_state (int or None, default None) – Random seed for reproducibility (K-fold shuffle and bootstrap).
show_progress (bool, default True) – Whether to display a tqdm progress bar.
sample_weight (np.ndarray of shape (n_samples,), optional) – Per-sample non-negative weights. Threaded into every sub-fit so each resample’s PLS uses the same weighting scheme as the parent model. Default None (all samples weighted equally).

Returns:

result – Dictionary-like object with the following keys:

Beta-coefficient uncertainty

beta_samplesnp.ndarray of shape (n_resamples, n_features, n_targets): Raw beta coefficients from every resample.
beta_meanpd.DataFrame of shape (n_features, n_targets): Mean beta across resamples.
beta_stdpd.DataFrame of shape (n_features, n_targets): Standard error of the beta coefficients.
beta_ci_lowerpd.DataFrame of shape (n_features, n_targets): Lower bound of the confidence interval.
beta_ci_upperpd.DataFrame of shape (n_features, n_targets): Upper bound of the confidence interval.
significantpd.DataFrame of shape (n_features, n_targets): True where the confidence interval excludes zero.

Prediction metrics

y_hat_cvpd.DataFrame of shape (n_samples, n_targets): Cross-validated predictions (out-of-fold). Only available for jackknife and K-fold; None for bootstrap.
pressfloat: Prediction Error Sum of Squares (sum over all Y elements). Only for jackknife / K-fold.
rmse_cvpd.Series of length n_targets: Root-mean-square error per Y variable (cross-validated). Only for jackknife / K-fold.
q_squaredpd.Series of length n_targets: Cross-validated R² (Q²) per Y variable. Only for jackknife / K-fold.

Metadata

n_resamplesint: Number of resamples performed.
methodstr: "jackknife", "kfold", or "bootstrap".
conf_levelfloat: The confidence level used.

Return type:

Bunch

Examples

>>> from process_improve.multivariate import PLS, MCUVScaler
>>> scaler_x = MCUVScaler().fit(X)
>>> scaler_y = MCUVScaler().fit(Y)
>>> X_s, Y_s = scaler_x.transform(X), scaler_y.transform(Y)
>>> pls = PLS(n_components=2).fit(X_s, Y_s)
>>> cv_results = pls.cross_validate(X_s, Y_s, cv="loo")
>>> cv_results.beta_mean          # mean beta across LOO resamples
>>> cv_results.significant        # which betas are significantly != 0
>>> cv_results.q_squared          # cross-validated R²

prediction_interval(X, *, conf_level=0.95, cv_result=None)[source]#

Prediction interval for the Y predictions of new observations.

The interval combines the residual error variance with the leverage of each new observation in the latent-variable space. For a new observation the prediction-interval half-width on target m is

t * s_E[m] * sqrt(1 + 1/N + T2_new / (N - 1))

where s_E is the residual error standard deviation, T2_new is the Hotelling’s T² of the new observation, N is the number of calibration samples, and t is the Student-t quantile.

Parameters:

X (array-like of shape (n_new, n_features)) – New observations, pre-processed the same way as the training data.
conf_level (float, default=0.95) – Confidence level for the interval, in (0.5, 1.0).
cv_result (sklearn.utils.Bunch or None, default=None) – The result of cross_validate(). When supplied, its cross-validated RMSE (rmse_cv) is used for the error variance, which is preferable to the optimistic calibration RMSE used otherwise.

Returns:

With keys y_hat (point predictions), lower and upper (prediction-interval bounds) - each a DataFrame of shape (n_new, n_targets) - and conf_level.

Return type:

set_fit_request(*, sample_weight='$UNCHANGED$')#

Configure whether metadata should be requested to be passed to the fit method.

Note that this method is only relevant when this estimator is used as a sub-estimator within a meta-estimator and metadata routing is enabled with enable_metadata_routing=True (see sklearn.set_config()). Please check the User Guide on how the routing mechanism works.

The options for each parameter are:

True: metadata is requested, and passed to fit if provided. The request is ignored if metadata is not provided.
False: metadata is not requested and the meta-estimator will not pass it to fit.
None: metadata is not requested, and the meta-estimator will raise an error if the user provides it.
str: metadata should be passed to the meta-estimator with this given alias instead of the original name.

The default (sklearn.utils.metadata_routing.UNCHANGED) retains the existing request. This allows you to change the request for some parameters and not others.

Added in version 1.3.

Parameters:

sample_weight (str, True, False, or None, default=sklearn.utils.metadata_routing.UNCHANGED) – Metadata routing for sample_weight parameter in fit.
self (PLS)

Returns:

self – The updated object.

Return type:

set_score_request(*, sample_weight='$UNCHANGED$')#

Configure whether metadata should be requested to be passed to the score method.

The options for each parameter are:

True: metadata is requested, and passed to score if provided. The request is ignored if metadata is not provided.
False: metadata is not requested and the meta-estimator will not pass it to score.
None: metadata is not requested, and the meta-estimator will raise an error if the user provides it.
str: metadata should be passed to the meta-estimator with this given alias instead of the original name.

The default (sklearn.utils.metadata_routing.UNCHANGED) retains the existing request. This allows you to change the request for some parameters and not others.

Added in version 1.3.

Parameters:

sample_weight (str, True, False, or None, default=sklearn.utils.metadata_routing.UNCHANGED) – Metadata routing for sample_weight parameter in score.
self (PLS)

Returns:

self – The updated object.

Return type:

PLS-DA#

PLS discriminant analysis: PLS regression against a one-hot class indicator, with the decision rule, the classifier diagnostics and the permutation test on top. Everything PLS offers is inherited, so a fitted PLSDA also has scores, loadings, VIP, Hotelling’s T2 and SPE.

class process_improve.multivariate.methods.PLSDA(n_components, *, decision_rule='max', priors='empirical', scale=True, max_iter=1000, tol=np.float64(1.4901161193847656e-08), copy=True, missing_data_settings=None, warn_on_uncentred=True)[source]#

Bases: ClassifierMixin, PLS

PLS discriminant analysis: PLS regression against a class-indicator matrix.

The labels are one-hot encoded into an N x G indicator Y, a PLS model is fitted on it, and a new sample’s class comes from the G predicted indicator values. Two decision rules are offered, and they answer different questions:

"max"

Take the largest predicted indicator. Every sample is assigned to exactly one class. This is the usual default and is what most PLS-DA software does.

"bayes"

Take the largest class posterior, built from Gaussians fitted to the training indicator values of each class (in-class and out-of-class) and weighted by the class priors.

This is the rule to prefer when the classes are badly unbalanced, and the reason is not the one that first suggests itself. A rare class’s indicator column is pulled toward zero, because nine rows in ten want it there, so "max" hands almost everything to the common class: on a 1:9 fixture in the test suite it finds two of eight rare samples while reporting 92.5% accuracy, which is the classic accuracy trap. "bayes" reads each column against its own in-class and out-of-class spread instead of against the other columns, so a value that is high for the rare column still counts. On the same fixture it finds seven of eight, and comes out ahead on accuracy too.

Both rules return exactly one label per sample. The per-class Bayesian thresholds, which do allow “in no class” and “in more than one class” readings, are exposed separately as thresholds_.

Parameters:

n_components (int) – Number of latent variables. As for PLS, more is not better: PLS-DA on wide data will separate anything given enough components, which is what permutation_test() exists to check.
decision_rule ({"max", "bayes"}, optional) – Which rule predict() applies. Default is "max".
priors ("empirical", "uniform", or array-like of shape (n_classes,), optional) – Class priors used by the "bayes" rule and by thresholds_. "empirical" (the default) takes them from the training class frequencies; "uniform" gives every class 1 / G, which is what you want when the training set was deliberately balanced but the population is not.
scale (bool, optional) – Passed to PLS. Default True, which mean-centres and unit-variance-scales both blocks. Scaling the indicator block is standard for PLS-DA: without it a rare class contributes less variance and is fitted less well.
max_iter (int) – Passed through to PLS unchanged.
tol (float) – Passed through to PLS unchanged.
copy (bool) – Passed through to PLS unchanged.
missing_data_settings (dict | None) – Passed through to PLS unchanged.
warn_on_uncentred (bool) – Passed through to PLS unchanged.

classes_#

Sorted unique labels seen in fit, in the column order of everything below.

Type:: np.ndarray of shape (n_classes,)

n_classes_#

len(classes_).

Type:: int

priors_#

The resolved priors, summing to 1.

Type:: np.ndarray of shape (n_classes,)

thresholds_#

The per-class Bayesian decision threshold on the predicted indicator; see _bayes_threshold().

Type:: pd.Series indexed by classes_

class_statistics_#

One row per class, columns mean_in, sd_in, mean_out, sd_out, prior, threshold: the fitted Gaussians the "bayes" rule uses.

Type:: pd.DataFrame

confusion_matrix_#

Training-set confusion matrix, rows the true class and columns the predicted one.

Type:: pd.DataFrame

accuracy_#

Training-set accuracy. Optimistic by construction; use score() on held-out data, or permutation_test(), to learn anything about generalisation.

Type:: float

sensitivity_, specificity_

Per-class true-positive and true-negative rates on the training set.

Type:: pd.Series indexed by classes_

Every fitted attribute of :class:`PLS` is also present (``scores_``, ``x_loadings_``,

``x_weights_``, ``r2_cumulative_``, ``hotellings_t2_``, ...), along with its

convenience methods (``score_plot``, ``loading_plot``, ``vip``, ``spe_limit``,

``t2_contributions``, ...), because this class is a :class:`PLS`.

Examples

>>> model = PLSDA(n_components=2).fit(X, labels)
>>> model.predict(X_new)
array(['good', 'bad', 'good'], dtype=object)
>>> model.confusion_matrix_
>>> model.permutation_test(X, labels, n_permutations=99).p_value

References

Barker, M. & Rayens, W. (2003). Partial least squares for discrimination. Journal of Chemometrics 17:166-173.

Brereton, R.G. & Lloyd, G.R. (2014). Partial least squares discriminant analysis: taking the magic away. Journal of Chemometrics 28:213-225.

Westerhuis, J.A. et al. (2008). Assessment of PLSDA cross validation. Metabolomics 4:81-89.

confusion_matrix_plot(matrix=None, settings=None, fig=None)#

Generate a confusion-matrix heat map for a fitted PLSDA model.

Rows are the true class, columns the predicted one, so the diagonal is what the model got right and every off-diagonal cell names a specific confusion: which class this one is mistaken for, which is the question a classification report cannot answer.

Parameters:

model (PLSDA object) – A fitted PLS-DA model generated by this library.
matrix (pd.DataFrame, optional) – A confusion matrix to plot instead of the model’s training-set one, indexed and labelled by class. Pass model.confusion(X_test, y_test).matrix to see the held-out picture, which is the one worth acting on: confusion_matrix_ is fitted on the same rows it is scored on and will always look better.

settings (dict) –

Default settings:

{
    "normalize": False,             # bool: show row fractions, not counts
    "title": "Confusion matrix",    # str: overall plot title
    "colorscale": "Blues",          # str: any Plotly colorscale name
    "show_values": True,            # bool: print the value in each cell
    "html_image_height": 500,       # int: image height in pixels
    "html_aspect_ratio_w_over_h": 1.0,  # float: width as ratio of height
    "template": "pi_journal",       # str: registered Plotly theme name
}

fig (go.Figure, optional) – An existing figure to draw onto. A new figure is created if omitted.

Returns:

fig

Return type:

go.Figure

Raises:

ValueError – If the model is not fitted and no matrix is supplied.

Examples

>>> model.confusion_matrix_plot()
>>> held_out = model.confusion(X_test, y_test).matrix
>>> model.confusion_matrix_plot(held_out, {"normalize": True})

fit(X, y, sample_weight=None)[source]#

Fit the PLS model on a one-hot encoding of y, then the decision layer.

Parameters:

X (array-like of shape (n_samples, n_features)) – Training data. May contain missing values, as for PLS.
y (array-like of shape (n_samples,)) – Class labels: strings, integers, or anything numpy.unique() can sort.
sample_weight (array-like of shape (n_samples,), optional) – Passed through to PLS.fit().

Returns:

self

Return type:

PLSDA

Raises:

ValueError – If y holds fewer than two distinct labels (there is nothing to discriminate), if its length does not match X, or if decision_rule / priors is not a recognised value.

decision_function(X)[source]#

Predicted class-indicator values, one column per class.

These are the raw PLS predictions of the one-hot Y, so they are centred near 1 for the class a sample belongs to and near 0 for the others, but they are not bounded to [0, 1] and do not sum to 1. Use predict_proba() for a normalised quantity.

Parameters:: X (array-like of shape (n_samples, n_features))
Returns:: scores – Indexed by X’s rows, columns in classes_ order.
Return type:: pd.DataFrame of shape (n_samples, n_classes)

predict_proba(X)[source]#

Class posteriors from the fitted per-class Gaussians.

In column g the training indicator values of class g’s members are modelled as N(mean_in, sd_in) and everyone else’s as N(mean_out, sd_out). The evidence for class g is prior_g times the ratio of those two densities, and the posteriors are those normalised across the classes. This is a genuine probability model rather than a rescaling of the indicators: it accounts for how tightly each class scores, for how well it separates from the rest, and for how common it is.

It is still only as good as the Gaussian assumption, which a bimodal class or one with three training members will not satisfy. The raw quantity the model actually computes is decision_function().

Parameters:: X (array-like of shape (n_samples, n_features))
Returns:: proba – Rows sum to 1; columns in classes_ order, as sklearn requires.
Return type:: np.ndarray of shape (n_samples, n_classes)

predict(X)[source]#

Predict a class label for every row of X.

The return type deliberately narrows PLS.predict(), which hands back a DataFrame of predicted Y values: a classifier’s predict returns labels, and sklearn’s classifier contract requires exactly that. The DataFrame that PLS.predict() would have returned is still available, as decision_function().

Parameters:: X (array-like of shape (n_samples, n_features))
Returns:: labels – Values drawn from classes_, chosen by decision_rule.
Return type:: np.ndarray of shape (n_samples,)

See also

decision_function: the indicator values the rule is applied to.
predict_proba: class posteriors.

score(X, y, sample_weight=None)[source]#

Return the mean accuracy on X against the true labels y.

Overrides PLS.score(), which returns R2 of the indicator predictions: a number that goes up when the indicators are fitted more tightly, not when more samples land in the right class. sklearn’s convention (higher is better) holds either way, but only accuracy answers the question a classifier is asked.

Parameters:

X (ndarray | DataFrame)
y (ndarray | DataFrame)
sample_weight (ndarray | None)

Return type:

confusion(X, y)[source]#

Confusion matrix and per-class rates on data the caller supplies.

The confusion_matrix_ / sensitivity_ / specificity_ attributes are the training-set versions and are optimistic; this is the same calculation on a held-out split.

Parameters:

X (array-like of shape (n_samples, n_features))
y (array-like of shape (n_samples,)) – True labels. Labels outside classes_ raise.

Returns:

result – matrix (DataFrame), sensitivity and specificity (Series indexed by class), and accuracy (float).

Return type:

Raises:

ValueError – If y holds a label the model was not fitted on.

roc_auc(X, y, *, positive_class=None)[source]#

Area under the ROC curve, from the continuous indicator rather than the label.

Parameters:

X (array-like of shape (n_samples, n_features))
y (array-like of shape (n_samples,)) – True labels.
positive_class (object, optional) – Which class counts as positive. Required when there are more than two classes, where the result is that class against all the others; for two classes it defaults to classes_[1].

Returns:

auc – Computed on decision_function(), not on the hard labels, so it measures the ranking the model produces and is unaffected by the decision rule.

Return type:

Raises:

ValueError – If positive_class is omitted with more than two classes, or names a class the model was not fitted on.

permutation_test(X, y, *, n_permutations=99, cv=5, random_state=None)[source]#

Test whether the model separates the classes better than shuffled labels do.

PLS-DA on wide data will separate almost anything: with more variables than samples there is always a direction that happens to line up with the labels. The permutation test is what tells a real effect from that, by refitting the same model on shuffled labels and asking how often chance does as well.

Parameters:

X (array-like of shape (n_samples, n_features))
y (array-like of shape (n_samples,)) – True labels.
n_permutations (int, optional) – Number of label shuffles. Default 99, which puts the smallest attainable p-value at 1 / 100.
cv (int, splitter, or None, optional) –
Cross-validation for the accuracy that is compared. Default 5, stratified. None uses training accuracy, which is far quicker and actively misleading. Measured on 40 samples of 30 pure-noise variables with two random labels, 49 permutations, four folds:

Statistic

Separable

Pure noise

Training accuracy

1.000

1.000

Cross-validated

1.000

0.575

p, cv=4

0.020

0.220

p, cv=None

0.020

0.020

The noise column is the whole argument. With more variables than samples the model fits random labels perfectly, so training accuracy says 1.000 either way and the cross-validated test correctly declines to call it significant, while the training-accuracy test reports p = 0.02 on data that has no signal in it at all. Westerhuis et al. (2008) is explicit that the comparison has to be made on cross-validated performance; cv=None is offered for a quick look and for well-conditioned data, not for a result to report.
random_state (int, np.random.Generator, or None, optional) – Seeds the shuffles, per the reproducibility contract.

Returns:

Statistic	Separable	Pure noise
Training accuracy	1.000	1.000
Cross-validated	1.000	0.575
p, `cv=4`	0.020	0.220
p, `cv=None`	0.020	0.020

result – observed (float), null (np.ndarray of accuracies), p_value (float), n_permutations (int, the number of draws the null actually holds) and n_failed (int).

A shuffled label set can leave a fold with nothing to fit at this component count, and NIPALS then goes singular. Those draws are dropped rather than scored as zero, which would push the null down and the p-value with it. When any are dropped a SpecificationWarning says how many, because a null that keeps collapsing means the component count is too high for the data rather than that the model is good.

The p-value is (1 + #{null >= observed}) / (1 + n_permutations). The observed statistic counts itself among the permutations, because no finite set of shuffles licenses a claim of exactly zero; the same convention is used by the Van der Voet and multiblock randomization tests in this package.

Return type:

set_fit_request(*, sample_weight='$UNCHANGED$')#

Configure whether metadata should be requested to be passed to the fit method.

The options for each parameter are:

True: metadata is requested, and passed to fit if provided. The request is ignored if metadata is not provided.
False: metadata is not requested and the meta-estimator will not pass it to fit.
None: metadata is not requested, and the meta-estimator will raise an error if the user provides it.
str: metadata should be passed to the meta-estimator with this given alias instead of the original name.

The default (sklearn.utils.metadata_routing.UNCHANGED) retains the existing request. This allows you to change the request for some parameters and not others.

Added in version 1.3.

Parameters:

sample_weight (str, True, False, or None, default=sklearn.utils.metadata_routing.UNCHANGED) – Metadata routing for sample_weight parameter in fit.
self (PLSDA)

Returns:

self – The updated object.

Return type:

set_score_request(*, sample_weight='$UNCHANGED$')#

Configure whether metadata should be requested to be passed to the score method.

The options for each parameter are:

True: metadata is requested, and passed to score if provided. The request is ignored if metadata is not provided.
False: metadata is not requested and the meta-estimator will not pass it to score.
None: metadata is not requested, and the meta-estimator will raise an error if the user provides it.
str: metadata should be passed to the meta-estimator with this given alias instead of the original name.

The default (sklearn.utils.metadata_routing.UNCHANGED) retains the existing request. This allows you to change the request for some parameters and not others.

Added in version 1.3.

Parameters:

sample_weight (str, True, False, or None, default=sklearn.utils.metadata_routing.UNCHANGED) – Metadata routing for sample_weight parameter in score.
self (PLSDA)

Returns:

self – The updated object.

Return type:

PRM#

Partial Robust M-regression: PLS with a bounded influence per observation, so that a handful of outlying rows move the fit a little instead of a lot. Each row carries the product of a residual weight and a leverage weight, and the fit and the weights are recomputed from each other until they settle. Everything PLS offers is inherited, so a fitted PRM also has scores, loadings, VIP, Hotelling’s T2 and SPE.

class process_improve.multivariate.methods.PRM(n_components, *, cutoff=4.0, max_weight_iter=100, weight_tol=0.0001, scale=True, max_iter=1000, tol=np.float64(1.4901161193847656e-08), copy=True, missing_data_settings=None, warn_on_uncentred=True)[source]#

Bases: PLS

Partial Robust M-regression: PLS with a bounded influence per observation.

Each row carries a weight $w_i = w_i^r \cdot w_i^x$, the product of a residual weight (how badly the current model predicts that row) and a leverage weight (how far its scores sit from the middle of the score cloud). Both come from the Fair function, so the two kinds of outlier that break a least-squares fit are handled by the same mechanism:

a vertical outlier has an ordinary position in X but a y that does not follow the relationship, and gets a small $w^r$;
a bad leverage point sits far out in X as well, where least squares gives it the most influence of all, and gets a small $w^x$.

The weights are recomputed from the fit and the fit is recomputed from the weights until they settle.

Note

Centring and scaling are weighted too, which is not a refinement but the part that makes the method work. The column mean and standard deviation have a breakdown point of zero, so leaving them unweighted leaves the outliers setting the coordinate system that the weighted fit then runs in: measured on a fixture with 15% vertical outliers, weighting only the fit recovers none of the damage, and weighting the scaling as well recovers essentially all of it. See _make_scalers().

Parameters:

n_components (int) – Number of components to extract.
cutoff (float, optional) – Tuning constant $c$ of the Fair weight function, default 4.0 as recommended by Serneels et al. Smaller is more aggressive: an observation at $c$ standardised units gets weight 0.25. This trades robustness against efficiency, and 4.0 keeps roughly 95% of the efficiency of ordinary PLS on clean Gaussian data.
max_weight_iter (int, optional) – Maximum reweighting iterations, default 100. Each one is a full PLS fit.
weight_tol (float, optional) – Convergence tolerance, default 1e-4, on the largest absolute change in any row’s weight between iterations.
scale (bool) – As for PLS.
max_iter (int) – As for PLS.
tol (float) – As for PLS.
copy (bool) – As for PLS.
missing_data_settings (dict | None) – As for PLS.
warn_on_uncentred (bool) – As for PLS.

robust_weights_#

The converged row weights, in (0, 1]. Small entries are this model’s statement about which rows it declined to be led by; see outlier_summary().

Type:: np.ndarray of shape (n_samples,)

n_weight_iter_#

Reweighting iterations actually run.

Type:: int

weights_converged_#

Whether the loop met weight_tol. False is not necessarily a failure: read weight_shift_ before treating it as one.

Type:: bool

weight_shift_#

The largest change in any row’s weight on the final iteration. This is what makes weights_converged_ = False actionable rather than a bare flag, because the reweighting can settle into a small limit cycle rather than a point; a shift of 0.03 says the weights are stable to 3%, which for most purposes is settled.

Type:: float

All the fitted attributes of :class:`~process_improve.multivariate.methods.PLS`

are also present, and mean the same thing, computed on the final weighted fit.

References

S. Serneels, C. Croux, P. Filzmoser and P.J. Van Espen, “Partial Robust M-regression”, Chemometrics and Intelligent Laboratory Systems, 79 (2005), 55-64.

Examples

>>> model = PRM(n_components=2).fit(X, y)
>>> model.outlier_summary(threshold=0.1)
>>> model.predict(X_new)

fit(X, Y, sample_weight=None)[source]#

Fit by alternating between a weighted PLS fit and a reweighting step.

Parameters:

X (array-like of shape (n_samples, n_features)) – Training data.
Y (array-like of shape (n_samples, n_targets)) – Target values.
sample_weight (array-like of shape (n_samples,), optional) – Prior row weights, multiplied into the robust weights at every iteration. Use this to express knowledge the data cannot carry, such as a run you already know was compromised; it is not needed to downweight outliers, which is what the model is for.

Returns:

self, fitted.

Return type:

PRM

Raises:

ValueError – If cutoff, max_weight_iter or weight_tol is out of range, or if the data contain missing values (see Notes).

Notes

Missing data are refused rather than threaded through. The weights are built from residual and leverage distances, and a row with missing cells has a distance that is not comparable with a complete row’s, so it would be downweighted for being incomplete rather than for being wrong. The underlying PLS does handle missing data; use it, or impute first.

outlier_summary(threshold=0.1)[source]#

Rows the fit declined to be led by, weakest weight first.

Parameters:: threshold (float, optional) – Report rows whose final weight is below this, default 0.1. There is no distinguished value: the weights are continuous by design, so this is a reading aid and not a test.
Returns:: Indexed as the training X was, with one weight column. Empty if nothing falls below threshold.
Return type:: pd.DataFrame
Raises:: AttributeError – If the model is not fitted.

set_fit_request(*, sample_weight='$UNCHANGED$')#

Configure whether metadata should be requested to be passed to the fit method.

The options for each parameter are:

True: metadata is requested, and passed to fit if provided. The request is ignored if metadata is not provided.
False: metadata is not requested and the meta-estimator will not pass it to fit.
None: metadata is not requested, and the meta-estimator will raise an error if the user provides it.
str: metadata should be passed to the meta-estimator with this given alias instead of the original name.

The default (sklearn.utils.metadata_routing.UNCHANGED) retains the existing request. This allows you to change the request for some parameters and not others.

Added in version 1.3.

Parameters:

sample_weight (str, True, False, or None, default=sklearn.utils.metadata_routing.UNCHANGED) – Metadata routing for sample_weight parameter in fit.
self (PRM)

Returns:

self – The updated object.

Return type:

set_score_request(*, sample_weight='$UNCHANGED$')#

Configure whether metadata should be requested to be passed to the score method.

The options for each parameter are:

True: metadata is requested, and passed to score if provided. The request is ignored if metadata is not provided.
False: metadata is not requested and the meta-estimator will not pass it to score.
None: metadata is not requested, and the meta-estimator will raise an error if the user provides it.
str: metadata should be passed to the meta-estimator with this given alias instead of the original name.

The default (sklearn.utils.metadata_routing.UNCHANGED) retains the existing request. This allows you to change the request for some parameters and not others.

Added in version 1.3.

Parameters:

sample_weight (str, True, False, or None, default=sklearn.utils.metadata_routing.UNCHANGED) – Metadata routing for sample_weight parameter in score.
self (PRM)

Returns:

self – The updated object.

Return type:

TPLS#

class process_improve.multivariate.methods.TPLS(n_components, d_matrix, max_iter=500, skip_f_matrix_preprocessing=False, tol=np.float64(1.4901161193847656e-08))[source]#

Bases: 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 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. 3. Conditions in Z may be missing (turning it into an L-shaped data structure).

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.
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.
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 is listed last so that no existing positional argument changed position.

Notes

The input X passed to fit() and predict() 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.
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.

n_samples#

Number of samples (rows) in the training data.

Type:: int

n_substances#

Number of materials (columns) in the F matrix, summed across groups.

Type:: int

n_conditions#

Number of process-condition columns summed across Z blocks.

Type:: int

n_outputs#

Number of quality-indicator columns summed across Y blocks.

Type:: int

is_fitted_#

Set to True once fit() completes.

Type:: bool

fitting_statistics#

Per-component iterations, convergance_tolerance and milliseconds lists.

Type:: dict

preproc_#

Nested per-block per-group preprocessors, indexed as preproc_[block][group] where block is "D", "F", "Y" or "Z".

Type:: dict

sums_of_squares_#

Per-component per-block sums of squares.

Type:: list[dict]

r2_frac#

Per-component per-block cumulative R^2 fractions.

Type:: list[dict]

feature_importance#

Per-block per-group variable-importance (populated for "D" and "F").

Type:: dict

d_mats, f_mats, z_mats, y_mats

Deflated block matrices (per group for D/F, per Z-block / Y-block for Z/Y).

Type:: dict[str, np.ndarray]

not_na_d, not_na_f, not_na_z, not_na_y

Boolean observed-value masks matching the shapes above.

Type:: dict[str, np.ndarray]

observation_names#

Row index shared across the F/Z/Y blocks.

Type:: pd.Index

property_names, material_names

Column and row names of each D group.

Type:: dict[str, list[str]]

condition_names, quality_names

Column names of each Z-block and Y-block respectively.

Type:: dict[str, list[str]]

t_scores_super#

Super-block scores T.

Type:: pd.DataFrame, shape (n_samples, n_components)

r_loadings_f#

F-block weights per group.

Type:: dict[str, pd.DataFrame]

w_loadings_z#

Z-block weights per Z-block.

Type:: dict[str, pd.DataFrame]

w_loadings_super#

Super-block weights, rows ["Z", "F"] (or just ["F"] when there are no process conditions).

Type:: pd.DataFrame

s_loadings_d, v_loadings_d

D-block scores and loadings per group.

Type:: dict[str, pd.DataFrame]

p_loadings_f#

F-block loadings per group (used for deflation).

Type:: dict[str, pd.DataFrame]

p_loadings_z#

Z-block loadings per Z-block.

Type:: dict[str, pd.DataFrame]

q_loadings_y#

Y-block loadings per Y-block.

Type:: dict[str, pd.DataFrame]

hat_#

Per-Y-block predictions on the preprocessed (centred / scaled) scale.

Type:: dict[str, pd.DataFrame]

hat#

Per-Y-block in-sample predictions on the original scale.

Type:: dict[str, pd.DataFrame]

spe#

Nested spe[block][group] squared prediction error tables.

Type:: dict[str, dict[str, pd.DataFrame or pd.Series]]

spe_limit#

Nested spe_limit[block][group] callables. Each is a functools.partial() over 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.

Type:: dict[str, dict[str, Callable]]

hotellings_t2#

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 diagnose() (#502).

Type:: pd.DataFrame, shape (n_samples, n_components)

scaling_factor_for_scores#

Standard deviation of each super-score column, computed with the unbiased (N - 1) divisor; used by the ellipse / T^2 helpers.

Type:: pd.Series

.. 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 __init__ and refined during fit (for example is_fitted_, preproc_, required_blocks_, required_inputs_) 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))

property tolerance_: float#

Deprecated alias for tol, kept for one deprecation cycle.

It was assigned in __init__ and read at two unrelated sites: the convergence test, which tol now owns, and a zero-variance column test, which has its own floor (_ZERO_VARIANCE_FLOOR).

hotellings_t2_limit(conf_level=0.95)[source]#

Hotelling’s T2 limit at the given confidence level (see hotellings_t2_limit()).

Parameters:: conf_level (float)
Return type:: float

ellipse_coordinates(score_horiz, score_vert, conf_level=0.95, n_points=100)[source]#

Coordinates of the T2 confidence ellipse (see ellipse_coordinates()).

Parameters:

score_horiz (int)
score_vert (int)
conf_level (float)
n_points (int)

Return type:

tuple[ndarray, ndarray]

fit(X, y=None)[source]#

Fit the preprocessing parameters and also the latent variable model from the training data.

Parameters:

X ({dictionary of dataframes}, keys that must be present: "F", "Z", and "Y") – The training input samples. See documentation in the class definition for more information on each matrix.
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 – Returns self.

Return type:

Raises:

ValueError – If y is not None.
TypeError – If X is not a DataFrameDict.

predict(X)[source]#

Forward to diagnose(); emits a DeprecationWarning.

Deprecated since version 1.38.4: Use TPLS.diagnose() instead. predict matches the sklearn-convention name (regression-style ndarray return), but TPLS’s natural return is the rich per-group / per-block diagnostics Bunch and TPLS cannot be placed in a standard sklearn Pipeline anyway (its input is a nested DataFrameDict). The rename aligns with PLS.diagnose() and PCA.diagnose(). Will be removed in 2.0.0.

Parameters:: X (DataFrameDict)
Return type:: Bunch

diagnose(X)[source]#

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: new_data = {“Z”: …, “F”: …} # you need at least the F block for a new prediction. “Z” is optional. predictions = estimator.diagnose(new_data)

Parameters:

X (DataFrameDict) – The input samples.

Returns:

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).
hotellings_t2 (DataFrame, shape (n_new, n_components)) : cumulative Hotelling’s T² per new observation after each component.

Return type:

display_results(show_cumulative_stats=True)[source]#

Display the results of the model fitting.

Parameters:: show_cumulative_stats (bool)
Return type:: str

score(X, y=None, sample_weight=None)[source]#

Return the mean r2_score() across Y blocks 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 r2_score().

Returns:

score – The mean of r2_score() computed separately on every Y block in X["Y"], using self.diagnose(X).hat as the prediction. A single Y block gives its own $R^2$; more than one is an unweighted average across blocks (not a pooled multiblock $R^2$).

Return type:

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 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.

help()[source]#

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]

Return type:: str

vip(block=None, method='vip')[source]#

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 feature_importance.

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 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 feature_importance.

Returns:

If block is None: {"D": {group: pd.Series, ...}, "F": {group: pd.Series, ...}}. If block is "D" or "F": the inner dict {group: pd.Series, ...} for that block, where each pd.Series is indexed by feature names.

Return type:

dict

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

property d_block_scaling_: 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:: Mapping of D-block group name to its scalar block-scaling factor.
Return type:: dict[str, float]
Raises:: AttributeError – If the model has not been fitted yet.

set_score_request(*, sample_weight='$UNCHANGED$')#

Configure whether metadata should be requested to be passed to the score method.

The options for each parameter are:

True: metadata is requested, and passed to score if provided. The request is ignored if metadata is not provided.
False: metadata is not requested and the meta-estimator will not pass it to score.
None: metadata is not requested, and the meta-estimator will raise an error if the user provides it.
str: metadata should be passed to the meta-estimator with this given alias instead of the original name.

The default (sklearn.utils.metadata_routing.UNCHANGED) retains the existing request. This allows you to change the request for some parameters and not others.

Added in version 1.3.

Parameters:

sample_weight (str, True, False, or None, default=sklearn.utils.metadata_routing.UNCHANGED) – Metadata routing for sample_weight parameter in score.
self (TPLS)

Returns:

self – The updated object.

Return type:

A T-shaped model carries its response inside X["Y"], so cross_val_score() is called without a y and a scorer string such as scoring="r2" cannot be honoured: sklearn’s _Scorer needs a y_true it was never given, the call fails before TPLS is reached, and every fold is recorded as NaN. Build the scorer with make_tpls_scorer() instead; sklearn passes a callable scoring= through untouched.

process_improve.multivariate.methods.make_tpls_scorer(metric='r2', *, greater_is_better=True, **metric_kwargs)[source]#

Build a scoring= callable that works with 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 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 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 – A callable with signature scorer(estimator, X, y=None, sample_weight=None), usable as the scoring= argument of cross_val_score(), cross_validate() and the *SearchCV classes. Several Y blocks are combined the way TPLS.score() combines them: an unweighted mean over the blocks in X["Y"], not a pooled multiblock statistic.

Return type:

Callable[…, float]

Raises:

ValueError – If metric is a string that is not a known metric name. The returned scorer raises 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)
>>> cross_val_score(model, X=blocks, cv=5, scoring=scorer)
array([-1.07, -0.98, -1.11, -1.02, -1.05])

ASCA#

ANOVA-Simultaneous Component Analysis: partition a response matrix by its design terms, then give each term’s effect matrix its own PCA. This is the bridge between process_improve.experiments and the latent-variable models: it answers which factor owns which direction of multivariate variation, whether that is more than chance, and which variables carry it.

class process_improve.multivariate.methods.ASCA(n_components=2, *, model='interactions', add_residuals=False, scale=False)[source]#

Bases: BaseEstimator

ANOVA-Simultaneous Component Analysis: PCA of each design term’s effect.

The response matrix is written as

\[X = \mathbf{1} m^T + \sum_f X_f + E\]

with $m$ the grand mean, one effect matrix $X_f$ per design term, and $E$ the residual. Each $X_f$ then gets its own PCA, so a term’s scores and loadings describe the multivariate structure of that factor’s effect alone.

Parameters:

n_components (int, optional) – Components to extract per term. Capped per term at the rank its design columns can support, because an effect matrix for a two-level factor has rank 1 and asking for two components of it is asking for a component that does not exist. Default 2.
model (str, optional) – "interactions" (the default) adds every two-way interaction to the main effects; "main_effects" fits main effects only. Anything else is taken as a patsy right-hand side and used verbatim, in which case the caller is responsible for the coding.
add_residuals (bool, optional) – If True, add the residual matrix back onto each effect matrix before its PCA (APCA / ASCA+). The effect matrix alone has one distinct row per factor-level combination, so a score plot of it is a handful of points; adding the residuals back restores the scatter and shows whether the levels actually separate. Default False.
scale (bool, optional) – If True, unit-variance scale the columns of X after centring, so a variable measured in large units does not dominate the decomposition. Default False, which is the right choice when the columns are already on one scale (a spectrum, say) and the wrong one when they are not.

terms_#

Design term names, in the order patsy resolved them, with the coding wrapper stripped: ["A", "B", "A:B"].

Type:: list[str]

grand_mean_#

Column means removed before the decomposition.

Type:: pd.Series

column_scale_#

Column scaling applied, all ones when scale=False.

Type:: pd.Series

effect_matrices_#

One N x K effect matrix per term.

Type:: dict[str, pd.DataFrame]

residuals_#

What no term explains.

Type:: pd.DataFrame

ssq_#

Sum of squares per term, plus "residual" and "total".

Type:: pd.Series

ssq_percent_#

The same as a percentage of the total, which is the “factor effect” summary worth reading first.

Type:: pd.Series

models_#

The fitted per-term PCA, keyed by term.

Type:: dict[str, PCA]

pvalues_#

Permutation p-values per term. Absent until permutation_test() is called.

Type:: pd.Series

is_balanced_#

Whether every factor-level combination occurs equally often.

Type:: bool

Notes

Balance matters, and the model says so rather than assuming it. On a balanced design the terms are orthogonal, the per-term sums of squares add up to the model sum of squares, and the decomposition is unique. On an unbalanced design they are not orthogonal: the split of the shared variation between correlated terms depends on how you choose to attribute it, which is the Type I / II / III question. This implementation fits every term simultaneously by least squares and reads each term’s fitted contribution off that single fit, which is the usual ASCA treatment. It warns when the design is unbalanced, and ssq_ then no longer partitions the total exactly; the gap is reported as ssq_["total"] minus the sum of the parts.

Examples

>>> model = ASCA(n_components=2).fit(X, design)
>>> model.ssq_percent_
>>> model.permutation_test(n_permutations=999, random_state=0)
>>> model.models_["A"].scores_

References

Smilde, A. K., Jansen, J. J., Hoefsloot, H. C. J., Lamers, R.-J. A. N., van der Greef, J., & Timmerman, M. E. (2005). ANOVA-simultaneous component analysis (ASCA): a new tool for analyzing designed metabolomics data. Bioinformatics, 21(13), 3043-3048.

Zwanenburg, G., Hoefsloot, H. C. J., Westerhuis, J. A., Jansen, J. J., & Smilde, A. K. (2011). ANOVA-principal component analysis and ANOVA-simultaneous component analysis: a comparison. J. Chemometrics, 25(10), 561-567.

Camacho, J., Vitale, R., Morales-Jimenez, D., & Gomez-Llorente, C. (2022). Variable-selection ANOVA Simultaneous Component Analysis (VASCA). Bioinformatics, 38(1), 295-298.

effect_summary_plot(settings=None, fig=None)#

Generate the per-term effect summary for a fitted ASCA model.

One bar per design term, showing the share of the total sum of squares it carries, with the residual alongside for scale. This is the plot to read first: it says which factor the variation actually belongs to, before any score plot is opened.

Permutation p-values are annotated on the bars when permutation_test() has been run, because a term’s share and its significance answer different questions: a term can hold a large share simply by having many degrees of freedom.

Parameters:

model (ASCA object) – A fitted ASCA model generated by this library.

settings (dict) –

Default settings:

{
    "include_residual": True,       # bool: draw the residual bar too
    "title": "Variation by design term",  # str: overall plot title
    "bar_color": None,              # str|None: bar colour; None uses the theme
    "html_image_height": 500,       # int: image height in pixels
    "html_aspect_ratio_w_over_h": 16/9,  # float: width as ratio of height
    "template": "pi_journal",       # str: registered Plotly theme name
}

fig (go.Figure, optional) – An existing figure to draw onto. A new figure is created if omitted.

Returns:

fig

Return type:

go.Figure

Raises:

ValueError – If the model is not fitted.

Examples

>>> model.effect_summary_plot()
>>> model.permutation_test(random_state=0)
>>> model.effect_summary_plot()   # now annotated with p-values

fit(X, design, y=None)[source]#

Decompose X by the design terms and fit a PCA to each term’s effect.

Parameters:

X (array-like of shape (n_samples, n_features)) – The multivariate response. Must be complete: ASCA solves a least-squares problem over every column at once and has no missing-data path.
design (pd.DataFrame) – One row per sample and one column per experimental factor, holding that factor’s level for each row. Values may be strings or numbers; they are treated as categorical.
y (object, optional) – Ignored, for sklearn API compatibility.

Returns:

self

Return type:

ASCA

Raises:

ValueError – If X and design disagree on the number of rows, if either holds missing values, or if a factor has only one level (it can carry no effect).

permutation_test(*, n_permutations=999, random_state=None)[source]#

Test each term’s effect against the null of exchangeable rows.

For every term the rows of the response are permuted, the decomposition is refitted, and the term’s sum of squares is recomputed. A term whose observed sum of squares sits in the upper tail of that null is carrying more variation than the design’s shape alone would produce.

Parameters:

n_permutations (int, optional) – Number of permutations. Default 999, which puts the smallest attainable p-value at 1 / 1000.
random_state (int, np.random.Generator, or None, optional) – Seeds the permutations, per the reproducibility contract.

Returns:

pvalues – One p-value per term, also stored as pvalues_. Each is (1 + #{null >= observed}) / (1 + n_permutations): the observed statistic counts itself among the permutations, because no finite set of shuffles licenses a claim of exactly zero. The same convention is used by the Van der Voet test, the multiblock randomization test and PLSDA.permutation_test.

Return type:

pd.Series

Notes

One permutation refits every term at once, so the whole test costs n_permutations least-squares solves rather than one per term. The PCA step is not repeated: the statistic is the sum of squares of the effect matrix, which the decomposition produces directly.

vasca(term, *, n_permutations=999, alpha=0.05, random_state=None)[source]#

Variable-selection ASCA: which variables carry this term’s effect.

The ASCA permutation test asks one question of the whole response matrix, so an effect that lives in three variables out of two hundred is diluted by the other hundred and ninety-seven and can fail to register at all. VASCA ranks the variables by their contribution to the term, then tests each nested subset of the top-ranked ones. A subset that contains the effect and little else gives a far smaller p-value than the whole matrix does.

Parameters:

term (str) – Which design term to examine; one of terms_.
n_permutations (int, optional) – Permutations used to build the null for every subset at once. Default 999.
alpha (float, optional) – Target false-discovery rate for the Benjamini-Hochberg step. Default 0.05.
random_state (int, np.random.Generator, or None, optional) – Seeds the permutations.

Returns:

result – table (pd.DataFrame, one row per subset size: the variable added at that step, the subset’s cumulative sum of squares, how many standard deviations it sits above its own null, and its raw and FDR-corrected p-values), selected (list of variable names), ranking (the variables in contribution order) and p_value (the smallest corrected p-value found).

selected is the subset that clears alpha and stands furthest above its own null. The second half of that matters: with a few hundred permutations the smallest attainable p-value is 1 / (1 + n_permutations) and many subset sizes reach it at once, so choosing by p-value alone would return every variable that happened to tie at the floor. The z-score does not tie, and it peaks where the effect is concentrated.

selected is empty when no subset clears alpha, which is the honest answer for a term that carries nothing.

Return type:

Raises:

ValueError – If term is not one of the fitted terms.

Notes

The permutations are shared across subset sizes: one shuffle produces a per-variable sum of squares vector, and every subset’s null statistic is a partial sum of it. The whole walk therefore costs the same n_permutations solves as the single-term test, rather than one run per subset size.

Because one test is made per subset size, the raw p-values are corrected across those tests with benjamini_hochberg(), which controls the false-discovery rate rather than the family-wise error rate.

set_fit_request(*, design='$UNCHANGED$')#

Configure whether metadata should be requested to be passed to the fit method.

The options for each parameter are:

True: metadata is requested, and passed to fit if provided. The request is ignored if metadata is not provided.
False: metadata is not requested and the meta-estimator will not pass it to fit.
None: metadata is not requested, and the meta-estimator will raise an error if the user provides it.
str: metadata should be passed to the meta-estimator with this given alias instead of the original name.

The default (sklearn.utils.metadata_routing.UNCHANGED) retains the existing request. This allows you to change the request for some parameters and not others.

Added in version 1.3.

Parameters:

design (str, True, False, or None, default=sklearn.utils.metadata_routing.UNCHANGED) – Metadata routing for design parameter in fit.
self (ASCA)

Returns:

self – The updated object.

Return type:

MBPLS#

Multi-block PLS in the hierarchical / superblock formulation of Westerhuis, Kourti & MacGregor (1998). Each X-block is preprocessed independently and weighted by 1/sqrt(K_b) before the inner NIPALS loop, so blocks of unequal width contribute fairly to the consensus super-score.

class process_improve.multivariate.methods.MBPLS(n_components, *, max_iter=500, tol=None, algorithm='auto', missing_data_settings=None)[source]#

Bases: _HotellingsT2LimitMixin, RegressorMixin, BaseEstimator

Multi-block PLS (hierarchical / superblock formulation).

Generic multi-block PLS as described by Westerhuis, Kourti & MacGregor (1998) and Westerhuis & Smilde (2001). Each X-block is preprocessed independently (mean-centred and unit-variance scaled), then divided by sqrt(K_b) so that blocks of unequal width contribute fairly to the super-score.

Parameters:

n_components (int) – Number of latent variables to extract.
max_iter (int, default=500) – Maximum NIPALS iterations per latent variable.
tol (float or None, default=None) – Relative convergence tolerance on the change in the Y-block score u: the norm of the change between two successive iterations, divided by the norm of the current u vector (#504). If None, epsqrt (about 1.49e-8) is used, the same default as PCA / PLS / TPLS. The legacy absolute tolerance np.finfo(float).eps ** (6/7) (about 3.8e-14) sits below the floating-point oscillation floor of a relative criterion, so it would never be reached in practice.
algorithm (str) –
Algorithm to use for fitting the model.
- "auto": dense vectorised hierarchical NIPALS when every block (X and Y) is complete; mask-aware NIPALS when any block contains missing values.
- "dense": dense vectorised hierarchical NIPALS. Raises if any block contains missing values.
- "nipals": mask-aware hierarchical NIPALS. Always uses the NaN-tolerant inner-loop primitives, even when the data is complete (slower than "dense" but produces equivalent results).
With missing data the super score of each row is estimated by one masked regression of the whole row onto the stacked block weights, rather than by adding up a score per block. The two are the same number on a complete row, so no fitted model changes; on an incomplete row the pooled form lets every observed cell carry its share, wherever in the blocks it sits, instead of letting a block seen in one variable speak as loudly as a block seen in twenty. A row with nothing observed in one block is therefore still scored, from the blocks it does have; only a row observed in no block at all is refused.
missing_data_settings (dict or None, default=None) – Settings for the iterative "nipals" path. Keys: md_tol (convergence tolerance on the score-vector change between iterations), md_max_iter (maximum NIPALS iterations per component). Defaults to {"md_tol": epsqrt, "md_max_iter": 1000}.
fitting) (Attributes (after)
--------------------------
block_names (list[str]) – Ordered list of X-block names (the keys of the input dict).
block_widths (dict[str, int]) – Number of variables in each X-block.
n_samples (int) – Number of rows fitted.
n_targets (int) – Number of Y columns.
n_features_in (int) – Total number of X variables summed across blocks.
feature_names_in (np.ndarray) – Concatenated column names, one per feature, in block order.
preproc (dict[str, MCUVScaler]) – Per-block preprocessors used to mean-centre and unit-variance scale each X-block.
y_preproc (MCUVScaler) – Preprocessor used on Y.
super_scores (pd.DataFrame, shape (n_samples, n_components)) – Super-block (consensus) X-scores T. Finite for every row that has at least one observed cell, in any block.
super_y_scores (pd.DataFrame, shape (n_samples, n_components)) – Super-block Y-scores U.
super_weights (pd.DataFrame, shape (n_blocks, n_components)) – Super-block weights w_super; rows indexed by block name.
super_y_loadings (pd.DataFrame, shape (n_targets, n_components)) – Y-block loadings c.
super_hotellings_t2 (pd.DataFrame, shape (n_samples, n_components)) – Cumulative Hotelling’s T^2 on the super-scores per component.
super_vip (pd.Series) – Variable-importance in projection for each X-block, indexed by block name.
block_scores (dict[str, pd.DataFrame]) – Per-block X-scores t_b, each shape (n_samples, n_components). NaN for a row with nothing observed in that block: the block has no score of its own there, and reporting zero would place the row at the block’s average instead.
block_weights (dict[str, pd.DataFrame]) – Per-block X-weights w_b, each shape (K_b, n_components). Each column has unit norm.
block_loadings (dict[str, pd.DataFrame]) – Per-block X-loadings p_b (used for deflation), each shape (K_b, n_components).
block_spe (dict[str, pd.DataFrame]) – Per-block squared prediction error per sample and component. NaN where the row has nothing observed in that block, for the same reason as block_scores_.
block_hotellings_t2 (dict[str, pd.DataFrame]) – Per-block cumulative Hotelling’s T^2 per sample and component.
block_vip (dict[str, pd.Series]) – Per-block variable-importance in projection, indexed by variable name inside each block.
predictions (pd.DataFrame, shape (n_samples, n_targets)) – In-sample Y predictions on the original scale.
explained_variance (np.ndarray, shape (n_components,)) – Variance of the super-score per component (ddof=1).
scaling_factor_for_super_scores (pd.Series) – sqrt(explained_variance_) per component.
r2_x_per_block_cumulative (pd.DataFrame, shape (n_blocks, n_components)) – Cumulative R^2X per block and component.
r2_x_per_block_per_component (pd.DataFrame, shape (n_blocks, n_components)) – Incremental R^2X per block and component.
r2_x_per_variable (dict[str, pd.DataFrame]) – Cumulative R^2X per variable within each block.
r2_y_cumulative (pd.Series, shape (n_components,)) – Cumulative R^2Y per component.
r2_y_per_component (pd.Series, shape (n_components,)) – Incremental R^2Y per component.
r2_y_per_variable (pd.DataFrame, shape (n_targets, n_components)) – Cumulative R^2Y per Y-variable and component.
fitting_info (dict) – Per-component iteration count and timing.
has_missing_data (bool) – Whether any X-block or Y had NaN values.
algorithm – The resolved algorithm actually used for the fit. With algorithm="auto", this is "dense" for complete data and "nipals" for NaN-containing data.

Notes

Block weighting uses the convention $X_b / \sqrt{K_b}$ so that every block contributes the same total sum of squares to the super-score, regardless of how many variables it has.

Missing data#

When any X-block or Y contains NaN entries, the "auto" algorithm routes to a mask-aware NIPALS variant. The X-block weights, block scores, block loadings used for deflation, Y-block loadings and Y-block scores are each computed as a regression that uses only the observed entries; the masked sum-of-squares is used as the denominator so missing values neither bias the latent direction nor contribute to the score. The mask is preserved across components automatically because deflation propagates NaN through subtraction. This is the standard skip-NaN NIPALS update; see Walczak & Massart (2001) and Arteaga & Ferrer (2002).

The fit refuses to run if any X-block or Y has a column with all entries missing, or a row with all entries missing for that block; either case leaves the masked denominator at zero. Drop or impute such rows or columns before fitting. Predict-time score estimation for new observations with NaN (Trimmed Score Regression / Projection to the Model Plane) is a separate follow-up.

References

Westerhuis, J. A., Kourti, T. & MacGregor, J. F. Analysis of multiblock and hierarchical PCA and PLS models. Journal of Chemometrics, 12 (1998), 301-321.

Westerhuis, J. A. & Smilde, A. K. Deflation in multiblock PLS. Journal of Chemometrics, 15 (2001), 485-493.

Walczak, B. & Massart, D. L. Dealing with missing data: Part I. Chemom. Intell. Lab. Syst., 58 (2001), 15-27.

Arteaga, F. & Ferrer, A. Dealing with missing data in MSPC: several methods, different interpretations, some examples. J. Chemometrics, 16 (2002), 408-418.

fit(X, y)[source]#

Fit the multi-block PLS model.

Parameters:

X (dict[str, pd.DataFrame]) – X-blocks. Keys are block names; values are DataFrames sharing the same row index (and row count). Each block is preprocessed independently.
y (pd.DataFrame) – Y-block. Same row index / row count as the X-blocks.

Return type:

MBPLS

block_spe_limit(block, conf_level=0.95)[source]#

SPE limit for one X-block using the Nomikos & MacGregor chi-square approximation.

Operates on the same scale as block_spe_[block] (sqrt of row sum of squares), so the value can be drawn directly on a SPE plot.

Parameters:

block (str)
conf_level (float)

Return type:

super_spe_limit(conf_level=0.95)[source]#

SPE limit for the merged super-block (sum of per-block SPE squared).

Parameters:: conf_level (float)
Return type:: float

spe_contributions(X)[source]#

Per-variable squared residuals for each X-block (SPE contributions).

For each new observation and each X-block, reconstruct the block as T_super @ P_b^T (matching the deflation step used during fit) and return the squared per-variable residuals. Useful for fault diagnosis: the variable with the largest contribution is the most likely culprit for a high SPE.

Returns:: One DataFrame per block, shape (n_samples, K_b). Values are preprocessed-scale squared residuals (centred and scaled inside the model). Sum across columns equals block_spe_[b].iloc[:, -1] ** 2.
Return type:: dict[str, pd.DataFrame]
Parameters:: X (dict[str, DataFrame])

score_contributions(X, component=1, scaling='none')[source]#

Per-block per-variable contributions to a super-score.

The multi-block analogue of PLS.score_contributions(). A super score is a weighted sum of the (deflated, preprocessed) variables across every block, so it splits exactly into one term per variable:

\[c_{b,ij}^{(a)} = \tilde{x}_{b,ij}^{(a)}\, \frac{w_b[j, a]\, w_\mathrm{super}[b, a]}{\sqrt{K_b}}, \qquad \sum_b \sum_j c_{b,ij}^{(a)} = t_{\mathrm{super},ia},\]

where $\tilde{x}^{(a)}$ is the block data deflated through the first $a-1$ components, which is what the super score at component $a$ is actually formed from.

Parameters:

X (dict[str, pd.DataFrame]) – Raw (un-preprocessed) X-blocks, keyed by block name, exactly as passed to fit(). The stored per-block preprocessing is applied internally.
component (int, default=1) – 1-based component index whose super score is decomposed.
scaling ({"none", "maximum", "within"}, default="none") – Presentation scaling, as for PLS.score_contributions(). Under "none" the contributions sum across all blocks to the super score. "maximum" divides by the largest absolute contribution over every block; "within" divides each observation by the total absolute contribution it accumulates across every block, so both scalings are taken over the blocks jointly rather than one block at a time.

Returns:

One frame per X-block, of shape (n_samples, K_b).

Return type:

dict[str, pd.DataFrame]

Examples

>>> mbpls = MBPLS(n_components=2).fit(blocks, Y)
>>> contrib = mbpls.score_contributions(blocks, component=1)
>>> sum(frame.sum(axis=1) for frame in contrib.values())  # super score 1

group_contributions(X, group, reference=None, component=1)[source]#

Per-block per-variable contributions to a group’s average super score.

The multi-block analogue of PLS.group_contributions(). See that method for the definition; the only difference is that the result is returned one Series per X-block, and the sum over every block equals the group’s average super score (or the difference between the two groups’ average super scores when reference is given).

Parameters:

X (dict[str, DataFrame])
group (Sequence)
reference (Sequence | None)
component (int)

Return type:

dict[str, Series]

super_score_plot(pc_horiz=1, pc_vert=2)[source]#

Scatter plot of super-scores for two components.

Parameters:

pc_horiz (int)
pc_vert (int)

Return type:

Figure

super_weights_bar_plot(component=1)[source]#

Bar plot of super-weights w_super for a single component.

Parameters:: component (int)
Return type:: Figure

predictions_vs_observed_plot(y_observed, variable=None)[source]#

Scatter plot of predicted vs observed Y, with y=x reference and RMSEE annotation.

Parameters:

y_observed (pd.DataFrame) – The observed Y on the original scale, same columns as the training Y.
variable (str or None, default=None) – If given, plot only that Y-variable. If None, plot the first one.

Return type:

Figure

display_results(show_cumulative=True)[source]#

Format a short text summary of per-block R²X, overall R²Y, iterations and timing.

Parameters:: show_cumulative (bool)
Return type:: str

transform(X)[source]#

Project new data to super-scores using the fitted model.

Parameters:: X (dict[str, DataFrame])
Return type:: DataFrame

diagnose(X)[source]#

Project new data and return the full diagnostics Bunch.

Returns a sklearn.utils.Bunch with fields super_scores (DataFrame, n_samples x n_components), block_scores (dict[str, DataFrame]), predictions (DataFrame on original Y scale), block_spe (dict[str, Series], per-block SPE of the new observations) and hotellings_t2 (Series of cumulative Hotelling’s T² over all components, per new observation).

The rename (since 1.38.4, #395) matches PLS.diagnose() and PCA.diagnose; predict() is kept as a deprecation shim.

Parameters:: X (dict[str, DataFrame])
Return type:: Bunch

classmethod select_n_components(X, y, *, max_components=None, cv=5, n_repeats=None, random_state=None, selection_rule='1se', **mbpls_kwargs)[source]#

Select the number of multi-block PLS components by cross-validation.

Whole rows are held out. The super score of a held-out row is computed from its X-blocks alone and its Y is what the model predicts, so the value being predicted never enters its own prediction. That is the same argument that makes row-wise cross-validation sound for a single-block PLS.select_n_components(), and it is unaffected by there being several X-blocks.

Each block is centred and scaled inside fit(), on the training rows only, so the fold statistics never see the held-out rows.

One model is fitted per fold and per component count, because the hierarchical NIPALS deflation means an a-component model is not recoverable from an A-component one. The cost is cv * n_repeats * max_components fits.

Parameters:

X (dict[str, pd.DataFrame]) – X-blocks, keyed by block name, all sharing y’s row index.
y (pd.DataFrame) – Y-block, one row per observation.
max_components (int, optional) – Largest component count to evaluate. Defaults to the largest the smallest training fold supports, capped at the total width of the X-blocks.
cv (int or sklearn CV splitter, default 5) – An integer is used as the n_splits of a shuffled KFold, or of a RepeatedKFold when n_repeats > 1. A splitter object is used as given, and n_repeats is then ignored.
n_repeats (int, optional) – How many times to repeat the split with a fresh shuffle. Resolved to 10 when cv is an integer; pass 1 to disable repeats.
random_state (int, optional) – Seed for the shuffling. Ignored when cv is a splitter.
selection_rule ({"1se", "min", "q2_increment"}, default "1se") – How n_components is chosen from the curve. See SelectionRule. "randomization" is not offered here.
**mbpls_kwargs – Passed to every MBPLS fitted, for instance tol or algorithm.

Returns:

With n_components (int), rmsecv and se_rmsecv (Series indexed 1..A), per_fold_rmsecv (DataFrame, components by fold), press (Series), r2y_validated (DataFrame with one column per target, plus "total" on the original Y scale and "scaled_total" with every target weighted equally), cv_predictions (DataFrame of the held-out predictions of the recommended model, averaged over repeats) and selection_rule.

Return type:

Raises:

ValueError – If X is not a non-empty dict of frames sharing y’s index, or if no component count could be evaluated.

Examples

>>> import numpy as np, pandas as pd
>>> from process_improve.multivariate.methods import MBPLS
>>> rng = np.random.default_rng(0)
>>> t = rng.standard_normal((40, 2))
>>> blocks = {
...     "a": pd.DataFrame(t @ rng.standard_normal((2, 5)) + rng.standard_normal((40, 5)) * 0.3),
...     "b": pd.DataFrame(t @ rng.standard_normal((2, 4)) + rng.standard_normal((40, 4)) * 0.3),
... }
>>> Y = pd.DataFrame(t @ rng.standard_normal((2, 2)) + rng.standard_normal((40, 2)) * 0.3)
>>> out = MBPLS.select_n_components(blocks, Y, max_components=3, cv=5, n_repeats=2, random_state=0)
>>> 1 <= out.n_components <= 3
True

predict(X)[source]#

Forward to diagnose(); emits a DeprecationWarning.

Deprecated since version 1.38.4: Use MBPLS.diagnose() instead. The sklearn-convention predict name suggests a regression-style ndarray return, but the historical return is the rich diagnostics Bunch. The rename aligns with PLS.diagnose() and PCA.diagnose() and frees the name for a future contract that returns just the predictions field. Will be removed in 2.0.0.

Parameters:: X (dict[str, DataFrame])
Return type:: Bunch

set_score_request(*, sample_weight='$UNCHANGED$')#

Configure whether metadata should be requested to be passed to the score method.

The options for each parameter are:

True: metadata is requested, and passed to score if provided. The request is ignored if metadata is not provided.
False: metadata is not requested and the meta-estimator will not pass it to score.
None: metadata is not requested, and the meta-estimator will raise an error if the user provides it.
str: metadata should be passed to the meta-estimator with this given alias instead of the original name.

The default (sklearn.utils.metadata_routing.UNCHANGED) retains the existing request. This allows you to change the request for some parameters and not others.

Added in version 1.3.

Parameters:

sample_weight (str, True, False, or None, default=sklearn.utils.metadata_routing.UNCHANGED) – Metadata routing for sample_weight parameter in score.
self (MBPLS)

Returns:

self – The updated object.

Return type:

process_improve.multivariate.methods.randomization_test_mbpls(model, X, y, n_permutations=200, *, seed=None)[source]#

Randomization (permutation) test for component significance in MBPLS.

For each component a, the null hypothesis is “there is no real relationship between X and Y at this component”; the test permutes the rows of y, refits a fresh MBPLS with the same number of components, and recomputes the test statistic. The risk is the fraction of permutations whose statistic equals or exceeds the original model’s.

Statistic: per-component absolute correlation between the super X-score and the super Y-score, |t_super(:,a)' u_super(:,a)| / (||t|| * ||u||).

Parameters:

model (MBPLS) – A fitted MBPLS model.
X (dict[str, DataFrame], DataFrame) – The same training data used to fit model.
y (dict[str, DataFrame], DataFrame) – The same training data used to fit model.
n_permutations (int, default=200) – Number of Y-row permutations to evaluate.
seed (int or None, default=None) – Seed for the permutation RNG (None uses non-reproducible randomness).

Returns:

Indexed by component 1..A with columns:

observed : the actual model’s per-component statistic.
risk_pct : Monte-Carlo estimate (in %) of the right-tail probability, 100 * (n_exceed + 1) / (n_permutations + 1). Low values (e.g. < 5%) suggest the component is significant; values near 50% suggest the component is no better than chance.

The + 1 on each side counts the observed statistic among the permutations, which is what keeps the estimate a valid p-value: the uncorrected n_exceed / n_permutations can report exactly 0, and no finite permutation set can license the claim that the true tail probability is zero. The floor is 100 / (n_permutations + 1), so the default 999 permutations cannot resolve below 0.1%. This matches the convention already used by the Van der Voet test in _pls (#513).

Return type:

pd.DataFrame

References

Wiklund, S., Nilsson, D., Eriksson, L., Sjöström, M., Wold, S. & Faber, K. A randomization test for PLS component selection. J. Chemometrics, 21 (2007), 427-439.

MBPCA#

Multi-block PCA / consensus-PCA. Same dict-of-DataFrames API as MBPLS; no Y-block.

class process_improve.multivariate.methods.MBPCA(n_components, *, max_iter=500, tol=None, algorithm='auto', missing_data_settings=None)[source]#

Bases: _HotellingsT2LimitMixin, TransformerMixin, BaseEstimator

Multi-block PCA (hierarchical / consensus PCA).

Generic multi-block PCA following the consensus-PCA / hierarchical PCA formulation of Westerhuis, Kourti & MacGregor (1998). Each X-block is preprocessed independently (mean-centred and unit-variance scaled), then divided by sqrt(K_b) so blocks of unequal width contribute fairly to the consensus super-score.

The hierarchical NIPALS loop alternates: (i) regress each block on the super-score to get block loadings and block scores, (ii) collect block scores into a super-block, (iii) regress the super-block to get a new super-score / super-loading, repeat to convergence. After convergence, deflate every block by the super-score and the corresponding block loading scaled by the super-loading element.

Parameters:

n_components (int) – Number of super-components (consensus latent variables) to extract.
max_iter (int, default=500) – Maximum NIPALS iterations per component in the hierarchical outer loop.
tol (float or None, default=None) – Relative convergence tolerance on the super-score change: the norm of the change between two successive super-score iterations, divided by the norm of the current super-score vector (#504). None uses epsqrt (about 1.49e-8), the same default as PCA / PLS / TPLS. The legacy absolute tolerance np.finfo(float).eps ** (9/10) (about 8.2e-15) sits below the floating-point oscillation floor of a relative criterion, so it would never be reached in practice.
algorithm (str) –
Algorithm to use for fitting the model.
- "auto": dense vectorised hierarchical NIPALS when the data is complete; mask-aware NIPALS (NaN-tolerant) when any block contains missing values.
- "dense": dense vectorised hierarchical NIPALS. Raises if any block contains missing values.
- "nipals": mask-aware hierarchical NIPALS. Always uses the NaN-tolerant inner-loop primitives, even when the data is complete (slower than "dense" but produces equivalent results).
missing_data_settings (dict or None, default=None) – Settings for the iterative "nipals" path. Keys: md_tol (convergence tolerance on the score-vector change between iterations), md_max_iter (maximum NIPALS iterations per component). Defaults to {"md_tol": epsqrt, "md_max_iter": 1000}.
fitting) (Attributes (after)
--------------------------
block_names (list[str]) – Ordered list of X-block names (the keys of the input dict).
block_widths (dict[str, int]) – Number of variables in each X-block.
n_samples (int) – Number of rows fitted.
n_features_in (int) – Total number of X variables summed across blocks.
feature_names_in (np.ndarray) – Concatenated column names, one per feature, in block order.
preproc (dict[str, MCUVScaler]) – Per-block preprocessors used to mean-centre and unit-variance scale each X-block.
super_scores (pd.DataFrame, shape (n_samples, n_components)) – Super-block (consensus) scores T.
super_loadings (pd.DataFrame, shape (n_blocks, n_components)) – Super-block loadings p_super; rows indexed by block name.
super_hotellings_t2 (pd.DataFrame, shape (n_samples, n_components)) – Cumulative Hotelling’s T^2 on the super-scores per component.
block_scores (dict[str, pd.DataFrame]) – Per-block scores t_b, each shape (n_samples, n_components).
block_loadings (dict[str, pd.DataFrame]) – Per-block loadings p_b, each shape (K_b, n_components).
block_spe (dict[str, pd.DataFrame]) – Per-block squared prediction error per sample and component.
block_hotellings_t2 (dict[str, pd.DataFrame]) – Per-block cumulative Hotelling’s T^2 per sample and component.
block_vip (dict[str, pd.Series]) – Per-block variable-importance in projection, indexed by variable name inside each block.
r2_x_per_block_cumulative (pd.DataFrame, shape (n_blocks, n_components)) – Cumulative R^2X per block and component.
r2_x_per_block_per_component (pd.DataFrame, shape (n_blocks, n_components)) – Incremental R^2X per block and component.
r2_x_per_variable (dict[str, pd.DataFrame]) – Cumulative R^2X per variable within each block.
explained_variance (np.ndarray, shape (n_components,)) – Variance of the super-score per component (ddof=1).
scaling_factor_for_super_scores (pd.Series) – sqrt(explained_variance_) per component.
fitting_info (dict) – Per-component iteration count and timing.
has_missing_data (bool) – Whether any X-block had NaN values.
algorithm – The resolved algorithm actually used for the fit. With algorithm="auto", this is "dense" for complete data and "nipals" for NaN-containing data.

Notes

The deflation step is $X_b \leftarrow X_b - t_{\rm super}\, (p_b\,p_s[b]\,\sqrt{K_b})^\top$, derived in Westerhuis et al. 1998. An earlier implementation of this method had this step marked as broken by its author; this implementation re-derives it directly from the paper and is independently validated against the pure-numpy reference oracles in the test suite.

Missing data#

When any block contains NaN entries, the "auto" algorithm routes to a mask-aware NIPALS variant. Each per-block projection in the inner loop is computed as a regression that uses only the observed entries; the masked sum-of-squares is used as the denominator so missing values neither bias the loading direction nor contribute to the score. The mask is preserved across components automatically because deflation propagates NaN through subtraction. This is the standard skip-NaN NIPALS update; see Walczak & Massart (2001) and Arteaga & Ferrer (2002).

The fit refuses to run if any block has a column with all entries missing, or any block has a row with all entries missing for that block; either case leaves the masked denominator at zero. Drop or impute such rows or columns before fitting. Predict-time score estimation for new observations with NaN (Trimmed Score Regression / Projection to the Model Plane) is a separate follow-up.

References

Westerhuis, J. A., Kourti, T. & MacGregor, J. F. Analysis of multiblock and hierarchical PCA and PLS models. J. Chemometrics, 12 (1998), 301-321.

Walczak, B. & Massart, D. L. Dealing with missing data: Part I. Chemom. Intell. Lab. Syst., 58 (2001), 15-27.

Arteaga, F. & Ferrer, A. Dealing with missing data in MSPC: several methods, different interpretations, some examples. J. Chemometrics, 16 (2002), 408-418.

fit(X, y=None)[source]#

Fit the multi-block PCA model.

Parameters:

X (dict[str, pd.DataFrame]) – X-blocks. Keys are block names; values are DataFrames sharing the same row index (and row count). Each block is preprocessed independently.
y (None) – Ignored; accepted so the transformer plugs into a sklearn Pipeline.

Return type:

MBPCA

transform(X)[source]#

Project new data to super-scores using the fitted model.

Parameters:: X (dict[str, DataFrame])
Return type:: DataFrame

diagnose(X)[source]#

Project new data; return super_scores, block_scores, block_spe, hotellings_t2.

The rename (since 1.38.4, #395) matches PCA.diagnose() and PLS.diagnose(); predict() is kept as a deprecation shim.

Parameters:: X (dict[str, DataFrame])
Return type:: Bunch

predict(X)[source]#

Forward to diagnose(); emits a DeprecationWarning.

Deprecated since version 1.38.4: Use MBPCA.diagnose() instead. predict matches the sklearn-convention name but MBPCA isn’t a regressor; the historical return is a diagnostics Bunch. The rename aligns with PCA.diagnose(). Will be removed in 2.0.0.

Parameters:: X (dict[str, DataFrame])
Return type:: Bunch

block_spe_limit(block, conf_level=0.95)[source]#

SPE limit for one X-block (Nomikos & MacGregor chi-square approximation).

Parameters:

block (str)
conf_level (float)

Return type:

super_spe_limit(conf_level=0.95)[source]#

SPE limit for the merged super-block.

Parameters:: conf_level (float)
Return type:: float

spe_contributions(X)[source]#

Per-variable squared residuals for each X-block (SPE contributions).

Reconstruction matches the MBPCA deflation step: X_b = T_super @ (P_b * p_super[b] * sqrt(K_b))^T summed over components. Returns squared residuals on the preprocessed scale; sum across columns equals block_spe_[b].iloc[:, -1] ** 2.

Parameters:: X (dict[str, DataFrame])
Return type:: dict[str, DataFrame]

score_contributions(X, component=1, scaling='none')[source]#

Per-block per-variable contributions to a super-score (MBPCA).

The multi-block analogue of PCA.score_contributions(). A super score is a weighted sum of the (deflated, preprocessed) variables across every block, so it splits exactly into one term per variable:

\[c_{b,ij}^{(a)} = \tilde{x}_{b,ij}^{(a)}\, \frac{P_b[j, a]\, p_\mathrm{super}[b, a]} {(p_b^\top p_b)(p_\mathrm{super}^\top p_\mathrm{super})\sqrt{K_b}}, \qquad \sum_b \sum_j c_{b,ij}^{(a)} = t_{\mathrm{super},ia},\]

where $\tilde{x}^{(a)}$ is the block data deflated through the first $a-1$ components.

See MBPLS.score_contributions() for the parameter and return descriptions; the API is identical.

Parameters:

X (dict[str, DataFrame])
component (int)
scaling (str)

Return type:

dict[str, DataFrame]

group_contributions(X, group, reference=None, component=1)[source]#

Per-block per-variable contributions to a group’s average super score.

The multi-block analogue of PCA.group_contributions(). See that method for the definition; the result is returned one Series per X-block, and the sum over every block equals the group’s average super score (or the difference between the two groups’ average super scores when reference is given).

Parameters:

X (dict[str, DataFrame])
group (Sequence)
reference (Sequence | None)
component (int)

Return type:

dict[str, Series]

super_score_plot(pc_horiz=1, pc_vert=2)[source]#

Scatter plot of MBPCA super-scores for two components.

Parameters:

pc_horiz (int)
pc_vert (int)

Return type:

Figure

super_loadings_bar_plot(component=1)[source]#

Bar plot of MBPCA super-loadings for a single component.

Parameters:: component (int)
Return type:: Figure

display_results(show_cumulative=True)[source]#

Format a short text summary of per-block R²X, iterations and timing.

Parameters:: show_cumulative (bool)
Return type:: str

Analysis#

process_improve.multivariate.methods.rv_coefficient(X, Y)[source]#

Compute the RV coefficient between two data blocks.

The RV coefficient (Robert and Escoufier, 1976) measures how much common structure two matrices, measured on the same observations, share. It is a multivariate generalisation of the squared Pearson correlation: it compares the observation-by-observation configuration matrices $XX^T$ and $YY^T$ rather than individual variables.

Parameters:

X (array-like of shape (n_samples, n_features_x)) – First data block.
Y (array-like of shape (n_samples, n_features_y)) – Second data block. Must have the same number of rows as X; the number of columns may differ.

Returns:

The RV coefficient in the range [0, 1]. A value of 1 means the two blocks describe the same configuration of observations up to a rotation and an overall scaling; 0 means no shared structure. nan is returned if either block has no variance.

Return type:

Notes

Each column is mean-centred internally, since the RV coefficient is defined on centred data. The blocks are not scaled; scale the columns yourself (for example with MCUVScaler) when the variables have different units.

For high-dimensional data (many more variables than observations) the RV coefficient is biased upwards and tends towards 1 even for unrelated blocks. Use rv2_coefficient() in that regime.

References

Robert, P. and Escoufier, Y. (1976). A unifying tool for linear multivariate statistical methods: the RV-coefficient. Journal of the Royal Statistical Society, Series C, 25(3), 257-265.

See also

rv2_coefficient: Modified RV coefficient, unbiased for high-dimensional data.

Examples

>>> rv_coefficient(X, Y)
>>> rv_coefficient(X, X)  # 1.0: a block is perfectly correlated with itself

process_improve.multivariate.methods.rv2_coefficient(X, Y)[source]#

Compute the modified RV coefficient (RV2) between two data blocks.

The modified RV coefficient (Smilde et al., 2009) is a variant of rv_coefficient() that removes the diagonals of the configuration matrices $XX^T$ and $YY^T$ before comparing them. This removes the upward bias that makes the ordinary RV coefficient tend towards 1 for high-dimensional data, so RV2 stays near 0 for genuinely unrelated blocks.

Parameters:

X (array-like of shape (n_samples, n_features_x)) – First data block.
Y (array-like of shape (n_samples, n_features_y)) – Second data block. Must have the same number of rows as X; the number of columns may differ.

Returns:

The modified RV coefficient, in the range [-1, 1]. A value of 1 means the two blocks describe the same configuration of observations; values near 0 mean no shared structure, and small negative values can occur. nan is returned if either block has no variance.

Return type: