Partial least squares (PLS) regression, introduced by [1], is a well-known dimension-reduction method, notably in chemometrics and spectrometric modeling [2]. In this paper, we focus on the PLS univariate response framework, better known as PLS1. Let n be the number of observations and p the number of covariates. Then, y = ( y 1 , … , y n ) T ∈ ℝ n represents the response vector, with ( . ) T denoting the transpose. The original underlying algorithm, developed to deal with continuous responses, consists of building latent variables t k ,   1 ⩽ k ⩽ K , also called components, as linear combinations of the original predictors X = ( x 1 , … , x p ) ∈ ℳ n , p ( ℝ ) , where ℳ n , p ( ℝ ) represents the set of matrices of n rows and p columns. Thus, t k = X k − 1 w k ,   1 ⩽ k ⩽ K , (1) where X 0 = X , and X k − 1 ,   k ⩾ 2 represents the residual covariate matrix obtained through the ordinary least squares regression (OLSR) of X k − 2 on t k − 1 . Here, w k = ( w 1 k , … , w p k ) T ∈ ℝ p is obtained as the solution of the following maximization problem [3]: w k = argmax w ∈ ℝ p { Cov 2 ( y k − 1 , t k ) } (2) = argmax w ∈ ℝ p { w T X k − 1 T y k − 1 y k − 1 T X k − 1 w } , (3) with the constraint ‖ w k ‖ 2 2 = 1 , and where y 0 = y , and y k − 1 ,   k ⩾ 2 represents the residual vector obtained from the OLSR of y k − 2 on t k − 1 .

The final regression model is thus: y = ∑ k = 1 K c k t k + ϵ (4) = ∑ j = 1 p β j PLS x j + ϵ , (5) with ϵ = ( ϵ 1 , … , ϵ n ) T the n × 1 error vector and ( c 1 , … , c K ) the regression coefficients related to the OLSR of y k − 1 on t k , ∀ k ∈ 〚 1 , K 〛 , also known as y-loadings. More details are available, notably in Höskuldsson [4] and Tenenhaus [5].

This particular regression technique, based on reductions in the original dimension, avoids matrix inversion and diagonalization, using only deflation. It allows us to deal with high-dimensional datasets efficiently, and notably solves the collinearity problem [6].

With great technological advances in recent decades, PLS regression has been gaining attention, and has been successfully applied to many domains, notably genomics. Indeed, the development of both microarray and allelotyping techniques result in high-dimensional datasets from which information has to be efficiently extracted. To this end, PLS regression has become a benchmark as an efficient statistical method for prediction, regression and dimension reduction [3]. Practically speaking, the observed response related to such studies does commonly not follow a continuous distribution. Frequent goals with gene expression datasets involve classification problems, such as cancer stage prediction, disease relapse, and tumor classification. For such reasons, PLS regression has had to be adapted to take into account discrete responses. This has been an intensive research subject in recent years, leading globally to two types of adapted PLS regression for classification. The first, studied and developed notably by Nguyen and Rocke [7], Nguyen and Rocke [8] and Boulesteix [9], is a two-step method. The first step consists of building standard PLS components by treating the response as continuous. In the second step, classification methods are run, e.g., logistic discrimination (LD) or quadratic discriminant analysis (QDA). The second type of adapted PLS regression consists of building PLS components using either an adapted version of or the original iteratively reweighted least squares (IRLS) algorithm, followed by generalized linear regression on these components. This type of method was first introduced by Marx [10]. Different modifications and improvements, using notably ridge regression [11] or Firth’s procedure [12] to avoid non-convergence and infinite parameter-value estimations, have been developed, notably by Nguyen and Rocke [13], Ding and Gentleman [14], Fort and Lambert-Lacroix [15] and Bastien et al. [16]. In this work, we focus on the second type of adapted PLS regression, referred to from now on as GPLS.

