Cross-Validation

Cross-validation is used for two purposes in multivariate analysis:

1. Component selection - choosing the right number of components, for both PCA and PLS.

2. Coefficient uncertainty - obtaining error bars for PLS beta coefficients.

Selecting the Number of Components

Choosing the right number of components is critical. Too few components underfit (miss important structure), too many overfit (model noise).

Element-wise Cross-Validation (PCA)

PCA.select_n_components() evaluates every component count from 1 to max_components, measures the Predicted Residual Error Sum of Squares (PRESS) of each, and recommends one. The default scheme is the element-wise k-fold (ekf) algorithm of Bro et al. (2008):

1. Split the individual cells of X into K folds, so that every cell is held out exactly once.

2. For each fold and each component count: mask that fold’s cells, impute them EM-style from a model fitted on the cells that remain, and add the squared error of those predictions to PRESS.

3. Recommend a component count from the PRESS curve using selection_rule.

Holding out individual cells, rather than whole rows, is what keeps a prediction independent of the value being predicted. When a whole row is held out and projected back through transform(), its own values flow into its prediction, so PRESS shrinks monotonically and reaches zero once the components equal the variables. That measures compression rather than prediction and cannot select a component count, which is why the cv_scheme="row_wise" scheme that did this was removed in 2.0.

Choosing among the schemes

Four schemes keep the prediction independent of the value predicted. They differ in what they hold out and what they cost.

cv_scheme	What is held out	Cost	Use it when
"ekf"	Scattered cells, imputed by EM from a model that never saw them	n_folds * n_repeats * max_components decompositions	The default. The only one that takes a block with missing cells.
"ek"	Nothing directly: a score comes from a model without the cell’s column, a loading from a model without its row	2 * n_folds decompositions	A number has to line up with Simca-P or with pcaMethods::Q2.
"sacv"	Nothing. Each residual is inflated by the leverage of the cell that produced it, approximating leave-one-cell-out	One decomposition	The block is large and "ekf" is too slow.
"gcv"	Nothing. One averaged leverage instead of one per cell	One decomposition	Comparing against FactoMineR or missMDA, whose default it is.

The three new schemes 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 hold nothing out in folds that could disagree, so they report no per-fold spread and selection_rule="1se" has nothing to work with.

Both leverage schemes are first-order approximations that degrade as the component count approaches the number of variables, because a column’s leverage approaches one and the divisor approaches zero with it. Read them well below that ceiling; FactoMineR defaults to five components for the same reason. At the ceiling itself no cell has a defined leave-one-out residual at all, and both schemes report NaN rather than a number there.

They also lean on the residual having something left in it. On the LDPE data of 54 rows and 19 variables, whose fit reaches 99.98% by eleven components, "ekf" and "ek" both turn over at two components, which is where Simca-P turns over on the same data. Neither "sacv" nor "gcv" turns over at all within the first five: their curves rise at every count, so the number they return is the largest one evaluated rather than an optimum. Any scheme that does that now says so, through a SpecificationWarning. Treat the warning as the answer: the criterion failed on that data, and a scheme that holds values out should be used instead.

from process_improve.multivariate.methods import PCA

# Pass the raw, unscaled X: with the default ``scale_inside_folds=True``
# the centring and scaling are fit inside each fold, so nothing about the
# held-out cells leaks into the model that predicts them.
result = PCA.select_n_components(
    X,
    max_components=10,
    cv=7,  # 7 element-folds
)

print(f"Recommended components: {result.n_components}")
print(result.press)  # PRESS per component count
print(result.q2)  # cross-validated R2 of X

The result is a Bunch with:

• n_components: recommended number of components

• press: PRESS for each number of components

• press_input_units: the same curve in the units of the matrix that was passed in, for comparing the prediction error against a known instrument error

• q2: cross-validated \(R^2_X\) per component count, on the same scale as the calibration r2_cumulative_ of a fitted model

• q2_per_variable: that same quantity for each column on its own

