fastPLS User Guide

macOS

Install Apple’s command-line developer tools before a source build:

xcode-select --install

macOS uses Apple Accelerate for CPU BLAS/LAPACK. Metal is detected automatically on supported Apple systems. To require Metal and stop installation if it cannot be compiled, set the following before the source installation:

Sys.setenv(FASTPLS_USE_METAL = "1")

In the benchmarks conducted on the Apple M3 system used for package testing, the Metal backend did not provide a meaningful computational-speed improvement over the macOS CPU backend. Performance may differ with the matrix dimensions, model configuration, and Apple hardware generation.

Ubuntu and Debian

Install the compiler toolchain and OpenBLAS development files before building the package:

sudo apt update
sudo apt install build-essential gfortran pkg-config libopenblas-dev

Then require OpenBLAS during the source build. This prevents an unnoticed fallback to the BLAS/LAPACK supplied by R:

Sys.setenv(FASTPLS_USE_OPENBLAS = "1")
install.packages("fastPLS", type = "source")

Ubuntu and Debian releases may provide different OpenBLAS versions, compilation options, threading implementations, and CPU kernels. Consequently, identical fastPLS calls can have substantially different runtimes across distributions or OpenBLAS installations without indicating a change in the statistical model. For reproducible timing, use a current build compiled for the target processor and record the version and selected kernel reported by fastPLS_blas(). Detecting the OpenBLAS family alone does not establish that an architecture-appropriate kernel is active.

Fedora

Install the corresponding development packages, then use the same R command shown above:

sudo dnf install gcc gcc-c++ gcc-gfortran make pkgconf-pkg-config \
    openblas-devel

Windows x86-64

Install the Rtools release that matches the installed R version. A source build can use R’s BLAS/LAPACK, but OpenBLAS is recommended for the Linux and Windows performance routes evaluated with fastPLS. One option is to install an x86-64 OpenBLAS archive from an MSYS2 UCRT64 terminal:

pacman -S --needed mingw-w64-ucrt-x86_64-openblas

Point the package configuration to the matching prefix and require OpenBLAS:

Sys.setenv(
    FASTPLS_USE_OPENBLAS = "1",
    OPENBLAS_ROOT = "C:/msys64/ucrt64"
)
remotes::install_github(
    "tkcaccia/fastPLS",
    upgrade = "never",
    force = TRUE,
    build_vignettes = TRUE
)

The OpenBLAS library must match the target architecture. In particular, a Windows ARM64 build must not use an x86-64 Rtools or MSYS2 archive. If a compatible ARM64 OpenBLAS development package is unavailable, leave automatic detection enabled and allow the build to use R’s BLAS/LAPACK instead.

Windows OpenBLAS distributions can also differ in version, compiler options, threading implementation, and processor-specific kernel. These differences can materially affect runtime even when the fastPLS version and model settings are unchanged. For reproducible benchmarks, record the OpenBLAS version, selected kernel, thread count, and library path reported by fastPLS_blas().

If OpenBLAS is installed in a nonstandard location on Linux or Windows, set OPENBLAS_ROOT to its installation prefix. Set FASTPLS_USE_OPENBLAS = "0" only when an R-supplied BLAS build is explicitly desired.

Verify the Compiled Libraries

After installation, restart R and verify the selected CPU library and optional accelerators:

library(fastPLS)

fastPLS_blas()
has_cuda()
has_metal()

fastPLS_blas() reports the backend together with the OpenBLAS version, configuration, selected CPU core, parallel runtime, active thread count, and resolved library path when available. Linux and Windows benchmarks should proceed only when fastPLS_blas()$backend is "OpenBLAS" and the expected version and core are active. Use fastPLS_blas(details = FALSE) when only the former scalar backend name is required. The campaign tools in fastPLS-extra validate the detailed report before running fastPLS timing stages. The separate installation vignette contains CUDA toolkit requirements, environment variables, and troubleshooting for architecture or linker errors.

Choosing a Starting Model

The public classification choices are "argmax" and "lda"; regression ignores classifier. LDA uses the fixed scale-normalized Cholesky fallback sequence described above and does not expose a ridge-tuning argument.

Fit And Predict A Classifier

The method argument selects the PLS algorithm, backend selects the implementation, and classifier selects the classification head. The default classification head is classifier = "lda", which fits linear discriminant analysis in the retained PLS score space. Use classifier = "argmax" explicitly to predict the class with the largest PLS-DA response score.

fit_cls <- pls(
    Xtrain,
    Ytrain_cls,
    Xtest,
    Ytest_cls,
    ncomp = 1:2,
    fit = TRUE,
    return_variance = FALSE,
    seed = 101
)

fit_cls$accuracy
#> ncomp=1 ncomp=2 
#>       1       1

Classification models can also be fitted once and predicted later. predict() returns only the predicted class by default and can optionally return ranked classes by setting top to a positive integer. The ordinary prediction in Ypred is always the rank-1 class. When top > 1, Ypred_top contains the ordered candidate labels for each sample: rank1 is the predicted class, rank2 is the next most likely class, and so on. If available, Ypred_top_score contains the corresponding class scores used to create that ranking. For both float64 and float32 models, the default raw_scores = FALSE evaluates ranked predictions in bounded row blocks and retains only the requested ranks. Set raw_scores = TRUE only when the complete class-score arrays are needed, because those arrays can be substantially larger.

fit_cls_train_only <- pls(
    Xtrain,
    Ytrain_cls,
    ncomp = 1:2,
    classifier = "lda",
    fit = TRUE,
    return_variance = FALSE,
    seed = 101
)

pred_cls_later <- predict(
    fit_cls_train_only,
    Xtest,
    Ytest = Ytest_cls,
    top = 2,
    raw_scores = TRUE
)

pred_cls_later$accuracy
#> ncomp=1 ncomp=2 
#>       1       1
pred_cls_later$metrics$metrics
#>          n accuracy no_information_rate lift_accuracy
#> ncomp=1 30        1           0.3333333             3
#> ncomp=2 30        1           0.3333333             3
#>         balanced_accuracy macro_precision macro_recall macro_f1
#> ncomp=1                 1               1            1        1
#> ncomp=2                 1               1            1        1
#>         kappa
#> ncomp=1     1
#> ncomp=2     1
head(pred_cls_later$Ypred_top[["ncomp=2"]])
#>      rank1        rank2       
#> [1,] "virginica"  "versicolor"
#> [2,] "virginica"  "versicolor"
#> [3,] "setosa"     "versicolor"
#> [4,] "versicolor" "virginica" 
#> [5,] "versicolor" "virginica" 
#> [6,] "versicolor" "virginica"