As previously mentioned, a feature of datasets of interest is their high-dimensional setting, i.e., n ≪ p . Chun and Keleş [17] have shown that the asymptotic consistency of PLS estimators does not hold in this situation, so filtering or predictor selection become necessary in order to obtain consistent parameter estimation. However, all methods described above proceed to classification using the entire set of predictors. For datasets that frequently contain thousands of predictors, such as microarray or RNAseq ones, a variable filtering pre-processing thus needs to be applied. A commonly used pre-processing method when performing classification uses the BSS/WSS-statistic: BSS j / WSS j = ∑ q = 1 Q ∑ i : y i ∈ G q ( μ ^ j q − μ ^ j ) 2 ∑ q = 1 Q ∑ i : y i ∈ G q ( x i j − μ ^ j q ) 2 , (6) with μ ^ j the sample mean of x j , and μ ^ j q the sample mean of x j in class G q for q ∈ { 1 , … , Q } . Then, predictors associated with the highest values of this are retained, but no specific rule exists to choose the number of predictors to retain. A Bayesian-based technique, available in the R-package limma, has become a common way to deal with this, computing a Bayesian-based p-value for each predictor, therefore allowing users to choose the number of relevant predictors based on a threshold p-value [18].

This method cannot be considered parsimonious, but rather as a pre-processing stage for exclusion of uninformative covariates. Reliably selecting relevant predictors in PLS regression is of interest for several reasons. Practically speaking, it would allow users to identify the original covariates which are significantly linked to the response, as is done in OLS regression with Student-type tests. Statistically speaking, it would avoid the establishment of over-complex models and ensure consistency of PLS estimators. Several methods for variable selection have already been developed [19]. Lazraq et al. [20] group these techniques into two main categories. The first, model-wise selection, consists of first establishing the PLS model before performing a variable selection. The second, dimension-wise selection, consists of selecting variables on one PLS component at a time.

A dimension-wise method, introduced by Chun and Keleş [17] and called Sparse PLS (SPLS), has become the benchmark for selecting relevant predictors using PLS methodology. The technique is for continuous responses and consists of simultaneous dimension reduction and variable selection, computing sparse linear combinations of covariates as PLS components. This is achieved by introducing an L 1 constraint on the weight vectors w k , leading to the following formulation of the objective function: w k = argmax w ∈ ℝ p { w T X k − 1 T y k − 1 y k − 1 T X k − 1 w } , s . c .     ‖ w ‖ 2 2 = 1 ,   ‖ w ‖ 1 < λ , (7) where λ determines the amount of sparsity. More details are available in Chun and Keleş [17]. Two tuning parameters are thus involved: η ∈ [ 0,1 ] as a rescaled parameter equivalent to λ, and the number of PLS components K, which are determined through CV-based mean squared error (CV MSE). We refer to this technique as SPLS CV in the following.

Chung and Keles [21] have developed an extension of this technique by integrating it into the generalized linear model (GLM) framework, leading it to be able to solve classification problems. They also integrate Firth’s procedure, in order to deal with non-convergence issues. Both tuning parameters are selected using CV MSE. We refer to this method as SGPLS CV in the following.

A well-known model-wise selection method is the one introduced by Lazraq et al. [20]. It consists of bootstrapping pairs ( y i , x i • ) ,   1 ⩽ i ⩽ n , where x i • represents the ith row of X , before applying PLS regression with a preset number of components K on each bootstrap sample. By performing this method, approximations of distributions related to predictors’ coefficients can be achieved. This leads to bootstrap-based confidence intervals (CI) for each predictor, and so opens up the possibility of directly testing their significance with a fixed type I risk α. The advantages of this method are twofold. First, by focusing on PLS regression coefficients, it is relevant for both the PLS and GPLS frameworks. Second, it only depends on one tuning parameter K, which must be determined earlier.

Two important related issue should be mentioned. Firstly, the presence of variables unrelated to the phenomenon under study strongly increases the level of inherent noise in the response and in the predictors. This phenomenon occurs frequently in various contexts and has been studied by several authors both in the context of PLS regression and its extensions [22, 23] and in the more general context of machine learning where the use of techniques related to least squares have proven to be relevant and effective [24–28].