• per_fold_press, se_press and q2_se: the per-fold PRESS contributions and the standard error built from them, which is what the 1-SE rule needs

• press_ratio: the ratio PRESS_a / PRESS_{a-1}, for inspection

• cv_scores: per-fold scores (an alias of per_fold_press under ekf)

• cv_scheme and selection_rule: which scheme and rule were used

What PRESS is measured in

With scale_inside_folds=True each fold centres and scales the matrix before fitting it, and the error is measured in that same space. Every variable therefore contributes to press in proportion to how much of its own variation the model predicts, not in proportion to its units. This matters whenever the columns are on different scales: on the LDPE data used in the book, the Mw column carries 99.5% of the raw sum of squares, so a PRESS accumulated in the raw units would be Mw’s prediction error and almost nothing else.

q2 is that PRESS divided by what a null model would have got wrong on the same held-out cells, measured the same way. The null model predicts each held-out cell by the mean of the cells that were not held out, so \(Q^2 = 0\) is “no better than the column mean” and \(Q^2 = 1\) is exact prediction, the same reading as r2_cumulative_.

Two consequences are worth knowing. Re-expressing a column in different units (kilograms instead of grams, say) leaves the whole curve unchanged. And passing the raw block gives the same curve as passing a mean-centred, unit-variance block, so the recommendation in the paragraph above costs nothing.

When a single number in the original units is what you need, for instance to compare the prediction error against a known instrument error, read press_input_units instead. To see whether one column is carrying the pooled figure, read q2_per_variable, which splits q2 by variable.

scale_inside_folds=False is the opt-out for callers who have scaled their own block. There is then no in-fold scale, so press is in the units of whatever matrix was passed and the two PRESS fields coincide.

n_repeats runs the whole pass again with a fresh fold permutation. Each repeat still covers every cell exactly once; more repeats narrow se_press, which helps when the 1-SE rule sits on a borderline.

Selection Rules

selection_rule decides which count is recommended from the error curve. PCA defaults to "min"; PLS defaults to "1se".

Rule	Recommends
"min"	The component count with the lowest cross-validated error. This is the GlobalMin criterion that Bro et al. pair with ekf.
"1se"	The smallest count whose error is within one standard error of that minimum, so it is never less parsimonious than "min".
"q2_increment"	Keeps a component only while it lifts the cumulative \(Q^2\) by at least min_q2_increase (default 0.01), and stops at the first one that does not. A Wold’s-R-style heuristic: cheap, but the threshold is absolute and hand-tuned.
"randomization"	PLS only. Van der Voet’s (1994) permutation test: the smallest model whose predictive ability is statistically indistinguishable from the lowest-RMSECV one, at significance level alpha.

Note

The original Wold PRESS-ratio cutoff, the threshold argument of PCA.select_n_components, is deprecated: passing it emits a DeprecationWarning and the value is ignored. Use selection_rule="q2_increment", tuned with min_q2_increase, for a comparable preference for parsimony.

PLS Component Selection

PLS.select_n_components() cross-validates a PLS model and reports how it performs on unseen data, in contrast to the calibration statistics stored on a fitted model (rmse_, r2_cumulative_), which always improve as components are added.

from process_improve.multivariate.methods import PLS

# Raw, unscaled blocks: each training fold fits its own MCUVScaler and
# RMSECV is reported on the original Y scale.
result = PLS.select_n_components(X, Y, max_components=8, cv=5)

print(f"Recommended components: {result.n_components}")
print(result.rmsecv["total"])  # RMSECV per component count
print(result.r2y_validated["total"])  # Validated R2 of Y

Do not scale the blocks yourself before calling either selector. In-fold re-standardisation overwrites whatever scaling you applied, so two deliberately different choices (autoscale versus Pareto, say) become the same model and a comparison between them shows no difference. Both selectors emit a SpecificationWarning when they receive an X that is already centred and unit-variance scaled. If you must keep your own scaling, pass scale_inside_folds=False, and for PLS scale=False as well, since a PLS(scale=True) fold model autoscales its own rows; the scaling then leaks from the full dataset into every fold, and a warning says so.