Evaluate Classification Predictions

Use evaluate() to summarize predicted class labels, ranked labels, class-score matrices, or a complete object returned by predict(). The task and requested rank are inferred from these inputs. For classification, lift_accuracy is accuracy divided by the no_information_rate, the accuracy obtained by always predicting the most frequent observed class. Values above one therefore improve on this simple majority-class baseline. For classification, the complete output includes global metrics, per-class metrics, and the confusion matrix.

eval_cls_path <- evaluate(
    observed = Ytest_cls,
    predicted = pred_cls_later
)

eval_cls <- eval_cls_path$by_component[["ncomp=2"]]
eval_cls
#> $task
#> [1] "classification"
#> 
#> $metrics
#>    n accuracy no_information_rate lift_accuracy balanced_accuracy
#> 1 30        1           0.3333333             3                 1
#>   macro_precision macro_recall macro_f1 kappa
#> 1               1            1        1     1
#> 
#> $metric_definitions
#> $metric_definitions$accuracy
#> [1] "Proportion of observed labels predicted correctly."
#> 
#> $metric_definitions$balanced_accuracy
#> [1] "Unweighted mean of class-specific recalls."
#> 
#> 
#> $per_class
#>        class support precision recall f1
#> 1     setosa      10         1      1  1
#> 2 versicolor      10         1      1  1
#> 3  virginica      10         1      1  1
#> 
#> $confusion
#>             observed
#> predicted    setosa versicolor virginica
#>   setosa         10          0         0
#>   versicolor      0         10         0
#>   virginica       0          0        10
#> 
#> $topk
#>   k accuracy
#> 1 1        1
#> 2 2        1

The same complete evaluation is also stored automatically in a fitted PLS object when observed test responses are supplied. Results are grouped first by data role and then by component count.

fit_cls$metrics$test[["ncomp=2"]]$metrics
#>    n accuracy no_information_rate lift_accuracy balanced_accuracy
#> 1 30        1           0.3333333             3                 1
#>   macro_precision macro_recall macro_f1 kappa
#> 1               1            1        1     1

Rows of the confusion matrix are predicted labels and columns are observed labels. The confusion matrix is returned as an ordinary R table.

eval_cls$confusion
#>             observed
#> predicted    setosa versicolor virginica
#>   setosa         10          0         0
#>   versicolor      0         10         0
#>   virginica       0          0        10

When class-score matrices are available, evaluate() can also report top-k accuracy. Top-k accuracy asks whether the true class appears anywhere among the first k ranked labels. Thus top-1 accuracy is ordinary classification accuracy, while top-5 accuracy gives credit when the correct class appears among the five highest-scoring alternatives.

score_last <- pred_cls_later$LDA_scores[
    , , dim(pred_cls_later$LDA_scores)[3L]
]

evaluate(
    observed = Ytest_cls,
    predicted = score_last
)
#> $task
#> [1] "classification"
#> 
#> $metrics
#>    n accuracy no_information_rate lift_accuracy balanced_accuracy
#> 1 30        1           0.3333333             3                 1
#>   macro_precision macro_recall macro_f1 kappa
#> 1               1            1        1     1
#> 
#> $metric_definitions
#> $metric_definitions$accuracy
#> [1] "Proportion of observed labels predicted correctly."
#> 
#> $metric_definitions$balanced_accuracy
#> [1] "Unweighted mean of class-specific recalls."
#> 
#> 
#> $per_class
#>        class support precision recall f1
#> 1     setosa      10         1      1  1
#> 2 versicolor      10         1      1  1
#> 3  virginica      10         1      1  1
#> 
#> $confusion
#>             observed
#> predicted    setosa versicolor virginica
#>   setosa         10          0         0
#>   versicolor      0         10         0
#>   virginica       0          0        10
#> 
#> $topk
#>   k accuracy
#> 1 1        1
#> 2 2        1
#> 3 3        1

Classification Heads

The same pls() interface exposes two classification-specific heads for factor responses: argmax PLS-DA and latent-space LDA. They are decoders applied after the PLS model has produced class-response scores or latent scores. For regression, leave classifier at its default; numeric responses are predicted directly as continuous values.

fit_cls_plssvd <- pls(
    Xtrain,
    Ytrain_cls,
    Xtest,
    Ytest_cls,
    ncomp = 1:2,
    method = "plssvd",
    seed = 100
)

head(fit_cls_plssvd$Ypred)
#>      ncomp=1    ncomp=2
#> 1  virginica  virginica
#> 2  virginica  virginica
#> 3     setosa     setosa
#> 4 versicolor versicolor
#> 5 versicolor versicolor
#> 6 versicolor versicolor

evaluate(
    observed = Ytest_cls,
    predicted = fit_cls_plssvd$Ypred[["ncomp=2"]]
)$confusion
#>             observed
#> predicted    setosa versicolor virginica
#>   setosa         10          0         0
#>   versicolor      0         10         0
#>   virginica       0          0        10

For PLS-DA with an LDA prediction head, use classifier = "lda". On systems with GPU support, the backend is selected through backend = "cuda" or backend = "metal" where available:

fit_cls_lda_gpu <- pls(
    Xtrain,
    Ytrain_cls,
    Xtest,
    Ytest_cls,
    ncomp = 1:2,
    method = "plssvd",
    backend = "cuda",
    classifier = "lda"
)

An unavailable CUDA or Metal selection raises an error; fastPLS never changes the request silently to CPU. To run on the CPU, select it explicitly:

fit_cls_lda_cpu <- pls(
    Xtrain,
    Ytrain_cls,
    Xtest,
    Ytest_cls,
    ncomp = 1:2,
    method = "plssvd",
    backend = "cpu",
    seed = 100,
    classifier = "lda"
)

head(fit_cls_lda_cpu$Ypred)
#>      ncomp=1    ncomp=2
#> 1  virginica  virginica
#> 2  virginica  virginica
#> 3     setosa     setosa
#> 4 versicolor versicolor
#> 5 versicolor versicolor
#> 6 versicolor versicolor

Kernel PLS For Classification

Kernel PLS changes the representation of the samples before the inner PLS fit. The linear kernel is equivalent to an ordinary inner-product representation and is useful as a fast baseline. The rbf kernel uses a radial-basis similarity; the poly kernel uses polynomial interactions among features.

