Multivariate Analysis#

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

Estimate scores and diagnostics for rows that may contain missing values.

This is the “batch so far” primitive of online batch monitoring: the future part of an unfolded batch row is missing by construction, and the score estimate at each time sample is this projection (see Garcia-Munoz, Kourti and MacGregor, 2004).

Parameters:

X (array-like of shape (n_samples, n_features)) – New data in the model’s (centred and scaled) space; NaN marks a missing entry. A row with nothing observed sits at the model centre, as fit() places one: scores, T2 and SPE of 0, no observed feature, and a condition_number of inf.
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, n_samples x n_components), hotellings_t2 (Series; total over all components, computed with the training score variances), spe (Series; the square root of the residual sum of squares over the observed columns only), condition_number (Series; the conditioning diagnostic of each row’s estimator, 1.0 when nothing is missing) and n_observed (Series; observed features per row). SPE and T2 for a partially-observed row must be compared against limits built from the same missingness pattern, not against the full-observation limits; see process_improve.batch.BatchMonitor.

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. 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, ekf_max_iter=50, ekf_tol=1e-06, random_state=None, return_consensus=False)[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, so PRESS does not shrink monotonically with the component count the way it does when a held-out row is scored through its own values.

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. Ignored by "sacv" and "gcv", which hold nothing out.
cv_scheme ({"ekf", "ckf", "ek", "sacv", "gcv"}, default "ekf") –
How a held-out value is produced. Every one of these keeps the prediction independent of the value being predicted; they differ in what they hold out and what they cost. "row_wise" was removed in 2.0 and raises a ValueError: it scored a held-out row through its own values, so it measured compression rather than prediction and could not pick a component count.
- "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. A component count that leaves some group with fewer retained columns than components cannot be scored there, so its PRESS and Q2 are NaN and it is not recommended.
- "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.
"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, ekf_max_iter and ekf_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.
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 the per-fold spread that "ekf" and "ckf" report). "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.

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().
ekf_max_iter (int and float, default 50 and 1e-6) – EM iteration cap and convergence tolerance for the ekf imputation step. Named for the loop they bound: that EM loop is not the per-component NIPALS loop, which PCA(tol=..., max_iter=...) governs (#588). Their old names n_iter and tol were removed in 2.0.
ekf_tol (int and float, default 50 and 1e-6) – EM iteration cap and convergence tolerance for the ekf imputation step. Named for the loop they bound: that EM loop is not the per-component NIPALS loop, which PCA(tol=..., max_iter=...) governs (#588). Their old names n_iter and tol were removed in 2.0.
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). The parallel analysis runs with its default scale=True whatever scale_inside_folds is.

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).
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.
cv_scores - alias of per_fold_press.
cv_scheme - the scheme used, echoed back (one of "ekf", "ckf", "ek", "sacv" or "gcv").
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.

detect_outliers(conf_level=0.95, *, judged_rows=None)[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).
judged_rows ({"fitted", "new"} or None, optional) – Which rows the T² limit judges (see hotellings_t2_limit()). This method judges the rows the model was fitted to, so "fitted" gives their limit. None, the default, gives the "new" limit until 3.0 and the "fitted" one from 3.0, with a FutureWarning.

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.
hotellings_t2_limit_note - only when judged_rows is None and the limit changes in 3.0: the text of the warning.

Return type:

list of dict

Examples

>>> pca = PCA(n_components=3).fit(X_scaled)
>>> outliers = pca.detect_outliers(conf_level=0.95, judged_rows="fitted")
>>> 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=500, tol=np.float64(1.4901161193847656e-08), copy=None, 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=500) – Maximum number of iterations per component for the NIPALS algorithm. The same default as every other iterative estimator in the package; it was 1000 before 1.98.0 (#588).
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 (None) – Deprecated since 2.4.0 and removed in 3.0; it has no effect. PLS never modifies the X or Y it is given, so there is nothing to copy.
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) –
How to fit when the data has missing cells. Keys:
- md_method, one of:
  - "nipals" (the default): NIPALS skips each missing cell in its sums, so an incomplete row simply contributes less. Cheapest, and no imputation happens.
  - "tsr": trimmed score regression. The missing cells are filled by EM on a principal-component model of X and Y together, refitted every round from the completed data, and PLS is then fitted once to the result. Measured against the model the complete data would have given: on synthetic data with clear low-rank structure it lands 2.0 to 2.9 times closer than "nipals" (5-45% of cells missing at random, closer on 17-19 trials in 20); on the LDPE process data 1.1 to 1.5 times closer (9-12 in 12). The gain grows with how much of the data the components capture, because a missing cell can only be rebuilt from what they describe. Roughly five times the cost of a "nipals" fit.
  - "pmp": the same loop, estimating each row by least-squares projection onto the components instead of by regression. Close to "tsr" throughout, and slower to settle as more goes missing.
  Any other value is refused, including "scp": that is a projection-time method (see project()), and at fit time it is simply what NIPALS already does.
- md_tol and md_max_iter were removed in 2.0 and raise a ValueError; set tol and max_iter instead. Those two bound the NIPALS loop and, for "tsr" / "pmp", the imputation as well: the tolerance on the largest change in any imputed cell (in standard deviations) and the cap on imputation rounds.
With "tsr" / "pmp", R², SPE and the other diagnostics are still computed on the observed cells only: an imputed cell is fitted by construction, so counting it would flatter the model.
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. A row with nothing observed has no residual, and an SPE of NaN, as in diagnose.
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. A fit with sample_weight measures this and every other R², and rmse_, in its weights, so a row of weight 0 is left out of them as it is left out of the fit.
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()).

A missing cell (NaN) is skipped, as the NIPALS fit skips it: a row’s score is regressed on the observed part of $w_{\text{TP}}$, and a feature’s loading on the scores of the rows where it is observed. A row with nothing observed scores 0 and adds nothing, and a feature with nothing observed has a loading of 0.

Parameters:

model (PLS) – A fitted PLS model (must expose beta_coefficients_).
X (array-like of shape (n_samples, n_features)) – The data the model’s predict takes, in the units it was fitted on. A model that centres and scales X itself (PLS(scale=True), the default) does so here first; PLS(scale=False) does not, and is given X scaled as its training data was.
response (str or int or None, default=None) – Which response (Y column) to project onto. None is allowed only for a single-response model; otherwise pass the response label (or its integer position).

Returns:

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.

A feature the TP component reproduces exactly has no residual, and its SR is inf, so a descending sort ranks it first. The residual counts as none when its sum of squares is at most machine epsilon (about 2.2e-16) times the explained one, which is a residual norm of at most about 1.5e-8 times the feature’s: rounding alone leaves about 1e-32 of its sum of squares. Every finite SR is therefore below 1 / eps, about 4.5e15. A feature with no sum of squares at all, such as a constant column after MCUVScaler, has neither an explained nor a residual part (0 / 0), and its SR is 0.

A missing cell (NaN) is skipped, as in target_projection(), and both sums of squares of a feature run over its observed cells, where they add up to the feature’s sum of squares. A feature with nothing observed has none, and its SR is 0. The F critical value of a feature counts its observed rows too.

Parameters:

model (PLS) – A fitted PLS model (must expose beta_coefficients_).
X (array-like of shape (n_samples, n_features)) – The data the model’s predict takes, in the units it was fitted on. A model that centres and scales X itself (PLS(scale=True), the default) does so here first; PLS(scale=False) does not, and is given X scaled as its training data was.
response (str or int or None, default=None) – Which response to compute SR for. None returns a feature-by-response DataFrame when the model has several responses, or a Series for a single-response model.
conf_level (float, default=0.95) – Confidence level for the advisory F-based critical value, attached to the result’s .attrs["f_critical"]: a Series indexed by feature, with $N - 2$ and $N - 3$ degrees of freedom, $N$ being the number of rows the feature is observed in. So ratio > ratio.attrs["f_critical"] compares each feature with its own.

Returns:

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

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.

optimal_inversion(rules, *, alternatives=0, leave_out=1)#

Find the inputs a PLS or TPLS model vouches for that meet the rules, at the least cost or closest to the targets.

Parameters:

model (PLS, TPLS or AffineModel) – A fitted model, or its maps as affine_model() read them at rules.conf_level (pass the maps to solve many problems on one model, or to search with maps stored by AffineModel.to_dict()).
rules (InversionRules) – The objective (costs or targets), the windows, bounds, fixed inputs, linear rules, counting rules and limits.
alternatives (int) – How many further answers to find, each with a different set of inputs: the best answer that leaves out at least leave_out of the counted inputs used by the answer and by every alternative before it (a no-good cut on the master problem). Needs counted inputs: a group with a max_count or a min_amount. Give the min_amount, so that a trace of an input does not count as using it.
leave_out (int) – How many of each earlier answer’s counted inputs every alternative must leave out.

Returns:

result – With keys:

feasiblebool: Whether any inputs meet every rule. When False, message says which rule is in the way.
x_newpd.Series: The inputs, in their original units.
y_hatpd.Series: The predicted qualities.
scores, hotellings_t2, t2_limit, spe, spe_limit: How far the answer sits from the data the model was built from, against its limits.
objectivefloat: The cost, or the summed distance to the targets in training standard deviations.
useddict of str to list of str: The inputs used from each group.
constraintspd.DataFrame: Every rule, window and limit, and each input at a bound: its value, bounds, the side it binds on, and its shadow_price, the change in the optimal objective per unit increase of that bound (per_percent: per 1% of it). A rule that does not bind has a shadow price of 0. A binding constraint implied by others (a lot at 100% of its material, held both by its bound and by the material’s total) has NaN: its price is not its own.
kkt_residualfloat: How nearly the optimality conditions hold at the answer; near 0 at an optimum.
proven_optimalbool: The search proved no better answer exists, and the answer is optimal for its inputs: the optimality conditions hold, or a linear program over the limits’ tangent planes at the answer does no better. With no answer: the search proved none exists; False when it stopped before it could.
n_cutsint: The tangent planes the counting rules needed.
messagestr: Why there is no answer, or why it is not proven optimal; empty otherwise.
alternativeslist of sklearn.utils.Bunch: Only when alternatives is asked for: each alternative as a result with the keys above, best first. Its proven_optimal means no input set the earlier answers allow does better. The list ends early, with a result whose feasible is False and whose message says why, when no further input set meets the rules; it is empty when there is no answer at all.

Return type:

Raises:

ValueError – The rules name what the model does not have, or cannot be read; the maps given were read at another confidence level; or alternatives are asked for without counted inputs, or with a count that is not a whole number.

Warns:

FutureWarning – An input that is searched has no lower or no upper bound in rules.bounds (a TPLS formula’s fractions have 0 and 1 unless given). Deprecated since 2.4.0, and an error from 3.0.

See also

inversion_sensitivity: check the shadow prices by solving again, and trace a limit, counts and prices.

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.

M. A. Duran and I. E. Grossmann, “An outer-approximation algorithm for a class of mixed-integer nonlinear programs”, Mathematical Programming, 36 (1986): 307-339.

Examples

>>> rules = InversionRules(
...     costs={"Acetic": 1.0, "H2S": 3.0, "Lactic": 8.0}, windows={"Taste": (20.0, 30.0)},
...     bounds={"Acetic": (4.0, 7.0), "H2S": (3.0, 10.0), "Lactic": (0.8, 2.1)},
... )
>>> result = optimal_inversion(pls, rules)
>>> result.x_new, result.objective, result.constraints

inversion_sensitivity(rules, *, front=None, prices=None, fraction=0.01)#

How the optimal inversion moves with its rules: shadow prices checked by solving again, a traced limit, the value of one input more or fewer per counted group, and how far costs can move before the answer changes.

Parameters:

model (PLS, TPLS or AffineModel) – A fitted model, or its maps, as for optimal_inversion().
rules (InversionRules) – The rules, as for optimal_inversion().
front ((str, list of float), optional) – A rule to trace, with the values its bound takes: "T²", an SPE limit ("SPE X"), a quality with a window, or a linear rule’s name. A two-sided rule moves the side that binds. The optimal objective at each.
prices (list of str, optional) – Inputs whose cost to range (costs objective only): how far, as a fraction of itself, each can fall or rise before the answer changes shape (another input used, another rule binding). Each costs about a dozen solves.
fraction (float, default 0.01) – How far each binding rule is moved to check its shadow price, as a fraction of its bound.

Returns:

result – With keys inversion (the optimal_inversion() result), checked (each binding rule: moved by, forecast change from its shadow price, actual change on solving again), front (value, feasible, objective), counts (group, max_count, feasible, objective, change) and prices (input, cost, used, falls_by, rises_by).

Return type:

Examples

>>> report = inversion_sensitivity(pls, rules, front=("T²", [4.0, 6.0, 8.0]), prices=["Lactic"])
>>> report.checked   # shadow prices against solving again
>>> report.front     # the cost as the T² limit moves

y_scores_: ndarray | DataFrame#

y_weights_: ndarray | DataFrame#

y_loadings_: ndarray | DataFrame#

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

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

Raises:

NotEnoughVarianceError – If every column of X is constant or has no observed value over the rows the fit uses, as across the family (#1144); or if every column of Y is, as the fit centres it (with scale=False, if Y is 0 throughout).

Warns:

UserWarning – If a column of X or Y has no observed value. It is left out of the model, as across the family (#1098): its weights and loadings are 0, a value a new row has in an X column of this kind adds nothing to that row’s scores or SPE, and a response of this kind has NaN predictions.

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. A row with a missing cell (NaN) is scored from the cells it has, as diagnose() scores it.
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): those of diagnose().

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)) – New rows. A row with a missing cell (NaN) is scored from the cells it has, by trimmed score regression, the estimator project() uses by default, and predicted from those scores; until 2.4.0 its scores, SPE and predictions were NaN and its T² 0 (#1144). A row with nothing observed is NaN throughout. So is a row with too few observed cells for its scores, typically fewer than n_components, with a SpecificationWarning. project() gives the other estimators, a ridge, and each row’s conditioning.

Returns:

result – With keys:

scores (pd.DataFrame, (n_samples, n_components)).
cumulative_hotellings_t2 (pd.DataFrame, (n_samples, n_components)): Hotelling’s T², cumulative, so column a holds the T² over the first a components and the last column is the full model’s, as in hotellings_t2_. Every diagnose of the family returns it in this shape (#983).
spe (pd.Series, (n_samples,)): SPE of the full model.
y_hat (pd.DataFrame, (n_samples, n_targets)): the predictions, on the original Y scale.
residuals (pd.DataFrame, (n_samples, n_features)): what is left of each row after the components, in the scaled space the model was fitted in, NaN at a missing cell; spe is the length of each row over its other cells.
hotellings_t2 (pd.Series, (n_samples,)): deprecated since 2.4.0, to be removed in 3.0; the T² of the full model, which is cumulative_hotellings_t2.iloc[:, -1]. Reading it warns.

Return type:

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.cumulative_hotellings_t2.iloc[:, -1]   # T² for each new observation

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

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

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. transform(), diagnose() and predict() score a row with a missing cell by this method’s default, through the same code; this method adds the choice of estimator, the ridge, and each row’s condition_number, and raises where they give NaN for a row the estimator cannot be inverted for (#1144).

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. A row with nothing observed sits at the model centre, as fit() places one: scores, T2 and SPE of 0, the training mean of Y as its prediction, no observed feature, and a condition_number of inf.
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. A DataFrame must have one row: one target is inverted at a time. A response the fit left out, having no observed value, cannot be targeted: its value is not used, and may be missing.
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; NaN for a column the fit left out, having no observed value.
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. NaN for a response the fit left out, as predict() gives it.
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

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 four 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, numpy.random.Generator or None, optional) – Seed forwarded to KFold / RepeatedKFold for reproducible shuffling; a Generator gives them its next draw of an int. 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. A fold model that scales (scale=True, the default) still autoscales its own rows, so pass scale=False as well to keep a scaling of your own. Each held-out row is predicted as its fold model predicts it, whichever scalings apply.

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 compare_cv_criteria(X, Y, *, max_components=None, cv=7, random_state=None, n_permutations=999, n_cv_permutations=None, alpha=0.05, angle_threshold=30.0, conf_level=0.95, **pls_kwargs)[source]#

Compare predictive and latent-structure validation criteria for a PLS model.

Runs one K-fold cross-validation (one fit of this class per fold at max_components, truncated for smaller counts) and reports, for every number of components, criteria that each answer a different question, together with the number of components each one recommends. Disagreement between them is the point: a component can predict Y without being stable, or be stable and real without improving the prediction.

Missing values (NaN) are handled as the NIPALS fit handles them. Every held-out row is scored the way NIPALS scores a training row, from its observed cells, and a missing Y cell is left out of PRESS, of the slope ratio and of the score correlation. The PRESS identity and s_a = 1 in sample stay exact.

Parameters:

X (array-like of shape (n_samples, n_features)) – Predictor block, on its raw scale, NaN where missing. Each training fold is autoscaled with its own MCUVScaler, as in PLS.select_n_components().
Y (array-like of shape (n_samples,) or (n_samples, n_targets)) – Response block, NaN where missing. A row with no observed response still helps fit X.
max_components (int, optional) – Largest number of components to evaluate. Capped at one fewer than the smallest training fold, and at the number of features.
cv (int or sklearn splitter, default=7) – Number of shuffled K-fold segments, or a splitter that holds out every row exactly once (KFold, LeaveOneOut, GroupKFold, …). With an integer random_state the folds are those of PLS.select_n_components(..., cv=cv, n_repeats=1, random_state=random_state).
random_state (int, numpy.random.Generator or None, default=None) – Seed for the fold assignment and for the permutation tests.
n_permutations (int, default=999) – Permutations for the van der Voet and the covariance tests (cheap: no refits).
n_cv_permutations (int or None, default=None) – Permutations of the rows of Y used to calibrate the held-out score correlation of each component. Each one repeats the cross-validation with kernel PLS refits from the folds’ X'X and X'Y (fast for tens to hundreds of features). For PRM, each fold is refitted by PRM itself, so that the correlation and its null come from the same estimator. That is one PRM fit per fold per permutation, about a hundred times slower per permutation. None uses 199 permutations, or 49 for PRM: there the smallest attainable p-value is 0.02 rather than 0.005, still enough to decide at alpha=0.05, and the test runs in about a quarter of the time 199 would take. With 0 the normal approximation $r \sim N(0, 1/N)$ is used instead; it is quick but too permissive, because fold models share rows and a strong low-rank X widens the null distribution of a pooled cross-validated correlation.
alpha (float, default=0.05) – Significance level for the tests and for the score-correlation threshold. PLS.select_n_components() uses 0.01 for van der Voet.
angle_threshold (float, default=30.0) – Jackknife-scaled subspace angle, in degrees, above which the leading weights are judged unstable.
conf_level (float, default=0.95) – Confidence level of the SPE and Hotelling’s $T^2$ limits whose out-of-sample alarm rates are checked.
**pls_kwargs – Passed to every fit of this class (for example max_iter, tol). scale=False is rejected: every fold is autoscaled.

Returns:

result – With these fields:

table (pandas.DataFrame, indexed by n_components), with columns:
- r2y, q2y, q2y_se: in-sample $R^2_Y$ (pooled in scaled units), cross-validated $Q^2_Y$ (every target weighted equally, which for one target is the ordinary $Q^2$), and its standard error across folds.
- vdv_p: van der Voet p-value against the minimum-PRESS model.
- cv_anova_p: CV-ANOVA p-value for the whole model, one response only (NaN otherwise). It is a monotone function of $Q^2$.
- r_train, r_cv, r_cv_threshold, r_cv_p: correlation of $t_a$ and $u_a$ on the training rows; the same correlation pooled over the held-out rows (about the training-fold centre, so it is unaffected by sign flips between folds); the $1-\alpha$ quantile of its null distribution under permuted Y; and its one-sided permutation p-value.
- slope_ratio: $s_a$, the held-out slope of the deflated response on $t_a$ relative to the training slope. PRESS falls exactly when $s_a > 1/2$; see the module notes.
- cov_perm_p: sequential permutation p-value for the covariance of component a.
- angle_component_deg, angle_subspace_deg: angle between each fold’s weight vector $w_a$ and the full-data one, and the largest principal angle between the spans of the first a weights, each scaled by the delete-d jackknife ($\tan\theta \to \sqrt{G-1}\,\text{rms}(\tan\theta_k)$) so that it estimates how far the full-data direction may lie from the population one. Raw fold angles shrink as the number of folds grows (fold models share most of their rows); the scaled ones do not. A large component angle with a small subspace angle means the components swap or mix inside a stable space: loadings are then not interpretable one by one, but SPE, $T^2$ and inversion, which depend on the span, are.
- d_ratio_median, d_ratio_min: Procrustes D ratios $c_{ka}^\top c_a / c_a^\top c_a$ across folds; a negative minimum means some fold reversed the component’s inner relation.
- pv_spe_alarm_rate, pv_t2_alarm_rate: fraction of the Procrustes pseudo-validation rows above the full-data model’s SPE and $T^2$ limits at conf_level; nominally 1 - conf_level.
- pv_spe_limit_ratio: the SPE limit fitted to the pseudo-validation (held-out) SPE, divided by the full-data limit. The full-data limit is fitted to training residuals, which are smaller than those of new rows, so the ratio is above 1, and more so the fewer rows per variable.
recommendations (pandas.DataFrame): one row per selection rule (q2_max, q2_1se, van_der_voet, score_correlation, covariance_permutation, subspace_stability) with the recommended n_components, the rule and the question it answers. The structural rules count leading components that pass, stopping at the first failure, and can return 0. For subspace_stability, a pair of components that swap inside a stable span (the span of the first fails, the span of both passes) counts as two.
press (pandas.Series): cross-validated PRESS in original Y units.
press_baseline (float): PRESS of predicting each held-out row by its training fold’s mean.
d_ratios (pandas.DataFrame): Procrustes D ratio per fold and component.
spe_limits (pandas.DataFrame, indexed by n_components): the SPE limit at conf_level fitted to the training residuals (full_data, what global_model reports) and to the pseudo-validation SPE (pseudo_validation). Use the second to monitor new rows. It is a limit for complete rows: a held-out row with missing cells has its squared SPE scaled by K / n_observed before the fit, and the SPE of a new row with gaps should be scaled the same way before the comparison.
cv_splits (list of (train, test) index arrays): the folds used.
global_model (PLS): the model fitted to all rows with max_components components.
alpha, conf_level, angle_threshold, alarm_rate_upper (float): the settings used, and the binomial upper bound for the alarm rates.

Return type:

Raises:

ValueError – For a row whose X is entirely missing, a column that is entirely missing, a training fold with fewer than two observed values in a column, a splitter that does not hold out every row exactly once, or an out-of-range setting.

See also

select_n_components: Choose the number of components from prediction alone.
pseudo_validation_set: Build the Procrustes pseudo-validation set explicitly.
process_improve.multivariate.plots.cv_criteria_plot: Plot the table.

Examples

>>> result = PLS.compare_cv_criteria(X, y, max_components=6, random_state=0)
>>> result.recommendations["n_components"]
>>> result.table[["q2y", "r_cv", "slope_ratio", "angle_subspace_deg"]]

classmethod pseudo_validation_set(X, Y, *, n_components, cv=7, scope='global', random_state=None, **pls_kwargs)[source]#

Build a Procrustes pseudo-validation set for a PLS model.

Procrustes cross-validation (Kucheryavskiy, Rodionova and Pomerantsev, 2023) turns the variation between cross-validation fold models into a data set, X_pv, of the same size as X, which the model fitted to all rows can be applied to like an independent test set. For every held-out row, the full-data model returns:

scores equal to the row’s scores in its fold model, multiplied per component by the D ratio $c_{ka}^\top c_a / c_a^\top c_a$;
an SPE equal to the row’s SPE in its fold model.

With one response and scope="global", the full-data model’s predictions for X_pv also equal the fold models’ predictions for the held-out rows. Y_pv is Y unchanged.

With missing values in X, the held-out rows are scored from their observed cells, X_pv has the same missing cells as X, and the properties above are exact for the complete rows.

Parameters:

X (array-like of shape (n_samples, n_features)) – Predictor block, raw scale, NaN where missing.
Y (array-like of shape (n_samples,) or (n_samples, n_targets)) – Response block, NaN where missing.
n_components (int) – Number of components of the model being validated.
cv (int or sklearn splitter, default=7) – Number of shuffled K-fold segments, or a splitter that holds out every row exactly once.
scope ({"global", "local"}, default="global") – "global" (the published default) fits the fold models on the rows as the model fitted to all rows scales them (autoscaled, for PLS), without re-centring them; the prediction property above is then exact. "local" autoscales each training fold separately, as PLS.select_n_components() does; scores and SPE are still reproduced exactly, predictions approximately. A fold model that centres its rows itself, as PRM always does, predicts approximately in either scope.
random_state (int, numpy.random.Generator or None, default=None) – Seed for the fold assignment and the random directions of the residual part.
**pls_kwargs – Passed to every fit of this class. scale=False is rejected.

Returns:

result – With these fields:

X_pv (pandas.DataFrame): pseudo-validation predictors, on the raw scale of X.
Y_pv (pandas.DataFrame): Y, unchanged.
scores (pandas.DataFrame): the scores the full-data model gives X_pv.
local_spe (pandas.Series): SPE of each held-out row in its fold model; the full-data model gives X_pv the same values. A limit fitted to these, spe_calculation(pv.local_spe, conf_level=0.95), is calibrated for new rows; the full-data model’s own limit is too tight for them.
d_ratios (pandas.DataFrame): D ratio per fold and component.
cv_splits (list of (train, test) index arrays): the folds used.
global_model (PLS): the model fitted to all rows; apply it to X_pv with diagnose.
scope (str): the scope used.

Return type:

Raises:

NotEnoughVarianceError – If n_components is more than the data support: a component fitted after there is nothing left in Y that X can explain has a Y-loading of 0, and so no D ratio.

See also

compare_cv_criteria: Uses the same construction to report alarm rates per component.

Examples

>>> pv = PLS.pseudo_validation_set(X, y, n_components=2, random_state=0)
>>> diagnostics = pv.global_model.diagnose(pv.X_pv)
>>> (diagnostics.spe > pv.global_model.spe_limit()).mean()   # out-of-sample SPE alarm rate

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, numpy.random.Generator or None, 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. A Generator gives its next draw of an int, which is then used as an int is.
**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, *, judged_rows=None)[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].
judged_rows ({"fitted", "new"} or None, optional) – Which rows the T² limit judges, as in PCA.detect_outliers(): "fitted" for the rows the model was fitted to, which this method judges. None gives the "new" limit until 3.0 and the "fitted" one from 3.0, with a FutureWarning.

Returns:

outliers – Sorted from most severe to least. Each dict contains observation, outlier_types, spe, hotellings_t2, spe_limit, hotellings_t2_limit, severity, and, only when judged_rows is None and the limit changes in 3.0, hotellings_t2_limit_note, the text of the warning.

Return type:

list of dict

Examples

>>> pls = PLS(n_components=3).fit(X_scaled, Y_scaled)
>>> outliers = pls.detect_outliers(conf_level=0.95, judged_rows="fitted")
>>> 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, numpy.random.Generator 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. A jackknife or K-fold result has one; a bootstrap result does not, and is refused.

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=500, tol=np.float64(1.4901161193847656e-08), copy=None, 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.
missing_data_settings (dict | None) – Passed through to PLS unchanged.
warn_on_uncentred (bool) – Passed through to PLS unchanged.
copy (None) – Deprecated since 2.4.0 and removed in 3.0; it has no effect. PLSDA never modifies the X or y it is given, so there is nothing to copy.

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`. Its cross-validation

(:meth:`select_n_components`, :meth:`compare_cv_criteria`, :meth:`pseudo_validation_set`,

:meth:`nested_cv` and :meth:`cross_validate`) takes the class labels and runs PLS's on

their class-indicator matrix. An integer ``cv`` shuffles each class into the folds on its

own, as :class:`~sklearn.model_selection.StratifiedKFold` does, with the seed PLS would

shuffle the rows with, so no training fold lacks a class.

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, if a row has no label, 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, or its length does not match X.

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:

classmethod select_n_components(X, Y, *, cv=5, **pls_kwargs)[source]#

Run PLS.select_n_components() on the class-indicator matrix of the labels Y.

Parameters:

X (ndarray | DataFrame)
Y (ndarray | DataFrame | Series)
cv (int | BaseCrossValidator)

Return type:

classmethod compare_cv_criteria(X, Y, *, cv=7, **pls_kwargs)[source]#

Run PLS.compare_cv_criteria() on the class-indicator matrix of the labels Y.

global_model is the PLS fitted to that matrix.

Parameters:

X (ndarray | DataFrame)
Y (ndarray | DataFrame | Series)
cv (int | BaseCrossValidator)

Return type:

classmethod pseudo_validation_set(X, Y, *, cv=7, **pls_kwargs)[source]#

Run PLS.pseudo_validation_set() on the class-indicator matrix of the labels Y.

Y_pv is that matrix, and global_model the PLS fitted to it.

Parameters:

X (ndarray | DataFrame)
Y (ndarray | DataFrame | Series)
cv (int | BaseCrossValidator)

Return type:

classmethod nested_cv(X, Y, *, outer_cv=5, inner_cv=5, **pls_kwargs)[source]#

Run PLS.nested_cv() on the class-indicator matrix of the labels Y.

Parameters:

X (ndarray | DataFrame)
Y (ndarray | DataFrame | Series)
outer_cv (int | BaseCrossValidator)
inner_cv (int)

Return type:

cross_validate(X, Y, **settings)[source]#

Run PLS.cross_validate() for a PLS with this model’s settings and its labels Y.

That PLS is fitted to the class-indicator matrix of Y, and its coefficients are the ones resampled.

Parameters:

X (ndarray | DataFrame)
Y (ndarray | DataFrame | Series)

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 (Serneels et al., 2005): PLS with a bounded influence per observation, so that a handful of outlying rows move the fit a little instead of a lot. X is centred at its L1 median and y at its median, 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 the regression on the scores changes by less than 1%. With scale=False it is the paper’s algorithm, and reproduces the R package chemometrics. 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, coefficient_tol=0.01, weight_tol=None, scale=True, max_iter=1000, tol=np.float64(1.4901161193847656e-08), copy=None, 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^t$, 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^t$.

The steps are those of Serneels et al. (2005):

X is centred at its L1 median and y at its median, once. Neither moves when an outlier does, so the outliers do not choose the coordinate system.
The starting weights come from the data alone: each centred row’s distance from the centre, over the median distance, and each centred y, over the median absolute y.
A PLS model is fitted with the weights, and new weights come from it: the residual, centred at its median and over its median absolute value; and the distance of the row’s scores from their L1 median, over the median distance, each component scaled to unit length in the weighted fit (as SIMPLS’s scores are), so the first components do not decide it alone.
Step 3 repeats until the norm of the regression of y on those unit-length scores changes by less than coefficient_tol, relative.
The intercept is the median of the residuals, $y - Xb$.

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.
max_weight_iter (int, optional) – Maximum reweighting iterations, default 100. Each one is a full PLS fit.
coefficient_tol (float, optional) – The relative change in the norm of the regression of y on the unit-length scores below which the reweighting stops, default 0.01 as in the paper.
weight_tol (None) – Deprecated since 2.2.2 and removed in 3.0; it has no effect. PRM stops on coefficient_tol, as the paper does, not on the change in the weights.
scale (bool, optional) – True (the default) also divides each column of X and of Y, after the robust centring, by 1.4826 times its median absolute deviation (by 1 where that is 0), computed once: columns in different units then count equally, as PLS’s scale makes them. False is the paper’s algorithm, which centres only.
max_iter (int) – As for PLS.
tol (float) – As for PLS.
missing_data_settings (dict | None) – As for PLS.
copy (None) – Deprecated since 2.4.0 and removed in 3.0; it has no effect. PRM never modifies the X or Y it is given, so there is nothing to copy.
warn_on_uncentred (bool, optional) – Accepted for PLS’s signature; it has no effect, since PRM always centres.

robust_weights_#

The weight the returned model gives each row, in [0, 1]: its residual weight times its leverage weight, from the final fit, as the paper reports them. Small entries are this model’s statement about which rows it declined to be led by; see outlier_summary(). A weight is 0 only when the product underflows, for a row so far out in both its residual and its leverage (1e100 standardised units, say) that it takes no part in the fit.

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

n_weight_iter_#

Reweighting iterations actually run.

Type:: int

weights_converged_#

Whether the loop met coefficient_tol within max_weight_iter.

Type:: bool

weight_shift_#

The largest change in any row’s weight on the final iteration.

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. The paper is written for one; with several, a row’s residual is the norm of its residuals, each standardised alone.
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. A row with a zero weight is left out of every median and centre, as if dropped.

Returns:

self, fitted.

Return type:

PRM

Raises:

ValueError – If cutoff, max_weight_iter or coefficient_tol is out of range, if the data contain missing values (see Notes), if X or Y has no columns, or if Y has an infinite value or a row count other than X’s.

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:

CA#

Correspondence analysis of a two-way contingency table: “PCA for a table of counts”. Each row’s profile is compared with the average profile by the chi-squared distance, so a rare category that departs strongly from independence counts for as much as a common one that departs a little. On the map, rows with similar profiles sit together, and a row lying out toward a column co-occurs with it more often than independence would predict.

class process_improve.multivariate.methods.CA(n_components=2)[source]#

Correspondence analysis of a two-way contingency table.

With $P = N / n$ the table of relative frequencies, $r$ and $c$ its row and column masses (the marginal proportions), the matrix of standardised residuals from independence,

\[S = D_r^{-1/2} \, (P - r c^T) \, D_c^{-1/2},\]

is decomposed as $S = U \Sigma V^T$. The squared singular values are the principal inertias; their sum is the total inertia, $\chi^2 / n$. The principal coordinates are $F = D_r^{-1/2} U \Sigma$ for the rows and $G = D_c^{-1/2} V \Sigma$ for the columns, and Euclidean distance between two rows on the full map is exactly the chi-squared distance between their profiles.

Parameters:: n_components (int, optional) – Number of axes to keep, default 2. A table with $I$ rows and $J$ columns has at most $\min(I, J) - 1$ axes; a larger request is clamped, with a SpecificationWarning.

n_components_#

Axes actually kept, after clamping.

Type:: int

eigenvalues_#

Principal inertia of each kept axis.

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

total_inertia_#

Total inertia of the table, over every axis, kept or not: $\chi^2 / n$.

Type:: float

explained_inertia_#

Share of the total inertia on each kept axis.

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

row_masses_, column_masses_

Marginal proportions of the table.

Type:: pd.Series

row_coordinates_, column_coordinates_

Principal coordinates, one row per category, one column per axis.

Type:: pd.DataFrame

row_contributions_, column_contributions_

Share of each axis’s inertia contributed by each category; each column sums to 1. A category that contributes little to an axis is not what defines it, however far out it plots.

Type:: pd.DataFrame

row_cos2_, column_cos2_

Squared cosine of each category with each axis: the share of that category’s own inertia the axis represents. A category that plots near the origin with a low cos2 is not “average”, it is poorly represented on these axes.

Type:: pd.DataFrame

References

Greenacre, “Correspondence Analysis in Practice”, 3rd edition, 2017.

Examples

>>> ca = CA(n_components=2).fit(defects_by_line)
>>> ca.explained_inertia_
>>> ca.map_plot()

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

Fit the correspondence analysis to a contingency table.

Parameters:

X (array-like of shape (n_rows, n_columns)) – The table of non-negative counts. Row and column labels are kept when it is a DataFrame.
y (ignored) – Accepted for Pipeline compatibility.

Returns:

self, fitted.

Return type:

Raises:

ValueError – If n_components is not positive, or the table is not a valid contingency table (see _as_table()).

transform(X)[source]#

Place rows on the fitted map, from their profiles.

A supplementary row is positioned by the transition formula: its profile (the row divided by its total) times the column standard coordinates. It takes no part in defining the axes, which is what lets a new batch, or a row held out of the fit on purpose, be shown against the existing structure. Transforming the fitted table itself returns row_coordinates_.

Parameters:: X (array-like of shape (n_rows, n_columns)) – Counts with the same columns, in the same order, as the fitted table.
Returns:: Principal coordinates, one row per input row.
Return type:: pd.DataFrame
Raises:: ValueError – If the column count differs from the fitted table’s, or the rows are not valid counts.

transform_columns(X)[source]#

Place supplementary columns on the fitted map.

The column counterpart of transform(): X has one row per fitted row category and one column per supplementary column, such as an outcome label (pass / fail) cross-tabulated against the fitted rows.

Parameters:: X (array-like of shape (n_rows, n_new_columns)) – Counts with the same rows, in the same order, as the fitted table.
Returns:: Principal coordinates, one row per supplementary column.
Return type:: pd.DataFrame
Raises:: ValueError – If the row count differs from the fitted table’s, or the columns are not valid counts.

map_plot(axis_x=1, axis_y=2, settings=None)[source]#

Draw the symmetric map: rows and columns together, both in principal coordinates.

Distances within the rows, and within the columns, are chi-squared distances. A distance between a row point and a column point has no such meaning on the symmetric map; what does is the direction: a row lying out toward a column co-occurs with it more than independence would predict.

Parameters:

axis_x (int, optional) – The axes to plot, counted from 1. Defaults 1 and 2.
axis_y (int, optional) – The axes to plot, counted from 1. Defaults 1 and 2.

settings (dict, optional) –

Default settings:

{
    "title": "Correspondence analysis map",  # str: plot title
    "row_color": None,        # str|None: marker colour for rows
    "column_color": None,     # str|None: marker colour for columns
    "html_image_height": 600, # int: image height in pixels
    "html_aspect_ratio_w_over_h": 1.0,  # float: width over height
    "template": "pi_journal", # str: registered Plotly theme name
}

Return type:

go.Figure

Raises:

ValueError – If an axis is outside 1 .. n_components_.

MCA#

Multiple correspondence analysis: correspondence analysis of the one-hot indicator matrix of several categorical variables, “PCA for categorical data”. Observations with similar combinations of levels plot together, and so do levels that tend to occur together. The raw eigenvalues of an indicator matrix understate the real association, so correction reports each axis’s share on the Benzecri or Greenacre scale.

class process_improve.multivariate.methods.MCA(n_components=2, correction=None)[source]#

Bases: CA

Multiple correspondence analysis of several categorical variables.

“PCA for categorical data”: each observation is described by several categorical variables (grade, supplier, line, shift). MCA is correspondence analysis of the one-hot indicator matrix, so observations with similar combinations of levels plot together, and levels that tend to occur together plot together.

The raw eigenvalues of an indicator matrix are notoriously pessimistic: with $Q$ variables, much of the inertia is an artefact of the coding, and the first axes look as if they explain a small share when they explain nearly all of the real association. correction reports that share honestly.

Parameters:

n_components (int, optional) – Number of axes to keep, default 2. At most $J - Q$ axes carry inertia, with $J$ the total number of levels.
correction ({None, "benzecri", "greenacre"}, optional) –
How to express each axis’s share of the inertia:
- None (default): the raw shares of the indicator matrix’s inertia.
- "benzecri": each eigenvalue above $1/Q$ becomes $(Q/(Q-1))^2 (\lambda - 1/Q)^2$, the rest zero, and shares are taken of their sum. Known to be optimistic.
- "greenacre": the same adjusted eigenvalues as a share of the adjusted total inertia $\frac{Q}{Q-1}(\sum \lambda^2 - (J-Q)/Q^2)$, the sum taken over every axis. More conservative; the shares need not sum to 100%.

corrected_eigenvalues_#

The eigenvalues after correction; the raw ones when it is None.

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

corrected_explained_inertia_#

Share of the (corrected) inertia on each kept axis.

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

variables_#

The categorical variables, in the order they were encoded.

Type:: list[str]

All the attributes of :class:`CA` are present as well. Its "rows" are the

observations and its "columns" the levels, named ``"variable=level"``.

References

M. Greenacre, “From simple to multiple correspondence analysis”, in M. Greenacre and J. Blasius (eds.), Multiple Correspondence Analysis and Related Methods, Chapman and Hall/CRC, 2006.

Examples

>>> mca = MCA(n_components=2, correction="greenacre").fit(batch_attributes)
>>> mca.corrected_explained_inertia_
>>> mca.transform_columns(batch_attributes[["outcome"]])  # a supplementary label

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

Fit to a table of categorical variables, one column per variable.

Parameters:

X (pd.DataFrame of shape (n_observations, n_variables)) – Each column a categorical variable; values of any type are compared as strings.
y (ignored) – Accepted for Pipeline compatibility.

Returns:

self, fitted.

Return type:

MCA

Raises:

ValueError – If correction is not recognised, fewer than two variables are given, or a cell is missing.

transform(X)[source]#

Place observations on the fitted map, from their levels.

Parameters:: X (pd.DataFrame of shape (n_observations, n_variables)) – The same variables as the fit, with levels the fit has seen.
Returns:: Principal coordinates, one row per observation.
Return type:: pd.DataFrame
Raises:: ValueError – If the variables differ from the fit’s, or a level is unseen or missing.

transform_columns(X)[source]#

Place the levels of supplementary categorical variables on the fitted map.

Each level lands at the barycentre of the observations that have it, so a supplementary outcome label (good / bad batch) shows which combinations of the fitted attributes go with each outcome, without the outcome shaping the axes.

Parameters:: X (pd.DataFrame of shape (n_observations, n_supplementary)) – One row per fitted observation, in the same order; one column per supplementary categorical variable.
Returns:: Principal coordinates, one row per supplementary level, named "variable=level".
Return type:: pd.DataFrame

FAMD#

Factor analysis of mixed data: one analysis of numeric and categorical columns together, with each variable given an equal say. It is exactly PCA of the correlation matrix when every column is numeric, and exactly MCA when every column is categorical.

class process_improve.multivariate.methods.FAMD(n_components=2, categorical=None)[source]#

Factor analysis of mixed data: one analysis of numeric and categorical columns together.

Real process records mix continuous readings (temperature, pressure, flow) with categorical context (grade, supplier, line, shift). A PCA cannot take the categorical columns without one-hot encoding them by hand, after which they are over- or under-weighted against the numeric ones, because PCA has no notion of balancing the two. FAMD gives every variable an equal say: each numeric column is standardised, so it contributes an inertia of 1, and each categorical column with $k$ levels contributes $k - 1$, exactly as it would in MCA.

Concretely the numeric columns are centred and scaled to unit (population) variance, each level’s indicator $\delta$ becomes $\delta / \sqrt{p} - \sqrt{p}$ for a level of frequency $p$, and the combined matrix $Z$ is decomposed as $Z / \sqrt{n} = U S V^T$. FAMD is exactly PCA of the correlation matrix when every column is numeric, and exactly MCA (eigenvalues scaled by $Q$) when every column is categorical.

Parameters:

n_components (int, optional) – Number of axes to keep, default 2.
categorical (list of str, optional) – Columns to treat as categorical even though they are numeric, such as a line or shift coded 1, 2, 3. Without it, numeric dtypes are numeric and everything else (strings, categoricals, booleans) is categorical.

n_components_#

Axes actually kept.

Type:: int

eigenvalues_#

Inertia of each kept axis.

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

total_inertia_#

Numeric columns plus, for each categorical column, its levels minus one.

Type:: float

explained_inertia_#

Share of the total inertia on each kept axis.

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

numeric_, categorical_

How each column was treated.

Type:: list[str]

row_coordinates_#

One row per observation, one column per axis.

Type:: pd.DataFrame

column_coordinates_#

One row per numeric column and per categorical level ("variable=level"): for a numeric column, its correlation with the axis.

Type:: pd.DataFrame

row_contributions_, column_contributions_

Share of each axis’s inertia from each observation, and from each numeric column or level; each column sums to 1.

Type:: pd.DataFrame

row_cos2_#

Share of each observation’s own inertia each axis represents.

Type:: pd.DataFrame

References

J. Pages, “Analyse factorielle de donnees mixtes”, Revue de Statistique Appliquee, 52 (2004), 93-111.

Examples

>>> famd = FAMD(n_components=2, categorical=["line"]).fit(batch_record)
>>> famd.explained_inertia_
>>> famd.map_plot()

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

Fit to a table mixing numeric and categorical columns.

Parameters:

X (pd.DataFrame of shape (n_observations, n_columns)) – The mixed table.
y (ignored) – Accepted for Pipeline compatibility.

Returns:

self, fitted.

Return type:

FAMD

Raises:

ValueError – If n_components is not positive, a cell is missing, a numeric column is constant, or categorical names a column that is not there.

transform(X)[source]#

Place new observations on the fitted map.

Parameters:: X (pd.DataFrame) – The same columns as the fit, with levels the fit has seen.
Returns:: One row per observation, one column per axis.
Return type:: pd.DataFrame
Raises:: ValueError – If the columns differ from the fit’s, or a level is unseen or missing.

map_plot(axis_x=1, axis_y=2, settings=None)[source]#

Draw observations and variables on one pair of axes.

The column points are numeric columns (their correlations with the axes, so inside the unit circle) and categorical levels; the observations spread wider. Read a numeric column’s direction, not its distance from an observation.

Parameters:

axis_x (int, optional) – The axes to plot, counted from 1. Defaults 1 and 2.
axis_y (int, optional) – The axes to plot, counted from 1. Defaults 1 and 2.
settings (dict, optional) – As for CA.map_plot(), with the title defaulting to “FAMD map”.

Return type:

go.Figure

MFA#

Multiple factor analysis: one balanced analysis of several groups of variables measured on the same observations. Each group is weighted by one over its first eigenvalue, so a wide block (a spectrum) cannot outvote a narrow one (a lab panel) by column count, and each observation’s position is the average of one partial position per group.

class process_improve.multivariate.methods.MFA(groups=None, n_components=2, scale=True)[source]#

Multiple factor analysis: one balanced analysis of several groups of variables.

Process data arrives in blocks (a spectrum, a lab panel, the process conditions) measured on the same observations. Concatenating them into one PCA lets the widest block win: a 600-wavelength spectrum outvotes a dozen lab results simply by having more columns. MFA first analyses each group on its own and divides it by the square root of its first eigenvalue, so that no group’s leading direction can carry more than an inertia of 1, then runs one PCA on the balanced whole.

Two properties make the result readable. The first eigenvalue lies between 1 and the number of groups: near the number of groups means every block agrees on the dominant direction, near 1 means only one does. And each observation’s position is the average of its partial positions, one per group, so the spread of an observation’s partial points shows how much its blocks disagree about it.

MBPCA answers a related question with a different weighting: it gives each block the same total inertia, where MFA gives each block’s leading direction the same inertia. The two agree exactly when every block’s first eigenvalue is the same share of its total, as it is for one-dimensional blocks.

Parameters:

groups (dict[str, list[str]], optional) – Group name to the columns in it, when the data comes as one frame. Every column must be numeric and appear in exactly one group; columns in no group are ignored. Leave it out to pass the data as a dict of blocks instead, the form MBPCA takes.
n_components (int, optional) – Number of axes to keep, default 2.
scale (bool, optional) – Scale every column to unit variance before the groups are formed, default True. Set False only when the columns within each group share units, a spectrum for instance, where scaling would inflate the noisy wavelengths.

n_components_#

Axes actually kept.

Type:: int

eigenvalues_#

Inertia of each kept axis; the first lies between 1 and the number of groups.

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

total_inertia_#

Sum over every axis.

Type:: float

explained_inertia_#

Share of the total inertia on each kept axis.

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

groups_#

The groups used: groups, or the blocks’ names and columns.

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

group_weights_#

The weight each group’s columns were multiplied by: one over its first eigenvalue.

Type:: pd.Series

row_coordinates_#

One row per observation, one column per axis.

Type:: pd.DataFrame

partial_row_coordinates_#

Per group, where that group alone would place each observation. Their average over groups is row_coordinates_.

Type:: dict[str, pd.DataFrame]

column_coordinates_#

Correlation of each column with each axis (for scale=True).

Type:: pd.DataFrame

group_coordinates_#

One row per group: the inertia it contributes to each axis. Each lies between 0 and 1, and a group near 1 on an axis has that axis as its own leading direction.

Type:: pd.DataFrame

column_contributions_, group_contributions_

Share of each axis’s inertia from each column, and from each group; each column sums to 1.

Type:: pd.DataFrame

References

B. Escofier and J. Pages, “Multiple factor analysis (AFMULT package)”, Computational Statistics and Data Analysis, 18 (1994), 121-140.

Examples

>>> mfa = MFA({"spectra": wavelengths, "lab": assays}).fit(batches)
>>> mfa = MFA().fit({"spectra": spectra, "lab": lab})     # as for MBPCA
>>> mfa.eigenvalues_[0]        # near 2 means both blocks agree
>>> mfa.group_coordinates_

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

Fit to a table whose columns are organised into groups.

Parameters:

X (pd.DataFrame of shape (n_observations, n_columns), or dict[str, pd.DataFrame]) – One frame whose column names match those in groups, or a dict of blocks with the same row index, one per group.
y (ignored) – Accepted for Pipeline compatibility.

Returns:

self, fitted.

Return type:

MFA

Raises:

ValueError – If the groups are invalid or missing, the blocks do not share their rows, a cell is missing, n_components is not positive, or a group has no variation.
TypeError – If groups is given but is not a dict.

transform(X)[source]#

Place new observations on the fitted map.

Parameters:: X (pd.DataFrame, or dict[str, pd.DataFrame]) – Observations with every column the groups name, in either of the forms fit() takes.
Returns:: One row per observation, one column per axis.
Return type:: pd.DataFrame
Raises:: ValueError – If a grouped column is absent or a cell is missing.

map_plot(axis_x=1, axis_y=2, settings=None)[source]#

Draw the observations and the columns’ correlations on one pair of axes.

Parameters:

axis_x (int, optional) – The axes to plot, counted from 1. Defaults 1 and 2.
axis_y (int, optional) – The axes to plot, counted from 1. Defaults 1 and 2.
settings (dict, optional) – As for CA.map_plot(), with the title defaulting to “MFA map”.

Return type:

go.Figure

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 on the blocks:

Matrices in F, Z and Y must all have the same number of rows.
Columns in F must be the same as the rows in D. They are paired by name, in any order; a material that no blend uses has a column of zeros in F.
Conditions in Z may be missing (turning it into an L-shaped data structure).
A blend with nothing observed in F, Z or Y, or a material with nothing observed in D, is fitted with a zero score or loading in that block and no weight in that block’s regressions, as PCA and PLS fit an all-missing row.

Parameters:

n_components (int) – A parameter used to specify the number of components.
d_matrix (dict[str, pd.DataFrame]) – A dictionary containing the properties of each group of materials. Keys are group names; values are DataFrames with properties as columns and materials as rows. This “D” matrix is provided once at construction and reused for fitting, prediction and cross-validation.
max_iter (int, optional) – The maximum number of iterations for the TPLS algorithm. Default is 500. A component that has not converged by then is the state of its last pass, and the fit warns.
skip_f_matrix_preprocessing (bool, optional) – If True, the F (formula) matrices are used as-is, skipping the internal centering and scaling of the F block. Default is False. The caller then owns F’s preprocessing: a column with no variance, which scaling would leave out with a warning, is used as given, without one. A missing cell is still left out, and a column with no observed value is still named.
tol (float, optional) – Relative convergence tolerance for the per-component super-score loop: the component is converged once the norm of the difference between two successive super-score vectors, relative to the norm of the previous one, falls below this. Default is epsqrt (about 1.5e-8), which is the value that was hard-coded before this became a parameter, so the default fit is unchanged. It must lie between machine precision and 10, as on every iterative estimator. It is listed last so that no existing positional argument changed position.

Notes

The input X passed to fit() and diagnose() is a dictionary with 3 keys:

F: Formula matrices (rows = blends, columns = materials). F = {"Group A": df_formulas_a, "Group B": df_formulas_b, ...}
Z: Process conditions - one row per blend, one column per condition, in one or more groups. Each Z group has its own weights and block score, as a block of a multi-block PLS.
Y: Product quality indicators - one row per blend, one column per indicator.

The D matrix (database of material properties) is supplied once at construction via the d_matrix argument; it is not part of X.

n_components_#

The number of components fitted. Everything after fit reads it, not the n_components parameter, so set_params without a refit changes nothing.

Type:: int

max_iter_#

The max_iter the model was fitted with, which display_results() reports.

Type:: int

skip_f_matrix_preprocessing_#

Whether the fit used the F block as given; diagnose() treats new F blocks the same way.

Type:: bool

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, milliseconds and converged lists. converged holds whether the component met tol within max_iter passes, as fitting_info_["converged"] does on MBPLS and MBPCA.

Type:: dict

preproc_#

Nested per-block per-group preprocessors, indexed as preproc_[block][group] where block is "D", "F", "Y" or "Z". With skip_f_matrix_preprocessing the "F" entries are learned but not applied.

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. A column excluded for having no variance is marked absent in every row.

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). Each Z group is a block of its own in the super-level, so several Z groups give one row each, labelled Z["<group>"], before the "F" row.

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 fit for the class’s own bookkeeping (is_fitted_, preproc_, required_blocks_) do carry the underscore.

Example

>>> import numpy as np
>>> import pandas as pd
>>> rng = np.random.default_rng()
>>>
>>> n_props_a, n_props_b = 6, 4            # Two groups of properties: A and B.
>>> n_materials_a, n_materials_b = 12, 8   # Number of materials in each group.
>>> n_formulas = 40                        # Number of formulas in matrix F.
>>> n_outputs = 3
>>> n_conditions = 2
>>>
>>> properties = {
>>>     "Group A": pd.DataFrame(rng.standard_normal((n_materials_a, n_props_a))),
>>>     "Group B": pd.DataFrame(rng.standard_normal((n_materials_b, n_props_b))),
>>> }
>>> formulas = {
>>>     "Group A": pd.DataFrame(rng.standard_normal((n_formulas, n_materials_a))),
>>>     "Group B": pd.DataFrame(rng.standard_normal((n_formulas, n_materials_b))),
>>> }
>>> process_conditions = {"Conditions": pd.DataFrame(rng.standard_normal((n_formulas, n_conditions)))}
>>> quality_indicators = {"Quality":    pd.DataFrame(rng.standard_normal((n_formulas, n_outputs)))}
>>> all_data = {"Z": process_conditions, "F": formulas, "Y": quality_indicators}
>>> estimator = TPLS(n_components=4, d_matrix=properties)
>>> estimator.fit(DataFrameDict(all_data))

property plot: Plot#

model.plot.scores() and model.plot.loadings().

Type:: Score and loading plots of the fitted model

optimal_inversion(rules, *, alternatives=0, leave_out=1)#

Find the inputs a PLS or TPLS model vouches for that meet the rules, at the least cost or closest to the targets.

Parameters:

model (PLS, TPLS or AffineModel) – A fitted model, or its maps as affine_model() read them at rules.conf_level (pass the maps to solve many problems on one model, or to search with maps stored by AffineModel.to_dict()).
rules (InversionRules) – The objective (costs or targets), the windows, bounds, fixed inputs, linear rules, counting rules and limits.
alternatives (int) – How many further answers to find, each with a different set of inputs: the best answer that leaves out at least leave_out of the counted inputs used by the answer and by every alternative before it (a no-good cut on the master problem). Needs counted inputs: a group with a max_count or a min_amount. Give the min_amount, so that a trace of an input does not count as using it.
leave_out (int) – How many of each earlier answer’s counted inputs every alternative must leave out.

Returns:

result – With keys:

feasiblebool: Whether any inputs meet every rule. When False, message says which rule is in the way.
x_newpd.Series: The inputs, in their original units.
y_hatpd.Series: The predicted qualities.
scores, hotellings_t2, t2_limit, spe, spe_limit: How far the answer sits from the data the model was built from, against its limits.
objectivefloat: The cost, or the summed distance to the targets in training standard deviations.
useddict of str to list of str: The inputs used from each group.
constraintspd.DataFrame: Every rule, window and limit, and each input at a bound: its value, bounds, the side it binds on, and its shadow_price, the change in the optimal objective per unit increase of that bound (per_percent: per 1% of it). A rule that does not bind has a shadow price of 0. A binding constraint implied by others (a lot at 100% of its material, held both by its bound and by the material’s total) has NaN: its price is not its own.
kkt_residualfloat: How nearly the optimality conditions hold at the answer; near 0 at an optimum.
proven_optimalbool: The search proved no better answer exists, and the answer is optimal for its inputs: the optimality conditions hold, or a linear program over the limits’ tangent planes at the answer does no better. With no answer: the search proved none exists; False when it stopped before it could.
n_cutsint: The tangent planes the counting rules needed.
messagestr: Why there is no answer, or why it is not proven optimal; empty otherwise.
alternativeslist of sklearn.utils.Bunch: Only when alternatives is asked for: each alternative as a result with the keys above, best first. Its proven_optimal means no input set the earlier answers allow does better. The list ends early, with a result whose feasible is False and whose message says why, when no further input set meets the rules; it is empty when there is no answer at all.

Return type:

Raises:

Warns:

See also

inversion_sensitivity: check the shadow prices by solving again, and trace a limit, counts and prices.

References

M. A. Duran and I. E. Grossmann, “An outer-approximation algorithm for a class of mixed-integer nonlinear programs”, Mathematical Programming, 36 (1986): 307-339.

Examples

>>> rules = InversionRules(
...     costs={"Acetic": 1.0, "H2S": 3.0, "Lactic": 8.0}, windows={"Taste": (20.0, 30.0)},
...     bounds={"Acetic": (4.0, 7.0), "H2S": (3.0, 10.0), "Lactic": (0.8, 2.1)},
... )
>>> result = optimal_inversion(pls, rules)
>>> result.x_new, result.objective, result.constraints

inversion_sensitivity(rules, *, front=None, prices=None, fraction=0.01)#

Parameters:

model (PLS, TPLS or AffineModel) – A fitted model, or its maps, as for optimal_inversion().
rules (InversionRules) – The rules, as for optimal_inversion().
front ((str, list of float), optional) – A rule to trace, with the values its bound takes: "T²", an SPE limit ("SPE X"), a quality with a window, or a linear rule’s name. A two-sided rule moves the side that binds. The optimal objective at each.
prices (list of str, optional) – Inputs whose cost to range (costs objective only): how far, as a fraction of itself, each can fall or rise before the answer changes shape (another input used, another rule binding). Each costs about a dozen solves.
fraction (float, default 0.01) – How far each binding rule is moved to check its shadow price, as a fraction of its bound.

Returns:

Return type:

Examples

>>> report = inversion_sensitivity(pls, rules, front=("T²", [4.0, 6.0, 8.0]), prices=["Lactic"])
>>> report.checked   # shadow prices against solving again
>>> report.front     # the cost as the T² limit moves

hotellings_t2_limit(conf_level=0.95, *, judged_rows=None)[source]#

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

judged_rows is "fitted" to judge the rows the model was fitted to (its hotellings_t2), or "new" to judge rows it was not fitted to. None gives the "new" limit until 3.0 and the "fitted" one from 3.0, with a FutureWarning.

Parameters:

conf_level (float)
judged_rows (Literal['fitted', 'new'] | None)

Return type:

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

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

Parameters:

score_horiz (int)
score_vert (int)
conf_level (float)
n_points (int)
judged_rows (Literal['fitted', 'new'] | None)

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 (DataFrameDict) – The training input samples. A plain dict raises TypeError. Keys that must be present: "F", "Z", and "Y". See the class docstring for the shape of each block.
y (object, optional) – Must be None. The signature exists only for sklearn API compatibility; a T-shaped model takes its response from X["Y"]. Anything else raises, rather than being silently ignored as it was before #565.

Returns:

self – Returns self.

Return type:

Raises:

ValueError – If y is not None, or if n_components, max_iter or tol is out of range, or if a material of a group of F is not a row of that group of d_matrix, or the reverse.
TypeError – If X is not a DataFrameDict, or d_matrix is not a dict of DataFrames.
NotEnoughVarianceError – If the first Y block has no variance, or no group has a column of F and a column of D that vary, as MBPLS raises for an X with no variance; if the first Y block has no variance left for a component, as when the components before it fit Y exactly, in the words of PLS and MBPLS; or if a component’s scores are not finite because a block has no variation for the iteration to follow. A group whose F or D has no variation, beside one whose F and D vary, is fitted with zero weight (#1144): its weights and loadings are 0, and a warning names it.

Warns:

SpecificationWarning – If a component did not converge within max_iter passes (#1163).

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. A plain dict raises TypeError, so wrap it. # You need at least the F block for a new prediction; “Z” is optional. new_data = DataFrameDict({“Z”: …, “F”: …}) predictions = estimator.diagnose(new_data)

Parameters:

X (DataFrameDict) – The input samples.

Returns:

A bunch with the following fields:

hat (dict[str, DataFrame]) : predicted Y per Y-group, on the original (un-scaled) scale, indexed by the observation names of X.
t_scores_super (DataFrame, shape (n_new, n_components)) : super-scores for the new observations.
spe (dict[str, dict[str, DataFrame]]) : per-component squared prediction error for the "Z" and "F" blocks, keyed first by block name and then by group; each inner DataFrame has shape (n_new, n_components).
cumulative_hotellings_t2 (DataFrame, shape (n_new, n_components)) : cumulative Hotelling’s T² per new observation after each component, as in hotellings_t2, the fitted attribute. Every diagnose of the family returns it in this shape (#983).
hotellings_t2 : deprecated since 2.4.0, to be removed in 3.0; the same DataFrame as cumulative_hotellings_t2. Reading it warns.
residuals (dict[str, dict[str, DataFrame]]) : what is left of each "Z" and "F" group after all the components, in the preprocessed (centred and scaled) units, keyed like spe. The last column of spe is the length of each row; a missing cell, and a column excluded for having no variance, leave 0.

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 r2_score() of the Y block on test data.

See RegressorMixin.score for the general contract.

Parameters:

X (DataFrameDict) – Test samples. The nested "Y" block supplies the actual response values; X["Z"] and X["F"] drive the prediction.
y (object, optional) – Must be None. The Y-data comes from X["Y"], not from a separate argument (the sklearn RegressorMixin.score signature is preserved only for API compatibility). Anything else raises, rather than being silently ignored as it was before #565.
sample_weight (np.ndarray or None) – Optional per-sample weight forwarded to r2_score().

Returns:

score – r2_score() of X["Y"], the one Y block a TPLS model fits, against self.diagnose(X).hat: the mean of the columns’ $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. Each component is weighted by its share of that block’s own variance (its per-column R2, averaged over the columns of each group and then over the groups with a column that varies), not of Y as a PLS VIP is. A group with no variation has zero weight, and a VIP of 0.

Parameters:

block (str or None, default=None) – Which block to return. Must be "D" or "F", or None to return all blocks.
method ({"vip", "deflated"}, default="vip") –
How importance is measured.
- "vip" (default): standard VIP on the raw block loadings, as reported by Garcia-Munoz (2014). These are the values stored in 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: a four-key dict keyed by block, of the shape {"D": {group: pd.Series, ...}, "F": {group: pd.Series, ...}, "Z": {}, "Y": {}}. VIP is defined only for the D- and F-blocks, so the "Z" and "Y" entries are always empty dicts; they are kept for shape-stability across every vip() caller. If block is "D" or "F": the inner dict {group: pd.Series, ...} for that block, where each pd.Series is indexed by feature names.

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])

Optimal inversion#

The inputs a PLS or TPLS model vouches for, at the least cost or closest to target qualities, under windows, bounds, linear rules (mixture totals, safety limits), the T² and SPE limits and counting rules; with the shadow price of every rule. See Model inversion and the orthogonal space.

process_improve.multivariate.methods.optimal_inversion(model, rules, *, alternatives=0, leave_out=1)[source]#

Find the inputs a PLS or TPLS model vouches for that meet the rules, at the least cost or closest to the targets.

Parameters:

model (PLS, TPLS or AffineModel) – A fitted model, or its maps as affine_model() read them at rules.conf_level (pass the maps to solve many problems on one model, or to search with maps stored by AffineModel.to_dict()).
rules (InversionRules) – The objective (costs or targets), the windows, bounds, fixed inputs, linear rules, counting rules and limits.
alternatives (int) – How many further answers to find, each with a different set of inputs: the best answer that leaves out at least leave_out of the counted inputs used by the answer and by every alternative before it (a no-good cut on the master problem). Needs counted inputs: a group with a max_count or a min_amount. Give the min_amount, so that a trace of an input does not count as using it.
leave_out (int) – How many of each earlier answer’s counted inputs every alternative must leave out.

Returns:

result – With keys:

feasiblebool: Whether any inputs meet every rule. When False, message says which rule is in the way.
x_newpd.Series: The inputs, in their original units.
y_hatpd.Series: The predicted qualities.
scores, hotellings_t2, t2_limit, spe, spe_limit: How far the answer sits from the data the model was built from, against its limits.
objectivefloat: The cost, or the summed distance to the targets in training standard deviations.
useddict of str to list of str: The inputs used from each group.
constraintspd.DataFrame: Every rule, window and limit, and each input at a bound: its value, bounds, the side it binds on, and its shadow_price, the change in the optimal objective per unit increase of that bound (per_percent: per 1% of it). A rule that does not bind has a shadow price of 0. A binding constraint implied by others (a lot at 100% of its material, held both by its bound and by the material’s total) has NaN: its price is not its own.
kkt_residualfloat: How nearly the optimality conditions hold at the answer; near 0 at an optimum.
proven_optimalbool: The search proved no better answer exists, and the answer is optimal for its inputs: the optimality conditions hold, or a linear program over the limits’ tangent planes at the answer does no better. With no answer: the search proved none exists; False when it stopped before it could.
n_cutsint: The tangent planes the counting rules needed.
messagestr: Why there is no answer, or why it is not proven optimal; empty otherwise.
alternativeslist of sklearn.utils.Bunch: Only when alternatives is asked for: each alternative as a result with the keys above, best first. Its proven_optimal means no input set the earlier answers allow does better. The list ends early, with a result whose feasible is False and whose message says why, when no further input set meets the rules; it is empty when there is no answer at all.

Return type:

Raises:

Warns:

See also

inversion_sensitivity: check the shadow prices by solving again, and trace a limit, counts and prices.

References

M. A. Duran and I. E. Grossmann, “An outer-approximation algorithm for a class of mixed-integer nonlinear programs”, Mathematical Programming, 36 (1986): 307-339.

Examples

>>> rules = InversionRules(
...     costs={"Acetic": 1.0, "H2S": 3.0, "Lactic": 8.0}, windows={"Taste": (20.0, 30.0)},
...     bounds={"Acetic": (4.0, 7.0), "H2S": (3.0, 10.0), "Lactic": (0.8, 2.1)},
... )
>>> result = optimal_inversion(pls, rules)
>>> result.x_new, result.objective, result.constraints

process_improve.multivariate.methods.inversion_sensitivity(model, rules, *, front=None, prices=None, fraction=0.01)[source]#

Parameters:

model (PLS, TPLS or AffineModel) – A fitted model, or its maps, as for optimal_inversion().
rules (InversionRules) – The rules, as for optimal_inversion().
front ((str, list of float), optional) – A rule to trace, with the values its bound takes: "T²", an SPE limit ("SPE X"), a quality with a window, or a linear rule’s name. A two-sided rule moves the side that binds. The optimal objective at each.
prices (list of str, optional) – Inputs whose cost to range (costs objective only): how far, as a fraction of itself, each can fall or rise before the answer changes shape (another input used, another rule binding). Each costs about a dozen solves.
fraction (float, default 0.01) – How far each binding rule is moved to check its shadow price, as a fraction of its bound.

Returns:

Return type:

Examples

>>> report = inversion_sensitivity(pls, rules, front=("T²", [4.0, 6.0, 8.0]), prices=["Lactic"])
>>> report.checked   # shadow prices against solving again
>>> report.front     # the cost as the T² limit moves

class process_improve.multivariate.methods.InversionRules(costs=<factory>, targets=<factory>, windows=<factory>, bounds=<factory>, fixed=<factory>, rules=(), closure=True, groups=<factory>, max_count=<factory>, min_amount=<factory>, conf_level=0.95, t2_limit=None, spe_limits=<factory>)[source]#

Bases: object

What the inputs must meet, and what to minimise.

Give costs or targets, not both. To keep a quality near its target while minimising cost, give it a window instead.

Parameters:

costs (Mapping[str, float])
targets (Mapping[str, float])
windows (Mapping[str, tuple[float | None, float | None]])
bounds (Mapping[str, tuple[float | None, float | None]])
fixed (Mapping[str, float])
rules (tuple[LinearRule, ...])
closure (bool)
groups (Mapping[str, list[str]])
max_count (Mapping[str, int])
min_amount (Mapping[str, float])
conf_level (float)
t2_limit (float | None)
spe_limits (Mapping[str, float])

costs#

Each input’s cost per unit: minimise their sum. Inputs not named cost nothing.

Type:: dict of str to float

targets#

Target qualities: minimise the summed distance to them, each in its training standard deviations.

Type:: dict of str to float

windows#

Each quality’s predicted value between a low and a high bound; None for no bound.

Type:: dict of str to (float or None, float or None)

bounds#

Each input’s lower and upper bound: its extreme limits, the values it can take. A TPLS formula’s fractions are 0 to 1 unless given. Give every other input both: the T² and SPE limits alone let an answer set an input to a value it cannot take, such as a negative amount. An input without both is deprecated since 2.4.0 (a FutureWarning names it) and an error from 3.0.

Type:: dict of str to (float or None, float or None)

fixed#

Inputs held at a value.

Type:: dict of str to float

rules#

Linear rules on the inputs: mixture totals, safety and legal limits.

Type:: tuple of LinearRule

closure#

For TPLS, each group of lots sums to 1, as formulas do.

Type:: bool

groups#

Groups of inputs for the counting rules, over and above the model’s own (a TPLS model’s materials).

Type:: dict of str to list of str

max_count#

The most inputs a group may use.

Type:: dict of str to int

min_amount#

The least amount of an input of the group, when used.

Type:: dict of str to float

conf_level#

The confidence level of the T² and SPE limits.

Type:: float

t2_limit#

A T² limit in place of the one at conf_level.

Type:: float or None

spe_limits#

SPE limits in place of those at conf_level, by block name ("X", "F/Group1").

Type:: dict of str to float

class process_improve.multivariate.methods.LinearRule(name, coefficients, lower=None, upper=None)[source]

Bases: object

A linear rule on the inputs: lower <= sum(coefficients[name] * x[name]) <= upper.

A mixture that sums to 1, a component’s content below its legal limit (each ingredient’s content as its coefficient), or a ratio held between two inputs, are each one rule.

Parameters:

name (str)
coefficients (Mapping[str, float])
lower (float | None)
upper (float | None)

name

How the rule is reported, and how inversion_sensitivity() names it.

Type:: str

coefficients

Each input’s coefficient; inputs not named have 0.

Type:: dict of str to float

lower, upper

The bounds; None for none.

Type:: float or None

process_improve.multivariate.methods.affine_model(model, conf_level=0.95)[source]#

Read a fitted PLS or TPLS model as exact affine maps of one input vector, checked against its diagnose.

The scores, predictions and residuals of a latent-variable model are affine in its inputs. One diagnose call, at the model’s centre and one step along each input, gives every map exactly. A second call, at random points, is the check: the maps must reproduce the model’s own predictions, Hotelling’s T² and SPE there.

Parameters:

model (PLS or TPLS) – A fitted model. For TPLS the inputs are the formulas’ fractions (each group’s lots) and the process conditions.
conf_level (float, default 0.95) – The confidence level of the T² and SPE limits.

Returns:

The maps, the limits, and each input’s training mean and spread.

Return type:

AffineModel

Raises:

TypeError – The model is not a PLS or TPLS model.
ValueError – The maps do not reproduce the model’s diagnose: the model is not one this reads exactly.

class process_improve.multivariate.methods.AffineModel(names, y_names, groups, y0, y_map, t0, t_map, score_sd, t2_limit, blocks, y_scale, centre, spread, conf_level=0.95)[source]#

Bases: object

A fitted latent-variable model read as exact affine maps of one vector of inputs, in their original units.

Parameters:

names (tuple[str, ...])
y_names (tuple[str, ...])
groups (Mapping[str, tuple[str, ...]])
y0 (ndarray)
y_map (ndarray)
t0 (ndarray)
t_map (ndarray)
score_sd (ndarray)
t2_limit (float)
blocks (tuple[ResidualBlock, ...])
y_scale (ndarray)
centre (ndarray)
spread (ndarray)
conf_level (float)

names#

The inputs, in the order of the maps’ columns.

Type:: tuple of str

y_names#

The predicted qualities.

Type:: tuple of str

groups#

Natural groups of inputs (a TPLS model’s materials, whose fractions sum to 1); empty for PLS.

Type:: dict of str to tuple of str

y0, y_map

Predictions: y = y0 + y_map @ x.

Type:: np.ndarray

t0, t_map

Scores: t = t0 + t_map @ x.

Type:: np.ndarray

score_sd#

Each score’s standard deviation, as the model’s Hotelling’s T² divides by it.

Type:: np.ndarray

t2_limit#

Hotelling’s T² limit at the confidence level read.

Type:: float

blocks#

Each block’s residual, whose length is its SPE.

Type:: tuple of ResidualBlock

y_scale#

Each quality’s training standard deviation: the unit of a distance to a target.

Type:: np.ndarray

centre, spread

Each input’s training mean and standard deviation (1 where unknown): the start of a search and its scale.

Type:: np.ndarray

conf_level#

The confidence level the limits were read at.

Type:: float

predict(x)[source]#

Return the predicted qualities of one input vector.

Parameters:: x (ndarray)
Return type:: ndarray

scores(x)[source]#

Return the scores of one input vector.

Parameters:: x (ndarray)
Return type:: ndarray

hotellings_t2(x)[source]#

Return Hotelling’s T² of one input vector.

Parameters:: x (ndarray)
Return type:: float

spe(x)[source]#

Return each block’s SPE for one input vector.

Parameters:: x (ndarray)
Return type:: ndarray

to_dict()[source]#

Return the maps as plain JSON: lists and numbers, to store beside a project and search later without the model.

Returns:: Every field; from_dict() reads it back.
Return type:: dict

classmethod from_dict(values)[source]#

Read the maps back from to_dict().

Parameters:: values (dict) – As to_dict() wrote it.
Returns:: The same maps, as arrays.
Return type:: AffineModel

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_#

The intercept removed before the decomposition: the column means on a balanced design, and on an unbalanced one the intercept of the sum-coded least-squares fit, which is not the column mean.

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, and equally often: a combination no row takes leaves the design unbalanced.

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 ssq_ does not partition the total, to rounding; the gap is reported as ssq_["total"] minus the sum of the parts. A design whose terms are orthogonal without being balanced, such as a fraction fitted for its main effects, partitions the total and does not warn.

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:

FPCA#

Functional PCA of curves: a time course, a profile or a spectrum per sample. Each curve is smoothed by one penalised B-spline smoother, and the components are found in the L2 metric of the basis, so they are orthonormal functions of the curve’s argument. Curves sampled at their own points share the union of the points as their grid, with NaN where a curve was not measured. The scores are ordinary numbers, ready to model against the factors of a designed experiment.

class process_improve.multivariate.methods.FPCA(n_components=None, *, n_basis=None, smoothing=None, degree=3, penalty_order=2)[source]#

Functional principal component analysis of curves measured on a grid.

Each curve is smoothed by a cubic penalised B-spline (P-spline): n_basis B-splines on equally spaced knots that run past the ends of the grid, with a penalty on the penalty_order-th differences of their coefficients. One penalty is shared by every curve, so the smoother is one linear map and cannot be mistaken for a difference between curves. The principal components of the smoothed curves are then found in the L2 metric of the basis: the coefficients are weighted by the square root of the Gram matrix W (the integral of the products of the B-splines), so the eigenfunctions are orthonormal as functions and do not depend on where the knots are placed (Ramsay and Silverman, Functional Data Analysis, 2005, section 8.4).

Parameters:

n_components (int or None, default=None) – Number of components. None fits as many as the curves support: one fewer than the number of curves, or n_basis, whichever is smaller. A larger request is clamped to that, with a SpecificationWarning.
n_basis (int or None, default=None) – Number of B-splines. None takes the number of grid points plus two, at most 40: the penalty, not the number of B-splines, sets how smooth the curves are.
smoothing (float or None, default=None) – The penalty’s weight. None chooses it by generalised cross-validation, pooled over all the curves.
degree (int, default=3) – Degree of the B-splines (3 is cubic).
penalty_order (int, default=2) – Order of the coefficient differences penalised. Raise it to 3 when rates (first derivatives) are what is wanted, so the rates are penalised for curvature, not the curves.

n_components_#

The number of components fitted.

Type:: int

grid_#

The grid the curves are measured on.

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

n_basis_#

The number of B-splines used.

Type:: int

knots_#

The knots of the B-splines.

Type:: np.ndarray

smoothing_#

The penalty’s weight, as given or as chosen by cross-validation.

Type:: float

mean_function_#

The mean curve on grid_.

Type:: pd.Series

eigenfunctions_#

The eigenfunctions (functional principal components) on grid_: orthonormal in L2, so the integral of each one squared is 1.

Type:: pd.DataFrame of shape (n_points, n_components)

scores_#

Each curve’s score on each component.

Type:: pd.DataFrame of shape (n_curves, n_components)

explained_variance_#

The variance of the scores of each component: the eigenvalues of the covariance of the smoothed curves.

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

explained_variance_ratio_#

The fraction of the smoothed curves’ variance each component explains.

Type:: pd.Series of length n_components

pca_#

The PCA of the weighted, centred coefficients the components come from. Its hotellings_t2_ and spe_ flag a curve unlike the others.

Type:: PCA

n_features_in_#

The number of grid points.

Type:: int

Examples

>>> import numpy as np, pandas as pd
>>> t = np.linspace(0, 1, 25)
>>> rng = np.random.default_rng(0)
>>> sizes = rng.uniform(1, 2, 12)
>>> curves = pd.DataFrame([a * np.sin(np.pi * t) + rng.normal(0, 0.01, t.size) for a in sizes], columns=t)
>>> model = FPCA(n_components=2).fit(curves)
>>> bool(model.explained_variance_ratio_.iloc[0] > 0.99)
True

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

Smooth the curves and find their principal components.

Parameters:

X (pd.DataFrame or np.ndarray of shape (n_curves, n_points)) – One curve per row, one grid point per column; NaN where a curve was not measured.
y (None) – Ignored; present for the scikit-learn API.
grid (np.ndarray of shape (n_points,), optional) – The grid, strictly increasing. Taken from the columns of a DataFrame when they are numbers (or text that reads as numbers) and grid is not given.

Returns:

The fitted model.

Return type:

FPCA

transform(X)[source]#

Return the scores of new curves, measured on the grid the model was fitted to.

Parameters:: X (pd.DataFrame or np.ndarray of shape (n_curves, n_points)) – One curve per row, on grid_; NaN where a curve was not measured.
Returns:: The scores.
Return type:: pd.DataFrame of shape (n_curves, n_components)

inverse_transform(scores, *, grid=None, derivative=0)[source]#

Return the curves that these scores describe: the mean plus each score times its eigenfunction.

Parameters:

scores (pd.DataFrame or np.ndarray of shape (n_curves, n_components)) – Scores on the fitted components.
grid (np.ndarray, optional) – Where to evaluate the curves, within the fitted grid’s range. grid_ when not given.
derivative (int, default=0) – Order of the derivative to return: 1 gives each curve’s rate of change, read from the B-splines exactly, not by differencing.

Returns:

One curve per row.

Return type:

pd.DataFrame of shape (n_curves, len(grid))

set_fit_request(*, grid='$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:

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

Returns:

self – The updated object.

Return type:

set_inverse_transform_request(*, derivative='$UNCHANGED$', grid='$UNCHANGED$', scores='$UNCHANGED$')#

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

The options for each parameter are:

True: metadata is requested, and passed to inverse_transform 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 inverse_transform.
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:

derivative (str, True, False, or None, default=sklearn.utils.metadata_routing.UNCHANGED) – Metadata routing for derivative parameter in inverse_transform.
grid (str, True, False, or None, default=sklearn.utils.metadata_routing.UNCHANGED) – Metadata routing for grid parameter in inverse_transform.
scores (str, True, False, or None, default=sklearn.utils.metadata_routing.UNCHANGED) – Metadata routing for scores parameter in inverse_transform.
self (FPCA)

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=np.float64(1.4901161193847656e-08), 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, default=``epsqrt`` (about 1.5e-8)) – 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). The same default, and the same meaning, as on PCA / PLS / TPLS. It used to default to None and be resolved to epsqrt inside fit; the literal is the default now, which removes a branch and lets get_params report the value that will actually be used (#588). The legacy absolute tolerance np.finfo(float).eps ** (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. A row observed in no block at all has a super score of 0 and no weight, as PCA and PLS fit an all-missing row; so does a row with no observed response in the regressions of the blocks on Y (#1098).
missing_data_settings (dict or None, default=None) – MBPLS reads no keys from it. The keys md_tol and md_max_iter were removed in 2.0 and raise a ValueError; set tol and max_iter instead.
fitting) (Attributes (after)
--------------------------
n_components – The number of components fitted. Everything after fit reads it, not the n_components parameter, so set_params without a refit changes nothing.
block_names (list[str]) – Ordered list of X-block names (the keys of the input dict).
block_widths (dict[str, int]) – Number of variables in each X-block.
block_informative_widths (dict[str, int]) – Number of variables in each X-block whose values vary (at least 1). The block weighting divides by its square root, so constant columns do not lower a block’s weight.
n_samples (int) – Number of rows fitted.
n_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, and 0 for a row with nothing observed 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. Computed from super_weights_, each component weighted by its share of the Y variance (r2_y_per_component_), as a PLS VIP is: it ranks the blocks by their part in predicting Y.
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. Each component is scaled by its block scores’ second moment about zero, as super_hotellings_t2_ is by explained_variance_.
block_vip (dict[str, pd.Series]) – Per-block variable-importance in projection, indexed by variable name inside each block. Computed from block_weights_, each component weighted by its share of that block’s X variance (r2_x_per_block_per_component_), not of Y as super_vip_ is: it ranks a block’s variables by their part in the components that describe the 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).

A column of an X-block or of Y with nothing observed is left out of the model, with a UserWarning that names it, as across the family (#1098): its weights and loadings are 0, a value a new row has in it adds nothing to that row’s scores or SPE, and a response with no observed value has NaN predictions.

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.

Raises:

NotEnoughVarianceError – If every column of Y, or every column of every X-block, is constant or has no observed value, or if Y has no variance left for a component, as when it is fitted exactly: as PLS.fit() raises. A block with no variance beside one that varies is fitted with zero weight: its block scores, weights and contributions are 0, and a warning names it.

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

A block a row has nothing observed in leaves no residual in that row: its NaN SPE counts as 0 in the sum.

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.

A row with missing cells is scored over the cells it has, as in diagnose(): its super score divides by the weights’ sum of squares over those cells, and so does each of its contributions, which then still add up to it (#1144). They are NaN at the missing cells, which add nothing to the sum.

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 cumulative_hotellings_t2 (DataFrame, n_samples x n_components: column a holds the T² over the first a components, as in super_hotellings_t2_; every diagnose of the family returns it in this shape, #983). A row with nothing observed in a block has NaN for its block scores and its block SPE there, as in block_scores_ and block_spe_.

hotellings_t2, the Series of the full model’s T², is deprecated since 2.4.0 and is removed in 3.0; it is cumulative_hotellings_t2.iloc[:, -1], and reading it warns.

The rename (since 1.38.4, #395) matches PLS.diagnose() and PCA.diagnose. The deprecated predict alias, which returned this Bunch, was removed in 2.0; predict() now returns the predictions only.

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

predict(X)[source]#

Predict Y for new observations, on its original scale.

Returns the predictions only, so the model meets scikit-learn’s regressor contract (score(), cross_val_predict(), a Pipeline step), as PLS.predict() does. For the super scores, the block scores, the block SPE and Hotelling’s T² as well, see diagnose().

Parameters:: X (dict[str, pd.DataFrame]) – The X-blocks of the new observations, under the names they were fitted with.
Returns:: The predicted Y, indexed by the rows of X and labelled with the fitted Y columns: diagnose(X).predictions.
Return type:: pd.DataFrame of shape (n_samples, n_targets)

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

Return the R² of predict() on X against y.

Parameters:

X (dict[str, pd.DataFrame]) – The X-blocks of the observations.
y (pd.DataFrame of shape (n_samples, n_targets)) – Their observed Y.
sample_weight (np.ndarray of shape (n_samples,) or None) – Weights of the observations.

Returns:

R² of the predictions, averaged over the Y columns with equal weights, as PLS.score() gives it.

Return type:

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, numpy.random.Generator or None, optional) – Seed for the shuffling; a Generator gives its next draw of an int, as PLS.select_n_components() takes it. 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), selection_rule, and n_splits (the per-repeat fold count: cv when it was given as an integer, else the number of splits produced by the supplied splitter).

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

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, *, random_state=None, **deprecated)[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 the model to them with its own settings (a clone of model, with as many components as it was fitted with), 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.
random_state (int, numpy.random.Generator or None, default=None) – Draws the permutations. An int draws the same ones on every call, a Generator draws from its own state, and None draws fresh ones.
**deprecated – Takes seed, the old name of random_state, deprecated since 2.4.0 and removed in 3.0. Passing it warns; passing both raises a TypeError, as does any other name.

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 200 permutations cannot resolve below ~0.5%; raise n_permutations for finer resolution. 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=np.float64(1.4901161193847656e-08), 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, default=``epsqrt`` (about 1.5e-8)) – Relative convergence tolerance on the super-score change: the norm of the change between two successive super-score iterations, divided by the norm of the current super-score vector (#504). The same default, and the same meaning, as on PCA / PLS / TPLS. It used to default to None and be resolved to epsqrt inside fit; the literal is the default now, which removes a branch and lets get_params report the value that will actually be used (#588). The legacy absolute tolerance np.finfo(float).eps ** (9/10) (about 8.2e-15) sits below the floating-point oscillation floor of a relative criterion, so it would never be reached in practice.
algorithm (str) –
Algorithm to use for fitting the model.
- "auto": dense vectorised hierarchical NIPALS when the data is complete; mask-aware NIPALS (NaN-tolerant) when any block contains missing values.
- "dense": dense vectorised hierarchical NIPALS. Raises if any block contains missing values.
- "nipals": mask-aware hierarchical NIPALS. Always uses the NaN-tolerant inner-loop primitives, even when the data is complete (slower than "dense" but produces equivalent results).
missing_data_settings (dict or None, default=None) – MBPCA reads no keys from it. The keys md_tol and md_max_iter were removed in 2.0 and raise a ValueError; set tol and max_iter instead.
fitting) (Attributes (after)
--------------------------
n_components – The number of components fitted. Everything after fit reads it, not the n_components parameter, so set_params without a refit changes nothing.
block_names (list[str]) – Ordered list of X-block names (the keys of the input dict).
block_widths (dict[str, int]) – Number of variables in each X-block.
block_informative_widths (dict[str, int]) – Number of variables in each block whose values vary (at least 1). The block weighting divides by its square root, so constant columns do not lower a block’s weight.
n_samples (int) – Number of rows fitted.
n_features_in (int) – Total number of X variables summed across blocks.
feature_names_in (np.ndarray) – Concatenated column names, one per feature, in block order.
preproc (dict[str, MCUVScaler]) – Per-block preprocessors used to mean-centre and unit-variance scale each X-block.
super_scores (pd.DataFrame, shape (n_samples, n_components)) – Super-block (consensus) scores T.
super_loadings (pd.DataFrame, shape (n_blocks, n_components)) – Super-block loadings p_super; rows indexed by block name.
super_hotellings_t2 (pd.DataFrame, shape (n_samples, n_components)) – Cumulative Hotelling’s T^2 on the super-scores per component.
block_scores (dict[str, pd.DataFrame]) – Per-block scores t_b, each shape (n_samples, n_components). NaN for a row with nothing observed in that block, as in MBPLS.
block_loadings (dict[str, pd.DataFrame]) – Per-block loadings p_b, each shape (K_b, n_components).
block_spe (dict[str, pd.DataFrame]) – Per-block squared prediction error per sample and component. NaN for a row with nothing observed in that block.
block_hotellings_t2 (dict[str, pd.DataFrame]) – Per-block cumulative Hotelling’s T^2 per sample and component. Each component is scaled by its block scores’ second moment about zero, as super_hotellings_t2_ is by explained_variance_.
block_vip (dict[str, pd.Series]) – Per-block variable-importance in projection, indexed by variable name inside each block. Computed from block_loadings_, each component weighted by its share of that block’s X variance (r2_x_per_block_per_component_), as a PCA VIP is.
r2_x_per_block_cumulative (pd.DataFrame, shape (n_blocks, n_components)) – Cumulative R^2X per block and component.
r2_x_per_block_per_component (pd.DataFrame, shape (n_blocks, n_components)) – Incremental R^2X per block and component.
r2_x_per_variable (dict[str, pd.DataFrame]) – Cumulative R^2X per variable within each block.
explained_variance (np.ndarray, shape (n_components,)) – Variance of the super-score per component (ddof=1).
scaling_factor_for_super_scores (pd.Series) – sqrt(explained_variance_) per component.
fitting_info (dict) – Per-component iteration count and timing.
has_missing_data (bool) – Whether any X-block had NaN values.
algorithm – 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).

A row with nothing observed in a block has no score for that block, and no say in the super loading; its super score comes from the blocks it does have. A row with nothing observed in any block is fitted with a super score of 0 and no weight, as PCA, PLS and TPLS fit an all-missing row. A column with nothing observed is left out of the model, with a UserWarning that names it: its loadings are 0, and a value a new row has in it adds nothing to that row’s scores or SPE (#1098).

transform() and diagnose() score a new row as the fit scores its own: over the cells the row has, so a training row is given its fitted scores again (#1076). Its block SPE is over those cells too, and NaN for a block it has nothing observed in. A new row with nothing observed at all has NaN scores, SPE and T², where the fit places such a row at 0.

References

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

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

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

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.

Raises:

NotEnoughVarianceError – If every column of every block is constant or has no observed value, as across the family (#1144). A block with no variance beside one that varies is fitted with zero weight: its block scores, loadings and contributions are 0, and a warning names it.

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

cumulative_hotellings_t2 is a DataFrame, n_samples x n_components: column a holds the T² over the first a components, as in super_hotellings_t2_; every diagnose of the family returns it in this shape (#983). hotellings_t2, the Series of the full model’s T², is deprecated since 2.4.0 and is removed in 3.0; it is cumulative_hotellings_t2.iloc[:, -1], and reading it warns.

A row with missing cells is scored over the cells it has, as the fit scores a row. A row with nothing observed in a block has NaN for its block scores and its block SPE there, as in block_scores_ and block_spe_; a row with nothing observed at all has NaN scores and T² as well.

The rename (since 1.38.4, #395) matches PCA.diagnose() and PLS.diagnose(). The deprecated predict alias was removed in 2.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.

A block a row has nothing observed in leaves no residual in that row: its NaN SPE counts as 0 in the sum, as in MBPLS.super_spe_limit().

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.

A row with missing cells is scored over the cells it has, as in diagnose(): $p_b^\top p_b$ is then taken over the row’s cells in block $b$, and $p_\mathrm{super}^\top p_\mathrm{super}$ over the blocks it has a score in, so its contributions still add up to its super score (#1144). They are NaN at the missing cells, which add nothing to the sum.

See 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

GPA#

Generalized Procrustes analysis: several configurations of the same objects (an assessor’s scores, an instrument’s measurements) are centred, rotated and scaled onto one consensus. The residuals show which configuration disagrees about which object, and GPA.consensus_test() shows whether the consensus is more than alignment alone would produce. Takes the same dict-of-DataFrames input as MBPCA; GPA.configurations_from_long() builds it from a sensory panel.

class process_improve.multivariate.methods.GPA(scale=True, tol=1e-10, max_iter=100)[source]#

Bases: BaseEstimator

Generalized Procrustes analysis: a consensus from several configurations of the same objects.

Each configuration places the same objects (products, parts, batches) as points, one column per dimension: an assessor’s scores on their attributes, an instrument’s measurements. GPA centres every configuration, rotates it (reflections allowed) and, with scale=True, stretches or shrinks it as a whole, so that together they sit as close as possible to their mean, the consensus (Gower, 1975; ten Berge, 1977). The configurations may have different numbers of columns, as in free-choice profiling where each assessor uses their own words; the narrower ones are padded with zeros.

GPA always finds a consensus, even among unrelated configurations: four of pure noise typically share around half their variation after alignment. Judge the consensus with consensus_test(), not with consensus_share_ alone.

Parameters:

scale (bool, optional) – Give each configuration its own isotropic scaling factor, default True. Keep it for panels, where assessors use different ranges of the scale; set False when the configurations share units and a difference in size is itself a finding.
tol (float, optional) – Stop when a round of rotations and scalings lowers the residual sum of squares by less than this share of the total, default 1e-10.
max_iter (int, optional) – Most rounds to run, default 100. Convergence usually takes under twenty.

consensus_#

The mean of the aligned configurations, turned to its principal axes, so the first column carries the most variation.

Type:: pd.DataFrame of shape (n_objects, n_dimensions)

aligned_#

Each configuration after translation, rotation and scaling, in the consensus’s axes.

Type:: dict[str, pd.DataFrame]

scale_factors_#

The factor each centred configuration was multiplied by. Above 1, it was stretched: the assessor used a narrower range than the others. The factors preserve the total sum of squares of the centred configurations.

Type:: pd.Series

rotations_#

The orthogonal matrix, of shape (n_dimensions, n_dimensions), applied to each centred and padded configuration.

Type:: dict[str, np.ndarray]

residuals_#

Squared distance of each object from the consensus, in each configuration. Summed down a column, how far a configuration disagrees; along a row, how contested an object is.

Type:: pd.DataFrame of shape (n_objects, n_configurations)

consensus_share_#

The share of the aligned configurations’ total sum of squares that the consensus accounts for: one when all agree exactly.

Type:: float

explained_variance_ratio_#

The share of the consensus’s own variation on each of its axes.

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

n_dimensions_#

Columns in the widest configuration.

Type:: int

n_iter_#

Rounds run.

Type:: int

References

J.C. Gower, “Generalized Procrustes analysis”, Psychometrika, 40 (1975), 33-51.

J.M.F. ten Berge, “Orthogonal Procrustes rotation for two or more matrices”, Psychometrika, 42 (1977), 267-276.

G.M. Arnold and A.A. Williams, “The use of generalised Procrustes techniques in sensory analysis”, in Statistical Procedures in Food Research, Elsevier (1986).

I.N. Wakeling, M.M. Raats and H.J.H. MacFie, “A new significance test for consensus in generalized Procrustes analysis”, Journal of Sensory Studies, 7 (1992), 91-96.

Examples

>>> panel_configs = GPA.configurations_from_long(panel)
>>> gpa = GPA().fit(panel_configs)
>>> gpa.residuals_.sum().sort_values()      # who disagrees
>>> gpa.consensus_test(random_state=0).p_value

static configurations_from_long(frame, *, configuration='panelist_id', objects='product', variables='attribute', values='score')[source]#

Pivot a long table into one configuration per configuration value: objects by variables.

Values are averaged over everything else, such as replicates and sessions. The defaults are the column names of the sensory subpackage’s descriptive_long schema, so a validated panel goes straight in. A configuration keeps only the variables it used, as free-choice profiling needs.

Parameters:

frame (pd.DataFrame) – One row per value.
configuration (str, optional) – The columns naming the configuration (the assessor), the object (the product), the variable (the attribute) and holding the value (the score).
objects (str, optional) – The columns naming the configuration (the assessor), the object (the product), the variable (the attribute) and holding the value (the score).
variables (str, optional) – The columns naming the configuration (the assessor), the object (the product), the variable (the attribute) and holding the value (the score).
values (str, optional) – The columns naming the configuration (the assessor), the object (the product), the variable (the attribute) and holding the value (the score).

Returns:

Configuration name to a table of objects by variables, ready for fit().

Return type:

dict[str, pd.DataFrame]

Raises:

ValueError – If a named column is missing.

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

Align the configurations and find their consensus.

Parameters:

X (dict[str, pd.DataFrame]) – Configuration name to a table with one row per object and one column per dimension. Every table must have the same objects, in the same order, as its index; the columns may differ.
y (ignored) – Accepted for Pipeline compatibility.

Returns:

self, fitted.

Return type:

GPA

Raises:

TypeError – If X is not a dict of tables.
ValueError – If there are fewer than two configurations or three objects, the objects differ between configurations, a value is missing or not numeric, or a configuration puts every object at the same point.

Warns:

SpecificationWarning – If max_iter rounds end before the residual settles.

transform(X)[source]#

Place new objects in the consensus, through each configuration’s fitted alignment.

Parameters:: X (dict[str, pd.DataFrame]) – The same configurations, with the same columns, describing new objects (one or more rows each, the same objects in every configuration).
Returns:: One row per object: the mean of where each configuration puts it.
Return type:: pd.DataFrame
Raises:: ValueError – If a configuration is missing or unexpected, or its columns differ from the fit.

consensus_test(n_permutations=199, *, random_state=None)[source]#

Test whether the consensus is more than GPA would find among unrelated configurations.

The objects are shuffled independently within every configuration but the first, which destroys any agreement about which object is which while keeping each configuration’s own shape, and GPA is refitted. The p-value is the share of shuffles whose consensus accounts for at least as much of the total as the real one (Wakeling, Raats and MacFie, 1992).

Parameters:

n_permutations (int, optional) – Shuffles to run, default 199; the smallest p-value possible is one over this plus one.
random_state (int, np.random.Generator or None, optional) – Seed or generator for the shuffles; None draws fresh ones.

Returns:

consensus_share (the fitted one), null_shares (one per shuffle) and p_value.

Return type:

Bunch

Raises:

ValueError – If n_permutations is below 1.

map_plot(axis_x=1, axis_y=2, settings=None)[source]#

Draw the consensus, with a line from each object to where each configuration puts it.

Short spokes mean agreement; one configuration’s long spokes point at the assessor or instrument that sees the objects differently.

Parameters:

axis_x (int, optional) – The consensus axes to plot, counted from 1. Defaults 1 and 2.
axis_y (int, optional) – The consensus axes to plot, counted from 1. Defaults 1 and 2.
settings (dict, optional) – title (default “GPA consensus”), html_image_height (default 600), html_aspect_ratio_w_over_h (default 1) and template.

Return type:

go.Figure

Raises:

ValueError – If an axis is outside 1 to n_dimensions_.

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:

Notes

Each column is mean-centred internally but the blocks are not scaled; scale the columns yourself (for example with MCUVScaler) when the variables have different units.

References

Smilde, A. K., Kiers, H. A. L., Bijlsma, S., Rubingh, C. M. and van Erk, M. J. (2009). Matrix correlations for high-dimensional data: the modified RV-coefficient. Bioinformatics, 25(3), 401-405.

See also

rv_coefficient: The original RV coefficient.

Examples

>>> rv2_coefficient(X, Y)

Model selection#

Predictive cross-validation ($Q^2$, the 1-SE rule, van der Voet) placed beside tests of the latent structure itself; see Cross-Validation.

process_improve.multivariate.methods.compare_cv_criteria(X, Y, **kwargs)[source]#

Compare predictive and latent-structure validation criteria for a PLS model.

Function form of PLS.compare_cv_criteria(), which documents the settings and the returned fields: Q2, the 1-SE rule, van der Voet and CV-ANOVA beside the held-out score correlation, the out-of-sample slope ratio, a covariance permutation test, weight-subspace angles and Procrustes cross-validation, per component, with the number of components each rule recommends.

Parameters:

X (array-like of shape (n_samples, n_features)) – Predictor block, on its raw scale.
Y (array-like of shape (n_samples,) or (n_samples, n_targets)) – Response block.
**kwargs – The keyword arguments of PLS.compare_cv_criteria().

Returns:

See PLS.compare_cv_criteria().

Return type:

Examples

>>> from process_improve.multivariate import compare_cv_criteria
>>> result = compare_cv_criteria(X, y, max_components=6, random_state=0)
>>> result.recommendations["n_components"]

process_improve.multivariate.methods.pseudo_validation_set(X, Y, **kwargs)[source]#

Build a Procrustes pseudo-validation set for a PLS model.

Function form of PLS.pseudo_validation_set(), which documents the settings and the returned fields.

Parameters:

X (array-like of shape (n_samples, n_features)) – Predictor block, raw scale.
Y (array-like of shape (n_samples,) or (n_samples, n_targets)) – Response block.
**kwargs – The keyword arguments of PLS.pseudo_validation_set(); n_components is required.

Returns:

See PLS.pseudo_validation_set().

Return type:

Examples

>>> pv = pseudo_validation_set(X, y, n_components=2, random_state=0)
>>> pv.global_model.diagnose(pv.X_pv).spe

Containers#

class process_improve.multivariate.methods.BlockSet(blocks)[source]#

Bases: dict

A dict[str, pd.DataFrame] of equal-height blocks that can also be sliced by row (#193).

MBPCA.fit and MBPLS.fit take a plain dict[str, pd.DataFrame], which is convenient to build and impossible to resample: a dict has no notion of “row 7 of every block”. Any resampling or cross-validation pass needs exactly that. Resampler, for instance, asks its data only for len(x) and x[indices].

TPLS already has DataFrameDict for this, but it is hardwired to the Z/F/Y block names and to a nested dict[str, dict[str, DataFrame]] layout, so the multi-block models could not borrow it. BlockSet is the flat equivalent: a real dict subclass, so anything that already accepts the plain dict keeps working, plus row indexing.

Warning

len(blocks) is the number of rows, not the number of blocks. That is surprising for a dict, and it is deliberate: it is the convention DataFrameDict already set, and it is what the resampling code means by the length of a dataset. Use len(blocks.keys()) to count blocks.

Parameters:

blocks (dict[str, pd.DataFrame]) – One entry per block. Every block must be a DataFrame with the same number of rows; widths may differ.

Raises:

ValueError – If blocks is empty, or the blocks disagree on their row count.
TypeError – If any value is not a DataFrame.

Examples

>>> blocks = BlockSet({"a": df_a, "b": df_b})
>>> len(blocks)                                 # rows, not blocks
40
>>> blocks[[0, 1, 2]].keys()                    # a 3-row BlockSet
dict_keys(['a', 'b'])

class process_improve.multivariate.methods.DataFrameDict(datadict)[source]#

Bases: dict

Container for the partitionable (Z, F) and static (Y) data blocks used by TPLS.

Parameters:: datadict (dict[str, dict[str, pd.DataFrame]])

__init__(datadict)[source]#

Initialize a DataFrameDict to handle partitionable and static dataframes.

datadict: Dictionary with 3 keys, one for each block: Z, F and Y.: Each block is itself a dictionary of dataframes: dict[str, dict[str, pd.DataFrame]]

Parameters:: datadict (dict[str, dict[str, DataFrame]])

keys()[source]#

Return the keys of the DataFrameDict.

Return type:: KeysView[str]

Preprocessing#

class process_improve.multivariate.methods.MCUVScaler[source]#