The result is a Bunch with:

• n_components: recommended count, chosen by selection_rule (the 1-SE rule by default, not the lowest RMSECV; see Selection Rules)

• rmsecv: root-mean-square error of cross-validation, per Y variable and overall

• se_rmsecv / q2_se: the standard error of that curve, on the RMSECV and the \(Q^2\) scale respectively

• r2y_validated / r2x_validated: validated explained variance, per variable and overall. r2y_validated carries two overall columns: "total" pools the targets on the original Y scale, and "scaled_total" gives every target equal weight (see Comparing the fitted and the validated R2)

• press: overall Y prediction error sum of squares per component count

• cv_predictions: out-of-fold Y predictions at the recommended count

• selection_rule: the rule that produced n_components

The cv argument accepts an integer (K-fold) or any scikit-learn splitter object, such as KFold or LeaveOneOut. When it is an integer, the split is repeated n_repeats times (10 by default) with a fresh shuffle, which is what gives the 1-SE rule a usable standard error. A splitter object is used as-is and n_repeats is then ignored.

Comparing the fitted and the validated R2

A fitted model’s r2_y_cumulative_ is computed on the scaled Y, where every target carries equal weight. r2y_validated["total"] is computed on the original Y scale, where a target with a wide range carries more weight than a narrow one. Read side by side on targets of unequal spread, the two can differ by tens of percent, or disagree in sign, without the model having changed.

r2y_validated["scaled_total"] is the held-out number on the fitted model’s footing: every target weighted equally, which on autoscaled Y is the arithmetic mean of the per-target values. Compare fitted against held-out with those two, and keep "total" for the question it answers, which is how much of the Y variation in its own units the model predicts.

fitted = PLS(n_components=result.n_components).fit(X, Y)
print(fitted.r2_y_cumulative_)  # fitted, equal weight
print(result.r2y_validated["scaled_total"])  # held out, equal weight
print(result.r2y_validated["total"])  # held out, original Y scale

MBPLS.select_n_components returns the same two columns, on the same footing.

Validating the latent structure, not only the predictions

A PLS model can predict Y well while some of its components are not a reproducible property of the process, and a component can be real and stable without improving the prediction. \(Q^2\) only answers the first question. A model used for SPE and \(T^2\) monitoring, or inverted to find operating conditions, also depends on the second. compare_cv_criteria runs one cross-validation and reports both, component by component, with the number of components each criterion recommends.

from process_improve.multivariate import compare_cv_criteria, cv_criteria_plot

result = compare_cv_criteria(X, Y, max_components=6, cv=7, random_state=0)
print(result.recommendations[["n_components", "question"]])
print(result.table[["q2y", "r_cv", "r_cv_p", "slope_ratio", "cov_perm_p", "angle_subspace_deg"]])
cv_criteria_plot(result)

The same function is available as PLS.compare_cv_criteria, which fits with the class it is called on.

Each selection rule answers one question.

Rule	Question	Evidence
q2_max, q2_1se, van_der_voet	Does the model predict Y?	Held-out residuals of Y.
score_correlation	Does each component’s inner relation hold on new rows?	Correlation of the held-out \(t_a\) and \(u_a\) scores, against a null from permuted Y.
covariance_permutation	Is each component’s covariance larger than chance?	Largest singular value of the deflated \(\mathbf{X}_a^\top\mathbf{Y}_a\), against permuted rows of \(\mathbf{Y}_a\).
subspace_stability	Does the latent space survive a change of rows?	Largest principal angle between the fold and the full-data weights.

The predictive rules never return fewer than one component. The structural rules count the leading components that pass and return 0 when none does. Two components of nearly equal strength can swap from fold to fold: the span of the first is then unstable while the span of the pair is not, and subspace_stability counts the pair as two.