kernel_fits <- lapply(c("linear", "rbf", "poly"), function(k) {
    pls(
    Xtrain, Ytrain_cls, Xtest, Ytest_cls,
    ncomp = 1:2,
    method = "kernelpls",
    kernel = k,
    degree = 2,
    seed = 102
    )
})
names(kernel_fits) <- c("linear", "rbf", "poly")

kernel_accuracy <- vapply(kernel_fits, function(fit) {
    mean(fit$Ypred[["ncomp=2"]] == Ytest_cls)
}, numeric(1))
kernel_accuracy
#>    linear       rbf      poly 
#> 1.0000000 0.9333333 0.9666667

Classification Score Plots

plot() can visualize stored PLS score maps. Refit with fit = TRUE to store training scores, or predict with proj = TRUE to store test scores. Ellipses can be ordinary confidence ellipses or Hotelling T2 ellipses.

plot(
    fit_cls,
    groups = Ytrain_cls,
    ellipse = TRUE,
    ellipse.type = "confidence"
)

The following example fits OPLS to the three-class iris split defined at the beginning of this section, with two orthogonal components (north = 2). The same plot() interface is available for PLS-SVD, the SIMPLS-family estimator, OPLS, and kernel PLS.

opls_three <- pls(
    Xtrain,
    Ytrain_cls,
    Xtest,
    Ytest_cls,
    ncomp = 1:2,
    method = "opls",
    fit = TRUE,
    proj = TRUE,
    north = 2L,
    seed = 201
)

old_par <- par(mfrow = c(1, 2), mar = c(4.2, 4.2, 2.2, 0.8))
plot(
    opls_three,
    score.set = "train",
    groups = Ytrain_cls,
    ellipse = TRUE,
    ellipse.type = "hotelling",
    main = "OPLS training, north = 2"
)
plot(
    opls_three,
    score.set = "test",
    groups = opls_three$Ypred[["ncomp=2"]],
    xlim = c(-0.16, 0.42),
    main = "OPLS test prediction"
)

par(old_par)

Univariate Regression

For a univariate task, supply the response as a numeric vector. The simplest workflow is to provide both the training data and an independent test set directly to pls(). The fitted object then contains test-set predictions and, when Ytest is supplied, predictive metrics.

fit_reg <- pls(
    Xreg_train,
    Ytrain_reg,
    Xreg_test,
    Ytest_reg,
    ncomp = 1:3,
    fit = TRUE,
    return_variance = FALSE
)

fit_reg$Q2Y
#>   ncomp=1   ncomp=2   ncomp=3 
#> 0.7564006 0.7593961 0.8634626

Float64 and Float32

Standard R numeric matrices use float64, whereas a float::float32 object requests float32 execution. Float64 stores each value in eight bytes and retains approximately 15–16 decimal digits of precision. Float32 stores each value in four bytes and retains approximately seven decimal digits. Float32 can therefore halve the representation size of a matrix and reduce data-transfer and memory-bandwidth costs, but it also introduces more rounding error.

Reduced representation size does not guarantee a faster fit or lower peak process memory. Temporary workspaces, decomposition costs, host-device transfers, and backend-specific kernels can dominate the calculation. Runtime, incremental host memory, and device memory may therefore increase or decrease when float32 replaces float64.

The lower precision can affect nearly tied latent directions, selected component counts, predictions, or classifications close to a decision boundary. For a new scientific workflow, compare float32 and float64 using the same split, seed, component grid, and model settings, and examine prediction differences as well as the final performance metric. Float64 remains the baseline precision when small numerical differences could alter the scientific conclusion.

CUDA keeps PLS-SVD, the SIMPLS-family estimator, OPLS, nonlinear kernel construction, prediction, and LDA on the selected GPU. Metal accepts float32 only and uses the fixed CPU/Metal operation split described above; CUDA supports float32 and float64. Nonlinear kernel PLS requires an explicit n by n Gram matrix, so the package checks the estimated live device or unified-memory storage before fitting.

The package warns once for measured-risk regimes. These include precision-sensitive SIMPLS-family and linear kernel-PLS classification, nonlinear kernels, and multivariate regression with at least 10,000 response columns and 50 components. In the latter regime, runtime, memory, and numerical behavior must be checked against float64 before the float32 result is used for scientific interpretation. Unsupported combinations stop before allocation and are never silently promoted to float64.

On Windows, the standard R toolchain does not expose the single-precision BLAS/LAPACK symbols used by the Unix-like compiled kernels. The Windows CPU route therefore combines float-package rSVD and Cholesky operations with portable C++ float kernels. It supports PLS-SVD, the SIMPLS-family estimator, OPLS, linear and nonlinear kernel PLS, argmax, and latent-space LDA when backend = "cpu" and rSVD is selected automatically. Model factors, LDA buffers, kernel matrices, and predictions remain float32, but this portable route can be slower than the native Unix-like implementation. Windows float32 accelerator requests stop with an unsupported-combination error rather than silently converting data to float64.

Xreg32 <- float::fl(as.matrix(Xreg_train))
Yreg32 <- float::fl(matrix(Ytrain_reg, ncol = 1))
fit_reg32 <- pls(
    Xreg32,
    Yreg32,
    float::fl(as.matrix(Xreg_test)),
    float::fl(matrix(Ytest_reg, ncol = 1)),
    ncomp = 1:2
)
fit_reg32$Q2Y
#>   ncomp=1   ncomp=2 
#> 0.7564006 0.7593962

For standard double-precision regression, Ypred is a numeric prediction array with one slice for each requested number of components. For float32 input, Ypred is kept as a named list of float::float32 prediction matrices. In the standard example below, the last array slice corresponds to the largest requested component count (ncomp = 3). The metric vectors are named by component count, for example fit_reg$Q2Y["ncomp=3"].

reg_component <- dim(fit_reg$Ypred)[3L]
pred_mpg <- fit_reg$Ypred[, , reg_component]
plot(
    Ytest_reg,
    pred_mpg,
    pch = 21,
    bg = "#4E79A7",
    col = "black",
    xlab = "Observed mpg",
    ylab = "Predicted mpg",
    main = "Regression: observed vs predicted"
)
abline(0, 1, col = "#D55E00", lwd = 2)

The alternative workflow is to fit the model once without a test set, then call predict() later. This is useful when the same model must be applied to several independent datasets. predict() automatically applies the centering/scaling stored in the fitted object. When Ytest is supplied, it also passes the predictions to evaluate() and stores the complete evaluation under metrics.

fit_reg_train_only <- pls(
    Xreg_train,
    Ytrain_reg,
    ncomp = 1:3,
    fit = TRUE,
    return_variance = FALSE
)