Secondly, while performing this technique, approximations of distributions are achieved conditionally on the fixed dimension of the extracted subspace. In other words, it approximates the uncertainty of these coefficients, conditionally on the fact that estimations are performed in a K-dimensional subspace for each bootstrap sample. The determination of an optimal number of components is crucial for achieving reliable estimations of the regression coefficients [29]. Thus, since this determination is specific to the dataset in question, it must be performed for each bootstrap sample, in order to obtain reliable bootstrap-based CI. We have established some theoretical results which confirm this (Supplementary Section 1.3).

Determining tuning parameters by using q-fold cross-validation (q-CV) based criteria may induce important issues concerning the stability of extracted results (Hastie et al. [30], p.249; Boulesteix [31]; Magnanensi et al. [32]). Thus, using such criteria for successive choosing of the number of components should be avoided. As mentioned, amongst others, by Efron and Tibshirani [33], p.255 and Kohavi [34], bootstrap-based criteria are known to be more stable than CV-based ones. In this context, Magnanensi et al. [32] developed a robust bootstrap-based criterion for the determination of the number of PLS components, characterized by a high level of stability and suitable for both the PLS and GPLS regression frameworks. Thus, this criterion opens up the possibility of reliable successive choosing for each bootstrap sample.

In this article, we introduce a new dynamic bootstrap-based technique for covariate selection suitable for both the PLS and GPLS frameworks. It consists of bootstrapping pairs ( y i , x i • ) , and successive extraction of the optimal number of components for each bootstrap sample, using the previously mentioned bootstrap-based criterion. Here, our goal is to better approximate the uncertainty related to regression coefficients by removing the condition of extracting a fixed K-dimensional subspace for each bootstrap sample, leading to more reliable CI. This new method both avoids the use of CV, and features the same advantages as those previously mentioned related to the technique introduced by Lazraq et al. [20]. This method also has the advantage of behaving consistently with the variation of the noise level in the dataset. We refer to this new dynamic method as BootYTdyn in the following.

We also succeed in adapting the bootstrap-based criterion introduced by Magnanensi et al. [32] to the determination of a unique optimal number of components related to each preset value of η in both the SPLS and SGPLS frameworks. Thus, these adapted versions, by reducing the use of CV, improve the reliability of the hyper-parameter tuning. We will refer to these adapted techniques as SPLS BootYT and SGPLS BootYT, respectively.

There is some constructive attributes of our proposed approach as follows:

1. theoretical results are developed and presented to give indications on the number of components related to a bootstrap sample and showed the necessity for a dynamic algorithm for Sparse PLS;

The article is organized as follows. In Section 2, we introduce our new dynamic bootstrap-based technique, followed by the description of our adaptation of the BootYT stopping criterion to the SPLS and SGPLS frameworks. In Section 3, we present simulations related to the PLS framework, and summarize the results, depending notably on different noise levels in y. In Section 4, we treat a real microarray gene expression dataset with a binary response, here the original localization of colon tumors, by benchmarking our new dynamic bootstrap-based approach for GPLS regression. Lastly, we discuss results and conclusions in Section 5 as well as possible future work.

In this section, we deal with the predictors matrix introduced in Section 3.3.1 and the original binary response vector. Five approaches for variable selection, adapted for the GPLS framework, are considered for comparison.

1. BootYT. The bootstrap-based method, introduced by Lazraq et al. [20], combined with the bootstrap-based criterion [32] for pre-selecting the number of components.

Concerning bootstrap-based approaches, the incorporation of PLS methodology into GLM, developed by Bastien et al. [16], was used. Due to non-convergence issues for parameter estimations, some samples were excluded using a threshold for parameter estimations. Indeed, a model built using a bootstrap sample with at least one parameter estimate that is higher in absolute value than 10 4 times the one (in absolute value as well) estimated on the original dataset, leads to the exclusion of this bootstrap sample. Thus, to ensure a sufficient number of relevant bootstrap samples, the preset number of computed samples was increased to R = 4000 . As for the lasso, three different loss functions for the establishment of the sparsity parameter using 10-fold CV were used: the number of misclassified values, the MSE, and the deviance. We refer to them as Lasso.Cl, Lasso.MSE and Lasso.Dev, respectively. The set of values defined for the sparsity parameter is preset by the “glmnet” package. For the SGPLS and RSGPLS approaches, the number of components K varies from 1 to 10, the sparsity parameter η varies from 0.04 to 0.99 (by steps of 0.05), and the ridge parameter involved in the RSGPLS technique is selected from 31 points log 10 -linearly spaced in the range [ 10 − 2 , 10 3 ] , as in Durif et al. [38].