The monitoring question, whether the SPE limit is right for rows the model has not seen, is answered by a limit rather than by a component count. result.spe_limits holds, per component count, the SPE limit fitted to the training residuals (full_data, the one the model reports) and the limit fitted to the Procrustes pseudo-validation SPE (pseudo_validation). Training residuals are smaller than the residuals of new rows, so the first limit is too tight. On simulated data with two latent variables and 16 variables, the full-data limit at 95% flagged 17% of fresh rows with 30 training rows and 10% with 60, while the pseudo-validation limit flagged 4% to 5%. The pv_spe_limit_ratio column is the ratio of the two limits.

With missing values, a row’s SPE is summed over its observed cells only, so it is smaller than the SPE of the same row with nothing missing. The pseudo-validation limit is a limit for complete rows: each held-out row’s squared SPE is scaled by \(K / K_\text{observed}\) before the limit is fitted, and the SPE of a new row with gaps should be scaled the same way before the comparison. On the simulated data with 60 rows and 15% or 30% of the cells missing, the model’s own limit flagged 19% or 34% of fresh complete rows at a nominal 5%, and the refitted limit 3.5% or 2.8%.

How Q2 and the held-out correlation are related

The slope_ratio column, \(s_a\), is the slope of the held-out residual on the held-out score \(t_a\) relative to the slope the fold model fitted. It is exactly 1 on the training rows, and on held-out rows

\[\text{PRESS}_{a-1} - \text{PRESS}_a = (2 s_a - 1)\sum_k \|\tilde{\mathbf{c}}_{ak}\|^2\,\mathbf{t}_{ak}^\top\mathbf{t}_{ak}.\]

PRESS therefore falls, and \(Q^2\) rises, only when \(s_a > 1/2\). The held-out correlation is positive as soon as \(s_a > 0\). A component with \(0 < s_a < 1/2\) points the right way on new rows, but its fitted slope is more than twice too steep. Tested against its permutation null, the correlation does not keep such a component either: in simulations with a real but weak component, 4 of 83 components in this band reached \(p < 0.05\), no more than chance. The band marks a direction that is real but cannot yet be told from noise with this many rows.

Points to keep in mind when reading the table

• The held-out correlation is calibrated by permutation (n_cv_permutations). The normal approximation \(N(0, 1/N)\) is too narrow for a pooled cross-validated correlation, because every held-out score is a weighted sum of the other folds’ responses. On an X with two strong latent variables it fired about three times as often as its nominal 5% under a null Y.

• The angles are scaled by the delete-d jackknife, \(\tan\theta \to \sqrt{G-1}\,\mathrm{rms}(\tan\theta_k)\) for \(G\) folds. Raw fold angles shrink as the number of folds grows, because the fold models share most of their rows with the full-data model; the scaled angles agree between 7-fold and leave-one-out cross-validation.

• A stable direction is not the same as a direction related to Y. With a noise Y, the first weight vector \(\mathbf{w} \propto \mathbf{X}^\top\mathbf{y}\) is still drawn towards the dominant directions of X, so subspace_stability can pass where the two tests of a Y relationship do not.

• The SPE limits concern the X model only; they say nothing about whether a component predicts Y.

• Strong variation in X that Y does not see leaks into the first PLS weight vector in a finite sample, and a second component is then needed to cancel it. On simulated data with 60 rows, one Y-related latent variable and a stronger Y-orthogonal one, the first weight vector was tilted by a median of 9 degrees towards the Y-orthogonal direction, and two components predicted 4000 fresh rows better than one in 25 of 30 data sets. The predictive rules and the covariance test mostly chose two components, which is the right answer for prediction, though the process has one Y-related latent variable. This is the motivation for O-PLS, which moves the Y-orthogonal variation out of the predictive component.

• Missing values (NaN) are handled as the NIPALS fit handles them: each held-out row is scored from its observed cells the way NIPALS scores an incomplete training row, and a missing Y cell is left out of every sum. A row whose X is entirely missing, or a column that is, raises ValueError.