pred_reg_later <- predict(
    fit_reg_train_only,
    Xreg_test,
    Ytest = Ytest_reg,
    proj = TRUE
)

pred_reg_later$Q2Y
#>   ncomp=1   ncomp=2   ncomp=3 
#> 0.7564006 0.7593961 0.8634626
pred_reg_later$metrics$metrics
#>             n        R2        Q2     RMSD     RMSE      MAE
#> component=1 8 0.7505150 0.7564006 3.066116 3.066116 2.739049
#> component=2 8 0.7535829 0.7593961 3.047205 3.047205 2.765666
#> component=3 8 0.8601637 0.8634626 2.295494 2.295494 1.929671
#>                  bias MRE_percent MAPE_percent      RPD Pearson_r
#> component=1 0.4938232    16.70827     14.92070 2.140295 0.8700688
#> component=2 0.6250298    17.87218     15.26995 2.153578 0.8743339
#> component=3 0.6784321    12.48956     10.76821 2.858815 0.9340810
#>             Spearman_r
#> component=1  0.8072875
#> component=2  0.8072875
#> component=3  0.8072875
head(pred_reg_later$Ttest)
#>             [,1]        [,2]       [,3]
#> [1,]  0.09841041 -0.01046739 -0.1900955
#> [2,] -0.10644825  0.19820127  0.1563451
#> [3,]  0.03394295  0.21190609 -0.2208999
#> [4,] -0.36597585  0.27972780 -0.2105406
#> [5,]  0.27550586  0.09182286  0.3494802
#> [6,] -0.25506247  0.31745909  0.2627156

Evaluate Regression Predictions

For numeric regression, evaluate() reports R2, Q2, RMSD/RMSE, MAE, bias, median relative error percentage, RPD, and correlations. If ytrain is supplied, Q2 is calculated relative to the training-set response mean, which is the preferred setting for independent test-set evaluation. For a response matrix, each response column is centered on its corresponding training mean before the denominator sums are aggregated. If ytrain is omitted, Q2 is returned as NA; it is not silently replaced by an R2 calculation.

For wide multivariate responses, pls(), pls.single.cv(), and pls.double.cv() return aggregate evaluation metrics by default. Set bycol = TRUE to also calculate the full response-wise metric table; this is useful for inspecting individual spectral bins, but can be expensive for NMR data with many response columns.

eval_reg <- evaluate(
    observed = Ytest_reg,
    predicted = pred_mpg,
    ytrain = Ytrain_reg
)

eval_reg$task
#> [1] "regression"
eval_reg$metrics
#>   n        R2        Q2     RMSD     RMSE      MAE      bias
#> 1 8 0.8601637 0.8634626 2.295494 2.295494 1.929671 0.6784321
#>   MRE_percent MAPE_percent      RPD Pearson_r Spearman_r
#> 1    12.48956     10.76821 2.858815  0.934081  0.8072875
lapply(eval_reg$metric_definitions, strwrap, width = 56)
#> $R2
#> [1] "Observed-set R2; each response is centered on its"   
#> [2] "observed mean before sums of squares are aggregated."
#> 
#> $Q2
#> [1] "Independent-test Q2; each response is centered on its"
#> [2] "training-response mean before sums of squares are"    
#> [3] "aggregated."
eval_reg$per_response
#>   response n        R2        Q2     RMSD     RMSE      MAE
#> 1       Y1 8 0.8601637 0.8634626 2.295494 2.295494 1.929671
#>        bias MRE_percent MAPE_percent      RPD Pearson_r Spearman_r
#> 1 0.6784321    12.48956     10.76821 2.858815  0.934081  0.8072875

Multivariate Regression

For multivariate regression, supply a numeric response matrix with one column per outcome. The following example predicts mpg, qsec, and drat simultaneously from a common predictor matrix. bycol = TRUE requests both the aggregate multivariate metrics and a response-wise summary.

Xmulti <- as.matrix(mtcars[, c("disp", "hp", "wt", "gear", "carb")])
Ymulti <- as.matrix(mtcars[, c("mpg", "qsec", "drat")])
Xmulti_train <- Xmulti[-reg_test_id, , drop = FALSE]
Xmulti_test <- Xmulti[reg_test_id, , drop = FALSE]
Ymulti_train <- Ymulti[-reg_test_id, , drop = FALSE]
Ymulti_test <- Ymulti[reg_test_id, , drop = FALSE]

fit_multi <- pls(
    Xtrain = Xmulti_train,
    Ytrain = Ymulti_train,
    Xtest = Xmulti_test,
    Ytest = Ymulti_test,
    ncomp = 1:3,
    method = "plssvd",
    bycol = TRUE,
    return_variance = FALSE,
    seed = 102
)

fit_multi$metrics$test[["ncomp=3"]]$metrics
#>    n        R2        Q2     RMSD     RMSE      MAE      bias
#> 1 24 0.7873104 0.7938265 1.663514 1.663514 1.145866 0.2123291
#>   MRE_percent MAPE_percent      RPD Pearson_r Spearman_r
#> 1    5.694746     8.996205 4.965097 0.9791838  0.8773381
fit_multi$metrics$test[["ncomp=3"]]$per_response
#>   response n        R2        Q2      RMSD      RMSE      MAE
#> 1      mpg 8 0.8022678 0.8069325 2.7296376 2.7296376 2.399265
#> 2     qsec 8 0.2848075 0.4708414 0.7892148 0.7892148 0.688538
#> 3     drat 8 0.5247638 0.5410309 0.4775554 0.4775554 0.349794
#>           bias MRE_percent MAPE_percent      RPD Pearson_r
#> 1  0.610220685   14.166000    13.262266 2.404126 0.9012015
#> 2  0.028614187    4.359165     3.746797 1.264109 0.5951443
#> 3 -0.001847712    5.265462     9.979553 1.550748 0.7560230
#>   Spearman_r
#> 1  0.7470422
#> 2  0.6904762
#> 3  0.6666667

Component selection uses the same interface. With a numeric response matrix, the default criterion is aggregate held-out RMSD across all response columns.

cv_multi <- pls.single.cv(
    Xdata = Xmulti_train,
    Ydata = Ymulti_train,
    ncomp = 1:3,
    kfold = 3,
    method = "plssvd",
    fit = FALSE,
    bycol = TRUE,
    return_splits = TRUE,
    seed = 102
)