Each method was performed one hundred times in order to obtain relevant results. However, due to the high observed stability of results extracted with the BootYTdyn approach, and in order to save computational time, this was only performed twenty times instead of a hundred.

Stability results for tuning parameter selection, except for the bootstrap-based methods, are shown in Figure 4 and Table 7.

Based on results summarized in Table 7, the lasso methodology, using CV-based MSE or deviance values for the selection of its hyper-parameter, is the most stable. This can be explained by the fact that only one hyper-parameter is involved, while the other techniques are based on two or three tuning parameters. Using the number of misclassified values for CV has to be avoided for stable selection of the sparsity parameter. Our SGPLS adaptation with the BootYT criterion improves reliability in selecting the set of tuning parameters, as previously observed in the PLS framework (Sections 2.2, 3.3.4). Note that, concerning RSGPLS, three different sets of optimal parameters were extracted with maximal occurrence rate of five, all of them selecting the same set of predictors. Thus, the set of parameters which leads to the smallest number of misclassified values on the training dataset was retained. As already mentioned in Section 3.3.3, extracting the same support does not necessarily lead to the same model when sparsity or ridge parameters are involved. Therefore, the numbers of sparse supports and models retained are summarized in Table 8.

In light of these results, BootYTdyn and both the MSE- and deviance-based lasso techniques are the most stable in extracting supports and models. This could be due in part to the fact than all of them depend on only one tuning parameter. Even if this hyper-parameter for the BootYTdyn approach, i.e., the number of components, has to be chosen R times, the high stability of the bootstrap-based stopping criterion introduced by Magnanensi et al. [32] endows this approach with good stability in selecting the sparse support. As for the PLS framework, our new bootstrap-based SGPLS implementation improves the stability here. The lack of stability of the lasso based on misclassified values is directly induced by the discrete form of this statistic. This issue was also observed and mentioned in Magnanensi et al. [32].

The various selected supports are displayed in Figure 5. Note that both the MSE- and deviance-based lasso regressions select the same support and the same model, noted Lasso.MSE.Dev in the following.

In the following, two additional independent public datasets, stated by Marisa et al. [35] as being comparable with our original dataset, were included for comparison of predictive abilities. These datasets are named GSE18088 (n = 53) [39] and GSE14333 (n = 247) [40]. Both the MSE and number of misclassified (MC) value of the selected models were computed on both the training and test parts of the original dataset, as well as on the two additional datasets. These results are summarized in Table 9. Since independent datasets were available, we decided to follow the suggestion of Van Wieringen et al. [41, p.1596]: “the true evaluation of a predictor’s performance is to be done on independent data”.

In this real data study, BootYTdyn retains only one predictor, which is also retained by all other methods. Thus, as expected, this most sparse support induced the highest values of both the MSE and number of misclassified observations, based on the training subset of the original dataset. It also provided the highest values based on the test subset of the original dataset, which is to be expected since, as explained in Section 3.3.1, these two parts are well-balanced in terms of anatomo-clinical characteristics. Thus, this causes a bias in evaluation of model’s predictive abilities, by making comparable MSE and misclassified values based on both the training and testing subsets. Results in Table 9 confirm this property, and also highlight the usefulness of additional independent datasets for reliably comparing predictive abilities. While a well-designed model for predictive purposes will provide similar PMSE values on comparable independent additional datasets compared to MSE obtained on the training dataset, an over-fitted one would be related to higher PMSE values, due to its dependance on training data. Thus, the differences, noted Δ MSE, between the MSE obtained on the training subset of the original dataset, and those obtained on the three test datasets, are shown in Figure 6.