• CV-ANOVA (cv_anova_p, one response) is a monotone function of \(Q^2\) for fixed degrees of freedom, so it ranks models as \(Q^2\) does.

Checking the criteria on a process with a known structure

The simulated results quoted in this section come from LatentStructure, which draws data from a linear latent-variable process whose structure is set by the caller, and from fresh rows of the same process, which are the ground truth for any cross-validated estimate. The best that any linear model can do on those fresh rows is known in closed form.

from process_improve.multivariate import compare_cv_criteria
from process_improve.simulation import LatentStructure

# Two latent variables drive y, so a PLS model should need two components.
process = LatentStructure(x_sd=[3.0, 1.0], y_coefficients=[1.0, 1.0])
train = process.sample(60, random_state=0)
result = compare_cv_criteria(train.X, train.Y, max_components=5, random_state=0)
print(result.recommendations["n_components"])   # 2 from every rule
print(process.population_r2())                  # the Q2 a perfect linear model reaches

The script scripts/cv_criteria_recovery.py repeats this for every scenario quoted above, 30 data sets each, and prints how often each rule recovers the simulated count.

Procrustes pseudo-validation set

Ordinary cross-validation assesses the fold models, not the model fitted to all rows. pseudo_validation_set (Kucheryavskiy, Rodionova and Pomerantsev, 2023) converts the variation between the fold models into a data set, X_pv, that the full-data model can be applied to as if it were an independent test set.

from process_improve.multivariate import pseudo_validation_set

pv = pseudo_validation_set(X, y, n_components=2, random_state=0)
diagnostics = pv.global_model.diagnose(pv.X_pv)
alarm_rate = (diagnostics.spe > pv.global_model.spe_limit(0.95)).mean()

For every row, the full-data model gives X_pv the row’s fold-model SPE and its fold-model scores rescaled per component by the D ratio \(\mathbf{c}_{ka}^\top\mathbf{c}_a / \mathbf{c}_a^\top\mathbf{c}_a\). The SPE of X_pv is therefore the row-wise cross-validated X residual. That residual decreases almost monotonically with the number of components (see the note on row-wise cross-validation of X above), so the pseudo-validation set is used for alarm rates against the full-data limits, not to choose the number of components from its SPE. See Permutation nulls: is there anything here at all? for permutation tests of whether a model predicts anything at all.

PLS Beta Coefficient Error Bars

For PLS models, model.cross_validate() refits the model on data subsets and computes confidence intervals for the regression coefficients. This answers the question: “How reliable is each beta coefficient?”

Three resampling strategies are supported:

• Jackknife (cv="loo", default) - leave-one-out resampling. Uses the jackknife variance formula with t-distribution critical values.

• K-fold (cv=5) - K-fold cross-validation. Faster for large datasets.

• Bootstrap (n_bootstrap=200) - resample with replacement. Uses percentile confidence intervals.

from process_improve.multivariate.methods 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)

# Jackknife (leave-one-out) cross-validation
cv = pls.cross_validate(X_s, Y_s, cv="loo")

print(cv.significant)  # Which betas have CIs excluding zero
print(cv.beta_ci_lower)  # Lower 95% CI
print(cv.beta_ci_upper)  # Upper 95% CI
print(cv.q_squared)  # Cross-validated R² (Q²)
print(cv.rmse_cv)  # Cross-validated RMSE

The result is a Bunch with:

• beta_mean, beta_std: mean and standard error of betas across resamples

• beta_ci_lower, beta_ci_upper: confidence interval bounds

• significant: boolean mask - True where the CI excludes zero

• beta_samples: raw betas from every resample (n_resamples × K × M)

• y_hat_cv: out-of-fold Y predictions (jackknife / K-fold only)

• press: Prediction Error Sum of Squares

• rmse_cv: cross-validated RMSE per Y variable

• q_squared: cross-validated R² (Q²) per Y variable

See Projection to Latent Structures (PLS) for detailed documentation and additional examples.