cv_multi$best_ncomp
#> [1] 1
cv_multi$best_metric_value
#> [1] 2.236604
head(cv_multi$split_index)
#>   fold_1     fold_2     fold_3    
#> 1 "training" "training" "test"    
#> 2 "training" "training" "test"    
#> 3 "training" "test"     "training"
#> 4 "test"     "training" "training"
#> 5 "test"     "training" "training"
#> 6 "training" "test"     "training"

OPLS For Regression

OPLS is accessed through the same pls() function with method = "opls". It is often used to separate predictive variation from response-orthogonal variation. ncomp counts predictive components, separately from the orthogonal components removed by north. Each removed direction reduces the available predictor rank. For example, with eight independent predictors and one removed direction, at most seven predictive components remain. A request exceeding the remaining rank raises an error; the same restriction applies within each CV training fold.

fit_opls <- pls(
    Xreg_train,
    Ytrain_reg,
    Xreg_test,
    Ytest_reg,
    ncomp = 1:2,
    method = "opls",
    seed = 101
)

class(fit_opls)
#> [1] "fastPLSOpls" "fastPLS"
fit_opls$Q2Y
#>   ncomp=1   ncomp=2 
#> 0.7593961 0.8634626

Ordinary k-fold CV

In ordinary k-fold CV, samples are split directly into kfold folds. The example below validates a fixed two-component SIMPLS-family classifier with five folds.

cv_kfold <- pls.single.cv(
    Xdata = Xtrain,
    Ydata = Ytrain_cls,
    ncomp = 2,
    kfold = 5,
    return_splits = TRUE,
    seed = 103
)

cv_kfold$metrics$cross_validated[["ncomp=2"]]$metrics
#>     n  accuracy no_information_rate lift_accuracy balanced_accuracy
#> 1 120 0.8333333           0.3333333           2.5         0.8333333
#>   macro_precision macro_recall  macro_f1 kappa
#> 1       0.8422619    0.8333333 0.8320269  0.75
head(cv_kfold$split_index)
#>   fold_1     fold_2     fold_3     fold_4     fold_5    
#> 1 "training" "training" "training" "test"     "training"
#> 2 "training" "test"     "training" "training" "training"
#> 3 "training" "training" "training" "test"     "training"
#> 4 "test"     "training" "training" "training" "training"
#> 5 "training" "training" "test"     "training" "training"
#> 6 "training" "test"     "training" "training" "training"

metrics$cross_validated contains the complete evaluate() output for each tested component count, while selection_metrics is the compact internal table used to select the best setting. When fit = TRUE, metrics$fitted evaluates the corresponding model fitted to the full dataset. With return_splits = TRUE, split_index has one row per input sample and one column per fold; entries identify whether that sample was used for training or testing. Nested CV uses the same optional field for all outer and inner splits; outer_test marks samples that are unavailable to a given inner CV cycle.

Grouped k-fold CV with constrain

The constrain argument controls grouped splitting. It is a vector with one entry per sample; samples with the same value are assigned to the same fold. In practice, this prevents leakage when multiple rows come from the same patient, subject, batch, or technical replicate. For example, if two spectra come from the same patient, giving them the same patient identifier in constrain ensures that both spectra are placed either in the training set or in the test set, never one in each.

patient_id <- rep(
    seq_len(ceiling(nrow(Xtrain) / 2)), each = 2
)[seq_len(nrow(Xtrain))]

cv_grouped <- pls.single.cv(
    Xdata = Xtrain,
    Ydata = Ytrain_cls,
    constrain = patient_id,
    ncomp = 2,
    kfold = 4,
    return_splits = TRUE,
    seed = 104
)

c(n_folds = length(unique(cv_grouped$fold)),
    n_patient_groups = length(unique(patient_id)))
#>          n_folds n_patient_groups 
#>                4               60

head(data.frame(
    sample_index = seq_len(nrow(Xtrain)),
    patient_id = patient_id,
    cv_grouped$split_index,
    check.names = FALSE
), 8)
#>   sample_index patient_id   fold_1   fold_2   fold_3   fold_4
#> 1            1          1 training training training     test
#> 2            2          1 training training training     test
#> 3            3          2 training training     test training
#> 4            4          2 training training     test training
#> 5            5          3 training     test training training
#> 6            6          3 training     test training training
#> 7            7          4 training     test training training
#> 8            8          4 training     test training training

patient_rows <- split(seq_along(patient_id), patient_id)
patients_kept_together <- vapply(
    seq_len(ncol(cv_grouped$split_index)),
    function(fold) {
        all(vapply(patient_rows, function(rows) {
            length(unique(cv_grouped$split_index[rows, fold])) == 1L
        }, logical(1L)))
    },
    logical(1L)
)

stopifnot(all(patients_kept_together))
data.frame(
    fold = colnames(cv_grouped$split_index),
    patients_kept_together = patients_kept_together
)
#>     fold patients_kept_together
#> 1 fold_1                   TRUE
#> 2 fold_2                   TRUE
#> 3 fold_3                   TRUE
#> 4 fold_4                   TRUE

The displayed rows show that the two samples from each patient receive the same training or test label within a fold. The executable stopifnot() check applies this rule to every patient and every fold, so the vignette build fails if grouped splitting ever places samples from one patient on opposite sides of a split.

Leave-one-out CV and leave-one-group-out CV

Leave-one-out CV is requested with kfold = "loocv". When constrain is not supplied, each sample is held out once. When constrain is supplied, LOOCV becomes leave-one-constraint-group-out CV, so a whole patient, subject, batch, or replicate group is held out together at each iteration. Numeric kfold values greater than or equal to the number of constraint groups are also treated as leave-one-group-out CV.

cv_loocv <- pls.single.cv(
    Xdata = Xtrain,
    Ydata = Ytrain_cls,
    constrain = patient_id,
    ncomp = 1,
    kfold = "loocv",
    seed = 105
)

c(n_loocv_folds = length(unique(cv_loocv$fold)),
    n_patient_groups = length(unique(patient_id)))
#>    n_loocv_folds n_patient_groups 
#>               60               60

Component and hyperparameter optimization CV

pls.single.cv() repeats the same CV splitting strategy over several candidate component counts and returns the best value according to the predictive metric: accuracy for classification and RMSD for regression. Fold-aware Q2Y, full-data fitted R2Y, balanced accuracy, and the other task-specific metrics documented below can instead be requested explicitly through selection. Predictive arguments can also be supplied as vectors to tune a compact grid. Fold-control arguments such as kfold remain single settings for the whole CV run.