Our new dynamic bootstrap-based approach is the only one that exhibits this MSE stability property. In contrast, all other approaches provided higher PMSE values on the two other datasets than on the original. This led BootYTdyn to have only 48 (16%) misclassified values on these two additional datasets, representing the best result among all studied approaches. Thus, we can reasonably assume that our new method has helped to remove non-informative predictors, in the sense of not being relevant for improving predictive ability.

Lastly, the unique extracted probe set, named 230784_at, is already known to be related to the original location in the distal colon. The sign of the regression coefficient obtained with BootYTdyn is coherent with this state-of-the-art result, and thus strengthens our conclusion.

In this article, we developed a new dynamic bootstrap-based technique for variable selection, suitable for both the PLS and GPLS frameworks, and proposed a bootstrap-based version of the SPLS and SGPLS methods for selecting the number of components. While the first of these lets us completely avoid CV, the second lets us select the set of tuning parameters in a more reliable way.

In state-of-the-art approaches, the use of CV-based techniques for selection of hyper-parameters is common, and can lead to important stability issues, observed notably by Boulesteix [31] and Magnanensi et al. [32], and confirmed in our studies. Our new dynamic bootstrap-based technique represents a useful method for avoiding CV-related issues, since the unique hyper-parameter is successively selected using a bootstrap-based criterion. This new technique improves on the original bootstrap-based methodology introduced by Lazraq et al. [20], in that it permits us to approximate the distribution of covariates’ regression coefficients by removing the condition of working in a subspace of fixed dimension K. Theoretical results have been established that strengthen the usefulness of building subspaces spanned by a dynamic number of components for performing PLS regressions on bootstrap samples.

In the PLS framework, conclusions based on our simulations are twofold. First, for datasets with negligible random noise in y , SPLS BootYT is recommended. Indeed, both in terms of accuracy and predictive abilities, it outperforms all other techniques compared here. Second, for datasets with non-negligible random noise in y, which represents the more realistic case, our new dynamic bootstrap-based method is recommended. As for SPLS BootYT for the negligible random noise variability framework, BootYTdyn outperform all other methods for each property studied. Furthermore, it is the only method which concludes as to a decreasing number of significant predictors, while the random noise variability increases, which was expected here.

Results obtained from our classification study using real datasets match with previous conclusions. Indeed, BootYTdyn is the only one which leads to the expected PMSE values on two additional independent datasets, which allows us to suppose that over-fitting issues were avoided, and that these PMSE are induced by noise or information which cannot be modeled using only gene expression. Furthermore, the extracted probe set is already known to be linked to the relevant location in the distal colon, which strengthens our confidence in this new dynamic approach.

Lastly, our new bootstrap-based SPLS implementation improves the stability of this method. Indeed, in all cases studied, both the SPLS CV and SGPLS CV select a higher number of unique sets of hyper-parameters than our bootstrap-based versions do, leading to higher numbers of unique supports and models too.

In the future, our work may be extended in the following three directions. Firstly, Sun et al. [42] also worked on this subject, and proposed a methodology for selecting tuning parameters of penalized regressions in order to stabilize variable selection. Even if their method is not applicable to the selection of the number of components in PLS regression, it would be interesting to adapt it to both SPLS CV and our new SPLS BootYT implementation, in order to look for potential stability gains.

Secondly combining partial least squares based regular or sparse techniques with the existing techniques that were developed in [24–28] may speed up, specialize toward a specific target or improve the accuracy of those existing algorithms by taking into account some latent structure of the data.

Thirdly, simulations may be performed to strengthen the results obtained on real data sets for the logistic framework (Section 4). In addition, testing the performance of these new approaches for responses following other distributions -such as Poisson, Gamma-should also be done. However, based on all the results obtained in this paper, our new dynamic method seems to be the most efficient compared to the state of the art approaches for datasets with non-negligible noise variability, a common situation in daily practice.