For LDA classification, a rare class may be absent from an inner training fold. fastPLS fits LDA to the classes represented in that fold and maps the predictions back to the original factor levels. Held-out observations from an absent class are retained in the accuracy and balanced-accuracy calculations, although that class cannot be predicted by a model that did not observe it. When otherwise identical argmax and LDA configurations are requested together, fastPLS calculates the PLS component path and fold projection once. The two classification heads retain separate predictions, metric paths, and selected component counts. For sufficiently tall SIMPLS classification problems, full-data predictor and class moments are also computed once. Each training fold obtains its moments by subtracting the held-out contribution before applying fold-specific centering and scaling; the PLS model, LDA model, and predictions are still fitted independently within every fold. This avoids repeatedly materializing large training-score matrices without leaking held-out information into model fitting. The same sufficient statistics can assemble eligible PLS-SVD folds through the fold-specific predictor Gram and class cross-product, without materializing the training-score matrix. On the Metal route, small reduced products are evaluated with Accelerate while the large sample-matrix products assigned to Metal remain on the GPU.

Some grouped or highly imbalanced folds have lower effective rank than the requested component path. In that case, fastPLS fits only the estimable score prefix and repeats its prediction and metric for larger requested counts. The requested path is preserved, while effective_ncomp records the prefix used in each fold. For LDA, a fold with no estimable PLS direction uses its empirical training-class priors and finite log-prior scores; a training fold containing one class predicts that class. These fallbacks are reported through the fold status field rather than surfacing as a matrix-dimension error.

The same rule applies to a direct regression fit. If fewer directions are estimable than requested, effective_ncomp records the usable prefix and later paths repeat the last estimable prediction and coefficient matrix. If the training response is constant, no direction is estimable: predictions equal the training-response mean, coefficient paths are zero, and response-variance metrics with a zero denominator are NA.

cv_opt <- pls.single.cv(
    Xdata = Xreg_train,
    Ydata = Ytrain_reg,
    ncomp = 1:3,
    kfold = 5
)

cv_opt$best_ncomp
#> [1] 1

Interpreting R2Y, Q2Y, and RMSD

For regression models, R2Y, Q2Y, and RMSD answer different questions and should not usually be identical. In pls(), R2Y is training-set R2 and independent-test Q2Y uses the training-response mean in its denominator. For multivariate responses, each response is centered separately before the sums of squares are aggregated. In pls.single.cv(), each held-out prediction is evaluated relative to the corresponding fold-training response means. pls.double.cv() applies the same rule to the outer folds. RMSD is also calculated from held-out predictions and is reported on the response scale, so lower values are better. R2Y is a training-set explained-variance estimate from one additional model fitted on the full dataset; set fit = FALSE to skip this extra fit when only cross-validated performance is needed. For classification, Q2Y is calculated from held-out dummy-coded PLS-DA response scores using fold-training class proportions, accuracy reports decoded-label accuracy, and R2Y is calculated from the full-data PLS-DA fit on the dummy-coded response scores. Dummy-response Q2Y and R2Y are not classification accuracy. The exact convention used by each function is also reported in its metrics$definitions element. RMSD is not used for classification.

Use selection = "R2Y" to select from the full-data fitted-response path or selection = "Q2Y" to select from out-of-fold predictions standardized by each fold’s training-response mean. Selecting R2Y automatically enables and returns the fitted path even if fit = FALSE; because it is a training criterion, Q2Y or another held-out metric is usually preferable for choosing model complexity. The former names "r2" and "q2" are rejected because they did not identify these two different quantities unambiguously.

Classification settings can be selected by accuracy, balanced_accuracy, AUROC, lift_accuracy, macro_precision, macro_recall, macro_f1, kappa, R2Y, or Q2Y. Binary AUROC pools continuous held-out class scores across folds and treats the second factor level as positive. Use balanced accuracy when unequal class frequencies make majority-class accuracy misleading. Regression settings can be selected by R2Y, Q2Y, RMSD, MAE, MAPE_percent, RPD, Pearson_r, or Spearman_r. These additional aggregate metrics are useful for high-dimensional multivariate responses such as spectra. RMSD, MAE, and MAPE_percent are minimized; the other criteria are maximized. Signed bias and signed MRE_percent remain available from evaluate() but are not tuning criteria because they have no unambiguous one-sided optimization direction. The definitions match evaluate(), and an incompatible choice such as classification accuracy for a regression response raises an error before model fitting.

data.frame(
    ncomp = cv_opt$ncomp,
    training_R2Y = round(cv_opt$R2Y, 3),
    heldout_Q2Y = round(cv_opt$Q2Y, 3),
    heldout_RMSD = round(cv_opt$RMSD, 3)
)
#>         ncomp training_R2Y heldout_Q2Y heldout_RMSD
#> ncomp=1     1        0.742       0.687        3.471
#> ncomp=2     2        0.743       0.673        3.547
#> ncomp=3     3        0.774       0.461        4.555

Permutation-test p-values

fastPLS provides two permutation-test procedures. In pls(), the permutation test is a single train/test procedure: the rows of Xtrain are randomly permuted, the model is refitted, and the permuted test-set Q2Y values are compared with the observed Q2Y values component by component. The returned pval is the corrected Monte Carlo upper-tail value (b + 1) / (B + 1), where b counts successful null fits at least as extreme as observed and B counts successful null fits. It can therefore never be zero. Because pls() has no grouping argument, one training row is the permutation unit. The full permutation table is stored in permutation and can be visualized with plot.permutation(), where the x-axis is the correlation between the original and permuted response structure and the y-axis shows R2 and Q2.

perm_fit <- pls(
    Xtrain = Xreg_train,
    Ytrain = Ytrain_reg,
    Xtest = Xreg_test,
    Ytest = Ytest_reg,
    ncomp = 2,
    fit = TRUE,
    perm.test = TRUE,
    return_variance = FALSE,
    seed = 108
)

perm_fit$pval
#>    ncomp=2 
#> 0.00990099
plot.permutation(perm_fit, ncomp = 2)

In pls.double.cv(), the permutation test repeats the complete nested cross-validation workflow. Independent rows are exchanged individually. When constrain identifies repeated observations, complete constraint blocks are exchanged only between groups with the same number of rows. Thus group sizes, within-group response structure, and class frequencies are preserved exactly. The observed and permuted analyses use identical outer and inner folds and the same randomized-SVD seeds, so null variation reflects the exchange operation rather than new folds or sketches. Model selection and permutation inference use the same selection. For imbalanced classification, selection = "balanced_accuracy" therefore selects the model by mean class-specific recall and tests that statistic against its permutation distribution; it does not substitute dummy-response Q2. The output records the statistic in permutation_metric, permutation_observed, and permutation_sampled. Larger predictive metrics use the upper permutation tail, whereas losses such as RMSD use the lower tail. Both tails use (b + 1) / (B + 1). Failed null fits are stored in permutation_errors, omitted from B, and summarized by permutation_completed and permutation_failed.

dcv_perm <- pls.double.cv(
    Xdata = Xreg_train[1:20, ],
    Ydata = Ytrain_reg[1:20],
    ncomp = 1:2,
    kfold_inner = 2,
    kfold_outer = 2,
    perm.test = TRUE,
    seed = 109
)

data.frame(
    permutation_metric = dcv_perm$permutation_metric,
    observed = dcv_perm$permutation_observed,
    p_value = dcv_perm$p.value,
    completed = dcv_perm$permutation_completed,
    failed = dcv_perm$permutation_failed
)
#>   permutation_metric observed    p_value completed failed
#> 1               RMSD 3.293571 0.00990099       100      0

For an imbalanced classification analysis, the corresponding call is:

dcv_balanced <- pls.double.cv(
    Xdata = Xtrain,
    Ydata = Ytrain_cls,
    ncomp = 1:5,
    kfold_inner = 5,
    kfold_outer = 5,
    constrain = patient_id,
    selection = "balanced_accuracy",
    perm.test = TRUE,
    times = 100,
    seed = 109
)
dcv_balanced$balanced_accuracy
dcv_balanced$metrics$permutation

Cross-validation results do not retain the training matrices. To fit the final model, combine the selected configuration with best_parameters and pass the original training data explicitly to pls(). This keeps the CV object compact and makes the data used for refitting unambiguous.

cv_select <- pls.single.cv(
    Xdata = Xtrain,
    Ydata = Ytrain_cls,
    ncomp = 1:3,
    kfold = 5,
    seed = 106
)

selected <- utils::modifyList(
    cv_select$tuning_config,
    cv_select$best_parameters
)
svd_controls <- selected$svd_dots
selected$svd_dots <- NULL
fit_selected <- do.call(
    pls,
    c(
        list(
            Xtrain = Xtrain,
            Ytrain = Ytrain_cls,
            Xtest = Xtest,
            Ytest = Ytest_cls,
            return_variance = FALSE
        ),
        selected,
        svd_controls
    )
)

data.frame(
    best_ncomp = cv_select$best_ncomp,
    test_accuracy = mean(fit_selected$Ypred[[1]] == Ytest_cls)
)
#>   best_ncomp test_accuracy
#> 1          3     0.8333333

For example, kernelpls can select the best combination of component count and kernel setting. The selected values are returned in best_parameters and can be combined with tuning_config for the final explicit refit.

cv_kernel <- pls.single.cv(
    Xdata = Xtrain,
    Ydata = Ytrain_cls,
    ncomp = 1:3,
    kfold = 5,
    method = "kernelpls",
    kernel = c("linear", "rbf"),
    gamma = c(0.1, 1),
    seed = 107
)

selected_kernel <- utils::modifyList(
    cv_kernel$tuning_config,
    cv_kernel$best_parameters
)
svd_controls <- selected_kernel$svd_dots
selected_kernel$svd_dots <- NULL
fit_kernel <- do.call(
    pls,
    c(
        list(
            Xtrain = Xtrain,
            Ytrain = Ytrain_cls,
            Xtest = Xtest,
            Ytest = Ytest_cls,
            return_variance = FALSE
        ),
        selected_kernel,
        svd_controls
    )
)

data.frame(
    best_ncomp = cv_kernel$best_parameters$ncomp,
    best_kernel = cv_kernel$best_parameters$kernel,
    test_accuracy = mean(fit_kernel$Ypred[[1]] == Ytest_cls)
)
#>   best_ncomp best_kernel test_accuracy
#> 1          3         rbf             1

Nested or double CV

Double cross-validation is a nested validation design for separating model optimization from the final estimate of predictive performance. This separation is especially important for PLS-DA because the number of latent variables and other modelling choices can otherwise be tuned on the same samples used to report performance, producing optimistic accuracy or Q2 estimates. In the terminology of Szymanska et al. (2012), the inner CV loop (CV1) is used to optimize model complexity, such as the number of latent variables, whereas the outer CV loop (CV2) holds out samples that are not used during optimization and therefore provides the performance estimate of the complete modelling strategy. In fastPLS, pls.double.cv() follows this idea using a reproducible grouped- fold plan constructed in R. For the eligible selection criteria described above, CPU and Metal use a compiled coordinator for inner selection, outer refitting, prediction, and metric accumulation. CUDA uses an R coordinator around native CUDA single-CV and outer-fit kernels. Both kfold_inner and kfold_outer can be ordinary fold counts or "loocv", and both respect constrain.

cv_double <- pls.double.cv(
    Xdata = Xtrain,
    Ydata = Ytrain_cls,
    constrain = patient_id,
    ncomp = 1:3,
    kfold_inner = 3,
    kfold_outer = 3,
    method = "simpls",
    classifier = "lda",
    selection = "balanced_accuracy",
    backend = "cpu",
    return_splits = TRUE,
    seed = 104
)

data.frame(
    selected_ncomp_mode = cv_double$bcomp,
    outer_metric = cv_double$metric_name,
    outer_accuracy = cv_double$accuracy,
    outer_balanced_accuracy = cv_double$balanced_accuracy,
    outer_Q2Y = cv_double$Q2Y,
    outer_R2Y = cv_double$R2Y
)
#>   selected_ncomp_mode      outer_metric outer_accuracy
#> 1                   3 balanced_accuracy          0.975
#>   outer_balanced_accuracy outer_Q2Y outer_R2Y
#> 1                   0.975 0.5503804 0.6044152

data.frame(
    outer_fold = seq_along(cv_double$results[[1]]$best_ncomp),
    inner_selected_ncomp = cv_double$results[[1]]$best_ncomp
)
#>   outer_fold inner_selected_ncomp
#> 1          1                    3
#> 2          2                    3
#> 3          3                    3

head(cv_double$split_index)
#>   run_1_outer_1 run_1_outer_1_inner_1 run_1_outer_1_inner_2
#> 1 "test"        "outer_test"          "outer_test"         
#> 2 "test"        "outer_test"          "outer_test"         
#> 3 "training"    "training"            "training"           
#> 4 "training"    "training"            "training"           
#> 5 "training"    "test"                "training"           
#> 6 "training"    "test"                "training"           
#>   run_1_outer_1_inner_3 run_1_outer_2 run_1_outer_2_inner_1
#> 1 "outer_test"          "training"    "test"               
#> 2 "outer_test"          "training"    "test"               
#> 3 "test"                "training"    "training"           
#> 4 "test"                "training"    "training"           
#> 5 "training"            "test"        "outer_test"         
#> 6 "training"            "test"        "outer_test"         
#>   run_1_outer_2_inner_2 run_1_outer_2_inner_3 run_1_outer_3
#> 1 "training"            "training"            "training"   
#> 2 "training"            "training"            "training"   
#> 3 "training"            "test"                "test"       
#> 4 "training"            "test"                "test"       
#> 5 "outer_test"          "outer_test"          "training"   
#> 6 "outer_test"          "outer_test"          "training"   
#>   run_1_outer_3_inner_1 run_1_outer_3_inner_2 run_1_outer_3_inner_3
#> 1 "training"            "test"                "training"           
#> 2 "training"            "test"                "training"           
#> 3 "outer_test"          "outer_test"          "outer_test"         
#> 4 "outer_test"          "outer_test"          "outer_test"         
#> 5 "training"            "training"            "test"               
#> 6 "training"            "training"            "test"

Here, all rows from one illustrative patient remain together in both loops. If an inner training fold contains one class, fastPLS predicts that class and uses a finite constant score for every requested component. These predictions remain in the pooled inner metric. If an outer training partition contains one class, the same rule supplies its held-out prediction and the fold is marked as not estimable for discrimination. The returned degenerate_inner_folds, constant_classifier_fallback, minimum class counts, component_selection_informative, and outer_discrimination_estimable fields make these cases explicit.

The same nested workflow accepts a numeric vector for univariate regression or a numeric matrix for multivariate regression. The default inner-loop criterion for both tasks is held-out RMSD; for a response matrix it is aggregated across response columns.

cv_double_uni <- pls.double.cv(
    Xdata = Xreg_train,
    Ydata = Ytrain_reg,
    ncomp = 1:2,
    kfold_inner = 2,
    kfold_outer = 2,
    method = "simpls",
    backend = "cpu",
    seed = 105
)

cv_double_multi <- pls.double.cv(
    Xdata = Xmulti_train,
    Ydata = Ymulti_train,
    ncomp = 1:2,
    kfold_inner = 2,
    kfold_outer = 2,
    method = "plssvd",
    backend = "cpu",
    bycol = TRUE,
    seed = 105
)

data.frame(
    task = c("univariate regression", "multivariate regression"),
    selected_ncomp = c(cv_double_uni$bcomp, cv_double_multi$bcomp),
    outer_RMSD = c(cv_double_uni$RMSD, cv_double_multi$RMSD)
)
#>                      task selected_ncomp outer_RMSD
#> 1   univariate regression              1   4.177333
#> 2 multivariate regression              1   2.665607

For each outer fold, the inner loop selects the component count by balanced accuracy and the selected LDA model predicts samples that were not involved in that choice. Component counts may therefore differ among outer folds; bcomp reports their most frequent value, while balanced_accuracy and Ypred summarize the outer held-out predictions.

The same interface can tune more than the component count. In this example the inner loop chooses the best combination of ncomp, kernel, and gamma for kernelpls. The outer loop then uses the selected combination for each held-out outer fold.

cv_double_kernel <- pls.double.cv(
    Xdata = Xtrain,
    Ydata = Ytrain_cls,
    constrain = patient_id,
    ncomp = 1:2,
    kfold_inner = 3,
    kfold_outer = 3,
    method = "kernelpls",
    kernel = c("linear", "rbf"),
    gamma = c(0.1, 1),
    seed = 108
)

cv_double_kernel$results[[1]]$best_parameters
#> [[1]]
#> [[1]]$ncomp
#> [1] 2
#> 
#> [[1]]$kernel
#> [1] "rbf"
#> 
#> [[1]]$gamma
#> [1] 1
#> 
#> 
#> [[2]]
#> [[2]]$ncomp
#> [1] 2
#> 
#> [[2]]$kernel
#> [1] "rbf"
#> 
#> [[2]]$gamma
#> [1] 1
#> 
#> 
#> [[3]]
#> [[3]]$ncomp
#> [1] 2
#> 
#> [[3]]$kernel
#> [1] "rbf"
#> 
#> [[3]]$gamma
#> [1] 1

The most commonly used pls.double.cv() output fields are:

Important elements returned by pls.double.cv() are:

• results: one entry per repeated outer CV run. Each entry stores Ypred and pred for that run, the outer fold assignment, the best_ncomp selected inside each outer fold, the full best_parameters selected inside each outer fold, the complete inner-CV objects in inner, the run-level metric_name and metric_value, and the fitted backend and method.
• Ypred: the final cross-validated prediction for each sample. For classification, repeated runs are combined by voting; for regression, predictions are averaged across repeated runs.
• acc_tot: classification-only text summary of the total number and percentage of correctly classified samples.
• conf: classification-only confusion matrix. Entries are printed as counts and column percentages so that class-wise errors can be inspected.
• vote_counts: classification-only matrix with one row per sample and one column per class, showing how many repeated outer-CV runs voted for each class.
• accuracy, Q2Y, RMSD, and R2Y: one value per repeated outer CV run. For classification, accuracy reports decoded-label accuracy, Q2Y reports held-out Q2 on dummy-coded PLS-DA response scores, R2Y reports the mean training-fit R2 of the selected outer-fold PLS-DA models, and RMSD is not used. For regression, Q2Y reports held-out Q2, RMSD reports held-out RMSD, and R2Y reports the mean training-fit R2 of the selected outer-fold models.
• metric_name: the metric used for inner/outer model selection. It is based on held-out predictions except when selection = "R2Y" explicitly requests the fitted-response criterion.
• medianQ2Y, CI95Q2Y, medianR2Y, CI95R2Y, medianRMSD, and CI95RMSD: summaries across repeated outer CV runs, returned only when runn > 1. They are omitted for the default runn = 1 output.
• bcomp: the most frequently selected number of components across all outer folds and repeated runs.
• backend and method: the default backend and PLS method used by the call. If vector-valued methods or backends are tuned in the inner loop, the selected fold-level values are stored in results[[run]]$best_parameters.
• selection_metric: the criterion used by the inner CV loop. The default "auto" means accuracy for classification and RMSD for regression. All task-compatible selection criteria listed above can also drive the nested permutation test.