Hadrien Lorenzo
email: hadrien.lroenzo@inria.fr

Olivier Cloarec

SISTM Rodolphe Thiébaut

Data-Driven Sparse Partial Least Squares

Keywords: PLS regression, Supervised learning, Variable selection, Soft thresholding, Multi-block data

In the supervised high dimensional settings with a large number of variables and a low number of individuals, variable selection allows a simpler interpretation and more reliable predictions. That subspace selection is often managed with supervised tools when the real question is motivated by variable prediction. We propose a Partial Least Square (PLS) based method, called data-driven sparse PLS (ddsPLS), allowing variable selection both in the covariate and the response parts using a single hyper-parameter per component. The subspace estimation is also performed by tuning a number of underlying parameters. The ddsPLS method is compared to existing methods such as classical PLS and two well established sparse PLS methods through numerical simulations. The observed results are promising both in terms of variable selection and prediction performance. This methodology is based on new prediction quality descriptors associated with the classical R 2 and Q 2 and uses bootstrap sampling to tune parameters and select an optimal regression model.

Introduction

While recent biotechnologies offer more and more descriptors (denoted as p) for monitoring biological processes, the number of observations (denoted as n) does not increase, which implies the ill-posed problem n < < p in supervised settings (such as regression framework). Different methods have been developed in the past century to tackle this high dimensional problem and, Partial Least Squares (PLS), particularly gathers attention since its introduction by [START_REF] Wold | Estimation of principal components and related models by iterative least squares[END_REF]. Indeed, this methodology has firstly become a standard in chemometrics calibration [START_REF] Martens | Multivariate calibration[END_REF] and is now used in biology Zoong [START_REF] Zoong Lwe | Heat stress elicits remodeling in the anther lipidome of peanut[END_REF], animal genetics [START_REF] Kautt | Contrasting signatures of genomic divergence during sympatric speciation[END_REF], in human genetics [START_REF] Rechtien | Systems vaccinology identifies an early innate immune signature as a correlate of antibody responses to the ebola vaccine rvsvzebov[END_REF] and in many other areas. However, to keep the ease of regression models' interpretation as p grows further, variable selection has also started to be needed by PLS users. Following the L1-penalization of the regression coefficients introduced by [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF], models that handle variable selection solutions (denoted as sparse methodologies in the following) and produce model regularization allowing to keep confidence in the built models have been introduced by Lê Cao et al. (2008); [START_REF] Chun | Sparse partial least squares regression for simultaneous dimension reduction and variable selection[END_REF]. Even if [START_REF] Chun | Sparse partial least squares regression for simultaneous dimension reduction and variable selection[END_REF] proved that PLS consistency is lost as n shrinks and p grows, the PLS and sparse PLS methods have allowed authors to analyze real datasets meaningfully. The PLS estimation goes through R ∈ N identical processes which can be generally described as follows.

(a) Estimate the covariance matrix between covariate (x) and response (y) parts.

(b) Estimate the singular-space associated with the largest singular-value of the previous matrix.

(c) Project the covariate and response parts on the previously defined subspace.

(d) Estimate the linear regression matrix of the response part on the covariate part in the subspace.

(e) Remove the information carried by the current subspace from the covariate and response parts. The step (a) corresponds to the estimation of a covariance matrix between the response and the covariate parts. The chosen estimator is very often M = 1/n n i=1 (yi -ȳ)(xix) , for a data-set (xi, yi)i=1..n of size n with empirical means x and ȳ. This estimator is very sensitive to the curse of dimensionality, when the sample size is small and/or the number of variables is large, see for example [START_REF] Bickel | Regularized estimation of large covariance matrices[END_REF]. However, the works previously mentioned only apply regularization methods in step (b), trying to overcome the failure to estimate the covariance matrix.

Therefore, the present work focuses on a regularized estimation of the covariance matrix directly in step (a). More precisely, the considered solution and its associated optimization problem are

S λ (M) = arg min Σ∈R q×p ||M -Σ|| 2 + 2λ |Σ| , (1)
where the soft-thresholding operator, an element-wise matrix operator, is defined as ∀X ∈ R p,q , S λ (X) = (S λ (xi,j)) (i,j)∈[ [1,p]]×[ [1,q]] = (sign(xi,j) max(0, |xi,j| -λ)) (i,j)∈[ [1,p]]×[ [1,q]] ,

|| • || and | • | are respectively the Frobenius and the L1 norms. Note that if λ = 0, the chosen estimator is the empirical estimator of the variance-covariance matrix. So far, no theoretical development has been conducted to study the behavior of this estimator but for the case of a variance matrix, that is when x = y. Indeed, this case has been widely analyzed for the soft-thresholding operator but also for different thresholding operators. [START_REF] Johnstone | On consistency and sparsity for principal components analysis in high dimensions[END_REF] have explored the pertinence of thresholding the diagonal of the empirical covariance matrices in the objective of building sparse versions of the Principal Components Analysis (PCA). [START_REF] Bickel | Regularized estimation of large covariance matrices[END_REF] have put forward theoretical and numerical arguments to various types of thresholding operators. Also, [START_REF] Deshpande | Sparse PCA via covariance thresholding[END_REF] have shown the interest of such an idea where the whole empirical covariance matrix is thresholded, and the diagonal is removed.

In the PLS context, the number of components must be tuned, which is generally performed using prediction error. For sparse models, other parameters must also be adjusted by cross-validation and, more specifically, the leave-one-out method in the context of PLS. Recently, a sparse PLS approach (denoted as sgPLS for sparse group PLS and introduced in [START_REF] Sutton | Sparse partial least squares with group and subgroup structure[END_REF]; [START_REF] Liquet | Group and sparse group partial least square approaches applied in genomics context[END_REF]) has used bootstrap approach in [START_REF] Broc | Penalized partial least squares for pleiotropy[END_REF] to efficiently overcome the difficulty of the "large p, small n" multi-block framework to select optimal models. Recently, a meta-analysis of sparse PLS methodologies was performed by [START_REF] Mehmood | Comparison of variable selection methods in partial least squares regression[END_REF] over 16 different types of strategies divided in 3 classes which are "filter", "wrapper" and "embedded" methods. The first class "filter" corresponds to a filter which is applied on the output of a PLS model (on the regression coefficients for example) and allows to classify variables according to their importance. The second class "wrapper" considers methodologies where the result of a "filter" is used to refit a PLS model. This operation is produced a certain amount of times. The most classical methods in this context are the backward (and also forward even if this has not been studied in the PLS case yet) variable selection methodologies. It also takes into account randomized methodologies regarding to the studied variables or sub samples. This list is non exhaustive and we encourage the reader to refer to the section "3.2 Wrapper methods" of [START_REF] Mehmood | Comparison of variable selection methods in partial least squares regression[END_REF] to appreciate the richness of the solutions devised. The third class "embedded" corresponds to context where variable selection is performed in the PLS model building.

In the following, the proposed sparse PLS model, based on empirical covariance thresholding, softthresholding more specifically, is studied. The thresholds and the number of components are tuned thanks to bootstrap operations. Following the classification of sparse PLS methods described in [START_REF] Mehmood | Comparison of variable selection methods in partial least squares regression[END_REF], the proposed methodology is close to the third class "embedded". Yet this methodology modifies the PLS algorithm working on the covariance matrix before the PLS model is actually built, this closeness to the data gave its name: data-driven sparse PLS (ddsPLS). This ddsPLS methodology opens a new class of sparse PLS methods.

In Section 2, the underlying statistical model associated with latent variables is introduced. Section 3 describes the PLS (NIPALS-PLS/PLS2) algorithm and details the existing sparse PLS solutions. Section 4 introduces the proposed ddsPLS solution and the associated model selection algorithm. Based on numerical studies, Section 5 reports the quality of the proposed methodology and compares it to the above-presented methods through simulation designs. Finally, Section 6 provides some concluding remarks.

A latent variable model

Let x (resp. y) be a p-dimensional (resp. q-dimensional) random variable. Let assume that the variables of x and of y are centered and of unit variance. Let n ∈ N be the sample size. Let X and Y be n × p and n × q random matrices where, for each row i ∈ [ [1, n]], xi ∼ x and yi ∼ y. Let us further assume that the xi's, respectively the yi's, are independent to each other. Different methods such as PCA and PLS are based on the principle of generative latent vector common to x and y, denoted by φ hereafter, a R-dimensional random variable. The corresponding latent variable model can be written as

       x = A φ + where A = [a1, • • • , aR] ∈ MR×p(R) with ||ar|| = 0, y = D φ + ξ where D = [d1, • • • , dR] ∈ MR×q(R) with ||dr|| = 0, with var(φ) = IR, Eφ = 0R, E = 0p, Eξ = 0q, cov(φ, ) = 0R×q, cov(φ, ξ) = 0R×q, cov( , ξ) = 0p×q.
(2) This model implies that var(x) = A A + var( ) and var(y) = D D + var(ξ). Furthermore, if the components of the random errors and ξ are independent and share the same variance, denoted by σ 2 , it comes var(x) = A A + σ 2 I and var(y) = D D + σ 2 I. In all cases, the relation holds cov(y, x) = D A.

Remark 1 Helland and Almøy (1994) have studied the number R of latent variables shared by y and x. They specified in which measure R is different from R in the general case. They have shown that the dimension R is equal to the largest number of eigen-vectors of A A = var (x), which have non-null projections on A D = cov (x, y). In the PLS framework, R is the dimension of interest and is retrieved in the methodology as the estimated number of components.

The underlying latent variable model (2) intuitively implies the following regression model

y = B x + e, (3)
where B ∈ R p×q is a non-stochastic matrix, e ∈ R q is a centered residual q-dimensional random vector, and e and x are independent. Matrix B must satisfy

(AB) φ = D φ =⇒ E • φ AB = D. (4)
In the general case, without any additional hypothesis, B is not unique, but for the trivial cases or when AA is invertible, the ordinary least-squares estimator of B only uses Equation (3), minimizing the l2-norm of e, and provides B = X + Y, which is known to be unstable in the case of multicollinearity in the X matrix.

Let us notice that the latent variable model ( 2) is also not identifiable. Indeed, for any matrix 2) implies the following latent decomposition used by most of the authors.

G ∈ GL R (R), denoting by t = G φ, P = G -1 A and C = G -1 D, model (
   x = P t + , y = C t + ξ, with var(t , , ξ ) diagonal, ( 5
)
where t is called the score and P and C are called the loadings. Since the score and loadings are not identifiable, the PLS method concentrates on building a subspace that covers the space drawn by φ through an unknown matrix G. Note that, only when numerical simulations are carried out, Equation (4) allows to screen the quality of structural reconstruction of the current estimator B of B thanks to the metrics (for example)

||A B -D|| 2 ||D|| 2 , (6)
which is as close to 0 as the estimated model is close to the theoretical one.

The PLS algorithm is a multi-linear regression method able to estimate regression matrices even in the context of degenerate datasets, i.e. when X X is singular. The principle of this algorithm is to iteratively estimate, here for the r th -iteration, the first right and left singular vectors, ur and vr, of the empirical covariance matrix M (r) = Y (r) X (r) /(n -1) where Y (r) and X (r) are the residual matrices of the previous steps, see technical details in next section.

NIPALS-PLS and sparse PLS solutions

This section recalls three different PLS algorithms. First the original PLS-Nipals algorithm, denoted by NIPALS-PLS [START_REF] Wold | Estimation of principal components and related models by iterative least squares[END_REF] hereafter, is presented. Then the sparse PLS algorithm using L1-weights regularization, denoted by sPLS Lê [START_REF] Lê Cao | A sparse PLS for variable selection when integrating omics data[END_REF] hereafter, is detailed. Finally, a modification of the previous sparse solution to control the non-convex part of the optimization problem, denoted by SPLS [START_REF] Chun | Sparse partial least squares regression for simultaneous dimension reduction and variable selection[END_REF] hereafter, is provided.

NIPALS-PLS (PLS2) solution

The NIPALS-PLS [START_REF] Wold | Estimation of principal components and related models by iterative least squares[END_REF] algorithm, often referred to as PLS2 algorithm in the literature, solves the optimization problem given at the first row of Table 1 . The r-th iteration of this algorithm is described in (7): using the initializations X (1) = X and Y

(1) = Y, ∀r ∈ [[1, R]]              (a) ur = --→ RSV M (r) , vr = --→ RSV M (r) , (b) tr = X (r) ur, (c) pr = X (r) tr/t r tr, (d) cr = Y (r) tr/t r tr, (e) X (r+1) = X (r) -trp r , Y (r+1) = Y (r) -trc r , (7)
where --→ RSV (for Right-Singular-Vector) is the function that returns the first right singular vector of its argument. In this algorithm, step (a) estimates the weights, step (b) estimates the scores, steps (c) and (d) estimate the linear regression matrices of X (r) on tr and Y (r) on tr, and finally step (e) applies deflation. The regression matrix estimation is performed using the relation

B = U P U -1 C (8)
where U = [u1, . . . , uR], P = [p1, . . . , pR] and C = [c1, . . . , cR]. If one wants to build the scores from an x-matrix, potentially different from the one used to build the model, it must not use the matrix U directly but the matrix U (P U) -1 : XU (P U) -1 = (t1, . . . , tR) = (X (1) u1, . . . , X (R) uR).

Remark 2 In its original form, the PLS2/NIPALS-algorithm details the step (a) alternating estimations of ur and vr until convergence of ur such as ur = X (r) sr/||X (r) sr|| tr = X (r) ur/(u r ur) cr = Y (r) tr/(t r tr) sr = Y (r) cr/(c r cr)

,
where sr is initialised as the first column of Y (r) , see [START_REF] Tenenhaus | La régression PLS: théorie et pratique[END_REF] for example. This alternating algorithm and the --→ RSV operators both build ur and the form associated to the operator will be kept in the following due to its compact aspect.

Two sparse PLS solutions

Variable selection can be performed during the computation of the components used to construct the PLS models. Based on L1-penalization, two sparse PLS solutions have been developed, denoted sPLS Lê [START_REF] Lê Cao | A sparse PLS for variable selection when integrating omics data[END_REF] and SPLS Chun and Keleş (2010) hereafter. More precisely, the NIPALS-PLS maximization criterion is then modified such as detailed in the second and third rows, respectively, of Table 1 . sPLS Lê [START_REF] Lê Cao | A sparse PLS for variable selection when integrating omics data[END_REF] introduces sparsity both in the X and in the y parts through 2R+1 parameters. The underlying idea is to construct a solution using the soft threshold operator S λ (•). The r-th step of the associated algorithm is detailed in (9): using the initializations X (1) = X and Y (1) = Y,

∀r ∈ [[1, R]]                         
(a † ) Initialize vr normalized, iterate until convergence of the two steps :

(i † ) ur = S λ (r) u M (r) vr , ur = ur/||ur||, (ii † ) vr = S λ (r) v M (r) ur , vr = vr/||vr||, (b) tr = X (r) ur, (c) pr = X (r) tr/t r tr, (d) cr = Y (r) tr/t r tr, (e) X (r+1) = X (r) -trp r , Y (r+1) = Y (r) -trc r , (9)
where the regression matrix estimation is performed using the relation (8). One might notice that λ (r) v fixes the sparsity of vr. However, since Equation (e) in Algorithm (7) does not explicitly impose sparsity in the y part, then the estimated matrix B is not constrained to be sparse in Y. In the implementation, the search for parameters (λ

(r) u , λ (r)
v ) is driven, in an equivalent way, by two new parameters (keepX , keepY )r which respectively correspond to the number of selected variables in both x and y parts for the r th component.

SPLS [START_REF] Chun | Sparse partial least squares regression for simultaneous dimension reduction and variable selection[END_REF] concentrates on the x part sparsity and needs only 4 parameters to be tuned. Looking at the optimization problem of SPLS [START_REF] Chun | Sparse partial least squares regression for simultaneous dimension reduction and variable selection[END_REF], the sparse solution c, relying on a L1-penalization via the parameter λ1, is not directly sought. Instead the algorithm defines c as a direction relatively close to w, the "associated PLS-based solution". This distance is measured according to the S X,(r) -metric and the importance of each of these two terms in the overall optimization problem is tuned by a parameter κ ∈ [0, 1]. If κ is high (equal 1 at the limit), the problem yields to the PLS [START_REF] Wold | Estimation of principal components and related models by iterative least squares[END_REF] solution, while a low κ means that c and w must be very close to each other but not necessarily to the classical PLS solution. According to [START_REF] Chun | Sparse partial least squares regression for simultaneous dimension reduction and variable selection[END_REF], λ2 can be tuned to +∞ as to preserve the softthresholding solution and κ < 1/2 avoids local solution issues. There are now only 3 parameters (R, κ, λ1) left to be tuned. Algorithm (10) is the associated algorithm in the case of so-called NIPALS deflation, as described by the authors where NIPALS-PLS(.) is a function that takes a covariate and a response matrix for the two first arguments and a number of components in third and sends the NIPALS-PLS regression matrix such as described in Algorithm (7). In the following, the regression matrix is estimated as B = B (R) , also A is the index set of active variables, updated at each iteration, such as r) and S X,(r) = Z (r) Z (r) Initialize c (r) , iterate until convergence of the two steps :

∀r ∈ [[1, R]] :                              Z (r) = X Y (
w (r) = arg min w,w w=1 -κw S X,(r) w + (1 -κ)(c (r) -w) S X,(r) (c (r) -w), c (r) = arg min c Z (r) c -Z (r) w (r) Z (r) c -Z (r) w (r) + λ1|c|, A = i ∈ [[1, p]]|w (r) i = 0 ∨ B (r-1) i,. = 0 , B (r)
A = NIPALS-PLS XA, Y (r) , r and B

(r)

-A = 0, (e') Y (r+1) = Y (r) -XB (r) .

(10)

The deflation in sPLS Lê [START_REF] Lê Cao | A sparse PLS for variable selection when integrating omics data[END_REF] is performed as in the classical NIPALS-PLS algorithm (see row (e) of Algorithm ( 7)), while only the Y part is deflated in the SPLS [START_REF] Chun | Sparse partial least squares regression for simultaneous dimension reduction and variable selection[END_REF] algorithm (using row (e') of Algorithm (10), instead of row (e) of Algorithm ( 7)).

Table 1: Core optimization problem of each PLS-based, where M (r) = Y (r) X (r) /n -1 and S X,(r) = X Y (r) Y (r) X.

Method

Optimization Problem Constraints Parameters NIPALS-PLS [START_REF] Wold | Estimation of principal components and related models by iterative least squares[END_REF]

max u,v v M (r) u u u = v v = 1 R sPLS Lê Cao et al. (2008) min u,v (n -1)M (r) -vu 2 + λ (r) u |u| + λ (r) v |v| u u = v v = 1 R, λ (r) u , λ (r) v r
SPLS [START_REF] Chun | Sparse partial least squares regression for simultaneous dimension reduction and variable selection[END_REF]

min w,c -κwS X,(r) w + (1 -κ)(c -w) S X,(r) (c -w) + λ1|c| + λ2||c|| w w = 1 R, κ, λ1, λ2 ddsPLS arg max u,v v S λr M (r) u u u = v v = 1 R, (λr) r
The method proposed in the next section is a new sparse PLS methodology based on a soft-thresholding operation of the empirical covariance matrix, where the associated hyperparameters are ajusted via bootstrap.

Data-driven sparse PLS

This section details the ddsPLS method. First, Section 4.1 provides the underlying algorithm. Then Section 4.3 introduces the main quality descriptors generally used in PLS regression and their adaptations to the bootstrap context with slight modifications. Section 4.4 explains how to tune the different hyperparameters of the ddsPLS method, and an illustration based on an "easy" example is given in Section 4.6. Finally Section 4.2 provides the extension of the ddsPLS approach to the case of multi-block x datasets.

Remark 3 Since the underling latent variable model ( 2) assumes centered and scaled x and y variables, the data are (marginally) standardized on the overall dataset before applying the ddsPLS methodology. So data are standardized by subtracting the empirical mean and dividing by the empirical standard deviation. This standardization step makes the ddsPLS methodology sensitive to potential extreme observations' adverse effects and thus to outliers.

ddsPLS algorithm

As noticed in [START_REF] Mehmood | Comparison of variable selection methods in partial least squares regression[END_REF]), sPLS Lê Cao et al. (2008) and SPLS Chun and Keleş (2010) are embedded sparse methodologies. But those methodologies are based on the penalization of the x parts from the eigendecomposition of Y'X (and also of the y part for sPLS Lê Cao et al. ( 2008), but through a different regularization coefficient). The idea of the ddsPLS method is to soft-threshold directly the empirical covariance matrix right, using the S λ (.) operator, before building the weights. This operation naturally tends to remove non-interesting variables since their associated soft-thresholded empirical covariance coefficients are equal to 0, for a threshold λ large enough, and the associated weight goes to 0. This regularization of the "data itself" (actually the empirical covariance matrix) and not of the associated weights (which are first eigenvectors of the empirical covariance matrices) led to the name "data-driven sparse PLS (ddsPLS)" to the proposed method. The corresponding optimization problem is provided in the last row of Table 1 . One might notice that if λr = 0, it corresponds to the usual NIPALS-PLS [START_REF] Wold | Estimation of principal components and related models by iterative least squares[END_REF] solution. The associated algorithm of the ddsPLS method is detailed below. Note that, to introduce sparsity in the response in the regression matrix B, the step (d) of the Algorithm ( 7) is modified such as cr = Y (r) Πr tr/t r tr where Πr = diag({δ =0 (vr)j} j∈[ [1,q]] ) is sparse since "vr" is sparse for λ large enough. Here also, if λ = 0, the result is equal to the solution offered by the NIPALS-PLS [START_REF] Wold | Estimation of principal components and related models by iterative least squares[END_REF] solution.

The algorithm of ddsPLS is then

∀r ∈ [[1, R]]                      (a ) ur = --→ RSV S λ (r) M (r) , vr = --→ RSV S λ (r) M (r) , (b) tr = X (r) ur, (c) pr = X (r) tr/t r tr, (d )    Πr = diag({δ =0 (vr)j} j∈[[1,q]] ) cr = arg min V Y (r) Πr -trV 2 = Y (r) Πr tr/(t r tr). (e) X (r+1) = X (r) -trp r , Y (r+1) = Y (r) -trc r , (11)
where the regression matrix estimation is also performed using the relation ( 8).

Notice that R + 1 parameters need to be tuned to obtain an optimal solution: the number R of components and the associated regularization coefficient λ (r) for each component. Based on different well-known quality descriptors used in PLS and on new appealing ones introduced in our ddsPLS context (see next subsection), a way to calibrate these different parameters is provided in Section 4.4.

Adaptation to the multi-block structure

A multi-block generalization of the latent variable model ( 2) is

∀t ∈ [[1, T ]], xt = A t φ + t and y = C φ + ξ, ( 12
)
where T is the total number of blocks. The corresponding underlying regression model can be written as

y = T t=1 B t bt + e. (13)
Let's introduce the following notation:

x = (x 1 , . . . , x T ) , A = [A1, . . . , AT ] and B = [B 1 , . . . , B T ] . (14)
Model ( 12) can be rewritten in the form of Model (2), thus it is possible to use the ddsPLS methodology as previously in this multi-block framework.

Remark 4 The adaptation of the PLS method to the multi-block context has firstly been described in [START_REF] Wold | Three PLS algorithms according to SW[END_REF] and finally defined in [START_REF] Vivien | Approches PLS linéaires et non linéaires pour la modélisation de multi-tableaux[END_REF] to take the form

max w t ,v,c t v Y (r) T t=1 ctX (r) t wt , s.t. ||wt|| = ||v|| = 1 and T t=1 c 2 t = 1. ( 15
)
If u is associated to a solution of any of the previous algorithms (except for SPLS [START_REF] Chun | Sparse partial least squares regression for simultaneous dimension reduction and variable selection[END_REF]) in steps (a), (i † ) or (a ), and ut is the extracted matrix from u associated with block t, then by identification

wt = ut/||ut|| and ct = ||ut||. ( 16
)
Remark 5 For the deflation steps (see steps (e) and (e')), the "mbPLS" algorithm (for multi-block PLS) uses the super-score deflation where the super-score gathers information from the T different blocks:

tr = T t=1 ctX (r)
t wt. Hence, in the following, when NIPALS-PLS [START_REF] Wold | Estimation of principal components and related models by iterative least squares[END_REF] or sPLS Lê Cao et al. ( 2008) are used in the multi-block framework, the super-score deflation is always used. Moreover, for the multi-block version of the ddsPLS algorithm, the super-score deflation strategy is also used.

Note different authors consider that deflating each block on the super-score might mix block information and, therefore, might also drive to senseless conclusions. For example, [START_REF] Westerhuis | Deflation in multiblock PLS[END_REF] present specific designs where no deflation on the X blocks provides better prediction results. This corresponds to the solution chosen in the SPLS [START_REF] Chun | Sparse partial least squares regression for simultaneous dimension reduction and variable selection[END_REF] algorithm.

Remark 6 When the response variable y has multi-blocks, the proposed methodology remains the same where the different blocks are concatenated, as done in this section if x has multi-blocks, see (14), and the ddsPLS methodology can again be applied on the resulting data-set.

Model quality descriptors

The PLS methodology is known to provide the correct common subspaces between the x and the y parts if there is indeed information to catch. Whereas, if this cross structure is poor, the methodology sometimes over-fits focusing on the x part variance in particular when the dimension p of x is large, see for example [START_REF] Wakeling | A test of significance for partial least squares regression[END_REF]; [START_REF] Cloarec | Can we beat over-fitting[END_REF]. Its side effect artificially increases the R 2 statistics, also denoted as explained variance or determination coefficient, which only describes the goodness of fit on the train dataset:

R 2 = 1 - q j=1 n i=1 (yi,j -ŷi,j) 2 q j=1 n i=1 (yi,j -ȳj) 2 , ( 17
)
where ȳj is the empirical mean of the j th -variable of the y part and ŷi,j is the estimation of yi,j thanks to the current model. The numerator is often denoted as the RSS (Residual Sum of Squares) and the denominator as the TSS (Total Sum of Squares) and thus R 2 = 1 -RSS/TSS. This statistics can actually be interpreted as an estimator of

γ(p) = 1 - q j=1 var(yj -y (p) j ) q j=1 var(yj) , (18)
built on the (potentially stochastic) prediction model p where y (p) j is an estimator of yj based on the prediction model p. This metrics is as close to 1 as the prediction error is low. Its interest is to take into account the different scales of the variables yj. As an interpretation, asymptotically (when n → +∞), if this metric is equal to 0, then the associated prediction model performs equivalently to the mean prediction model, and if it is below 0, it performs worse. In the oracle context (a prediction model which knows the true relations between the informational variables), the estimator of y by the oracle predictor model p is y = A φ according to Model (2), and the previous metrics can be written as

γ = 1 - q j=1 var( j ) q j=1 var(yj) . ( 19
)
Remark 7 The two previous statistics are global in the sense of the q-dimensional response variable y and therefore do not allow for a marginal critique of the model. Another interesting statistics would be γj(p) = 1 -var(yj -y (p) j )/var(yj), for j = 1 . . . q, which exhibits the marginal effects of each variable on the response variable. However, the response variables have been standardized in order to be processed on the same scale (see Remark 3), which rationalizes the definition of the γ. Although this does not make it possible to identify marginally the variables of y that can potentially degrade the criterion, their deleterious effects remain observable. Moreover, in the case of marginal γj, it would have been necessary to define a transformation aggregating the q marginal criteria making it possible to define a decision algorithm on the conservation or not of the current component. The definition of such a transformation is not trivial and, therefore, will not be considered in this work. For all these reasons, the global statistics γ is used in the proposed methodology.

The R 2 estimator is subject to over-fitting, this is due to the fact that the quality of the model is assessed on the whole sample used to build the model.

To circumvent this over-fitting drawback, the criterion Q 2 has been introduced by [START_REF] Stone | Cross-validatory choice and assessment of statistical predictions[END_REF]; [START_REF] Geisser | A predictive approach to the random effect model[END_REF]. It is based on F -folds cross-validation such as

Q 2 = 1 - q j=1 F f =1 i / ∈cv f yi,j - ŷ(cv f ) i,j 2 q j=1 F f =1 i / ∈cv f (yi,j -ȳj) 2 , ( 20
) where ŷ(cv f ) i,j
is the estimation of yi,j thanks to the current model where the observations in "cv f " are used as the train sample for fold "f ". The numerator is often denoted as the PRESS (PRediction Error Sum of Squares) and thus Q 2 = 1 -PRESS/TSS. Thanks to R 2 and Q 2 , different rules have been used to select "optimal" models. In the general prediction context, the most popular maximizes Q 2 , which is equivalent to minimize the Mean Squared Error in Prediction (MSEP) since the denominator of the right-hand side term of Q 2 does not use the model's prediction (but the mean estimator). In the context of PLS, a rule of thumb is to favor models for which the difference between R 2 and Q 2 is minimum, which is supposed to correspond to a minimum of over-fitting, see for example [START_REF] Cloarec | Can we beat over-fitting[END_REF].

Those metrics R 2 and Q 2 are helpful to evaluate the quality of a complete model. Other metrics have been introduced to evaluate the quality of a current r th -component of the model and are denoted as R 2 r and Q 2 r , see for example [START_REF] Tenenhaus | La régression PLS: théorie et pratique[END_REF]. Let us suppose that the first (r -1) components have already been built, then the quality of the r th -component can be screened thanks to

R 2 r = 1 - q j=1 n i=1 yi,j - ŷ(r) i,j - ŷ(r-1) i,j -ȳj 2 q j=1 n i=1 (yi,j -ȳj) 2 , Q 2 r = 1 - q j=1 F f =1 i / ∈cv f yi,j - ŷ(cv f ,r) i,j 2 q j=1 F f =1 i / ∈cv f yi,j - ŷ(cv f ,r-1) i,j 2 ,
(21) where ŷ(r) i,j , resp. ŷ ,(r) i,j , is the prediction of the i th -observation and the j th -response variable thanks to model based on components the first r and all observations, resp. all observations except the i th one. By convention, ŷ(0) i,j = ȳj and ŷ

,(0) i,j = 1 n-1 i =i y i ,j . The metrics R 2
r represents the proportion of variability which is explained by the new component r and not by the (r -1) previous ones, while the Q 2 r does the same but based on out-of-fold samples from the cross-validation fold. Rules of thumb exist to validate the current component based on the value of Q 2 r , for example Q 2 r ≥ 0.0975 following for instance Tenenhaus (2005); Lê [START_REF] Lê Cao | A sparse PLS for variable selection when integrating omics data[END_REF], or Q 2 r ≥ 0.05 according to Pérez-Enciso and Tenenhaus (2003); Tenenhaus (2004)) or Q 2 r ≥ 0 (or ≥ 0.05 if n ≤ 100) according to [START_REF] Thévenot | ropls: Pca, pls (-da) and opls (-da) for multivariate analysis and feature selection of omics data[END_REF]. Based on theses descriptors Q 2 r , another metrics which describes the overall quality of a R-components model is the cumulative Q 2 criterion, denoted as Q 2 R,(cum) hereafter, and defined as an aggregation of the R criteria

(Q 2 1 , . . . , Q 2 R ) such as Q 2 R,(cum) = 1 - R r=1 1 -Q 2 r . This metric Q 2 R,(cum) is used in Tenenhaus (1998) with the rule that the current component R is retained if Q 2 R,(cum) >> Q 2 R-1,(cum) .
In the previous statistics Q 2 r (and thus in the metric Q 2 R,(cum) ), the performances of numerous prediction models are simultaneously mixed, over the F folds, in the numerator and in the denominator. One consequence of this is that model selection is driven by prediction models for which the errors are the most important. To circumvent this drawback, for F -folds cross-validation, one way is first to consider the F different "marginal" Q 2 -like statistics, denoted as 1 -

q j=1 i / ∈cv f yi,j - ŷ(cv f ) i,j 2
q j=1 i / ∈cv f yi,j - ȳ(cv f ) j 2 , for f = 1, . . . , F, (22)
and then to combine them a posteriori using the empirical mean. To the best of our knowledge, this approach has never been proposed in the PLS literature. However, when the sample size n is small, using F -folds cross-validation is not relevant and the use of bootstrap is interesting to enrich the underlying information of the available data. This is the strategy chosen for the continuation of this work.

Based on B bootstrap samples, the bootstrapped versions of R 2 and Q 2 are given by

R2 B = 1 B B b=1 R 2 b and Q2 B = 1 B B b=1 Q 2 b (23)
with, for the current bootstrap sample b,

R 2 b = 1 - q j=1 i∈IN(b) yi,j -ŷb i,j 2 q j=1 i∈IN(b) yi,j -ȳb j 2 , Q 2 b = 1 - q j=1 i∈OOB(b) yi,j -ŷb i,j 2 q j=1 i∈OOB(b) yi,j -ȳb j 2 ,
(24) where IN(b) and OOB(b) are respectively the set of the in-bag observations and the set of the out-of-bag observations, ŷb i,j is the prediction of the i th -observation and the j th -response variable and ȳb j is the empirical mean of the j th -response variable estimated on the bootstrap sample b.

In the same way, bootstrapped versions of R 2 r and Q 2 r are given by

R2 B,r = 1 B B b=1 R 2 b,r and Q2 B,r = 1 B B b=1 Q 2 b,r (25)
where

R 2 b,r = 1 - q j=1 i∈IN(b) yi,j - ŷb,(r) i,j - ŷb,(r-1) i,j -ȳb j 2 q j=1 i∈IN(b) yi,j -ȳb j 2 , Q 2 b,r = 1 - q j=1 i∈OOB(b) yi,j - ŷb,(r) i,j 2 q j=1 i∈OOB(b) yi,j - ŷb,(r-1) i,j 2 . ( 26
)
Remark 8 Those metrics, R2 B , Q2 B , R2 B,r and Q2 B,r , share the same philosophies as those discussed for (17), ( 20) and ( 21), but they introduce two innovative ideas (except for the bootstrap analysis itself ). Firstly, they are empirical means over all bootstrap replications. This allows not mixing bootstrap quadratic errors, which is a good feature since estimations can vary from one bootstrap sample to another. Secondly, the "R 2 "-spirit metrics are based directly on the b th bootstrap replication but use the associated in-bag sample, while the "Q 2 "-spirit metrics use the associated out-of-bag sample. This implies that the estimators R 2 b and R 2 b,r are even more subject to over-fitting and lead the difference " R2

B -Q2 B " to be more discriminant and informative in the presence of over-fitting case.

The following section details how these metrics are used to adjust the hyperparameters of the ddsPLS method by following the spirit of their cross-validation versions.

Selection of the hyperparameters

The ddsPLS algorithm detailed in (11) and the other sparse PLS methods require hyperparameters to be adjusted. In the case of ddsPLS, the following regularization parameters, R and ΛR = {λr, r = 1, . . . , R}, must be optimally tuned. Due to the variance constraints of the variables x and y, we know that λr ∈ [0, 1] where

• λr = 0, corresponds to the NIPALS-PLS solution,

• λr = 1, leads to no variable selected and thus to the prediction model equal to empirical mean estimation.

This optimal hyperparameter selection is based on the minimization of R 2 r -Q 2 r (that might correspond to a minimum in over-fitting for the r th -component). To this end, the proposed ddsPLS strategy considers R2 B,r -Q2 B,r , the bootstrapped version of the corresponding metrics defined in (23). More precisely, the selected model must satisfy that, for a current model built on (r -1) components with estimated values for regularization parameters Λr-1 = { λ1, • • • , λr-1}, the r th -component is acceptable if 1. the r th -component performs better than the mean estimation (i.e. Q2

B,r > 0), 2. the model built on r components performs better than the one built on (r -1) components (i.e. Q2 B is increasing with the number of components), 3. the bootstrapped explained variance is close to the bootstrapped cross-validated explained variance (i.e. R2 B -Q2 B is minimum on the set of hyperparameters). This hyperparameter selection procedure is illustrated on a toy example in Section 4.6.

A theoretical lower bound for the regularisation coefficients

Among the sparse covariance estimators the solution developed by [START_REF] Cai | Adaptive thresholding for sparse covariance matrix estimation[END_REF] uses a fully datadriven estimator for each element of the covariance matrix. Even if they consider the variance case (where y = x) and assume that random parts are normally distributed, they use the adaptive threshold

λi,j = δ θi,j log p n , (27)
where δ can be whether taken equal to 2, whether empirically chosen though cross-validation and θi,j = var ((xi -Exi)(xj -Exj)), the variance of the covariance coefficient between xi and xj.

As n decreases/p increases, the probability for a covariance coefficient to be larger than a given value increases. The value of λi,j estimates the variability of the covariance coefficient below which non null coefficient can be considered as null with a high probability.

Inspired from the previous work, let us introduce the following value

λ (0) r = 1 pq q j=1 p i=1 θ (r) (j,i) log (max(p, q)) n where θ (r) (j,i) = 1 n n k=1 x (r) (k,i) y (r) (k,j) -m (r) (j,i) 2 , where x (r) (k,i) , y (r)
(k,j) and m (r) (j,i) are the elements of the matrices X (r) , Y (r) and M (r) respectively at positions (k, i), (k, j) and (j, i) respectively. This value λ (0) r is the lowest accessible value by the algorithm detailed in the previous section. It prevents from building components in which a large number of variables would be selected. This lower bound for λ is implemented in the package.

A toy example

Consider the latent variable model (2) with the following "easy" data structure

A = √ 1 -σ 2 1 50 0 950 , D = √ 1 -σ 2 1 , (28)
where 1 -σ 2 = 0.9025 and ψ = (φ , 1...50 /σ, 51...1000 , ξ/σ) ∼ N (0, I1+1000+1) .

The dependent variable y is only associated with the first 50 variables of x while the last 950 variables of x are not linked to y. According to Remark 1, the true value for the number of components is R = 1. Thus, in this simulation framework, one of the objectives of the ddsPLS method is to properly retain one component for which the first 50 variables of x are selected in the final model. Moreover the corresponding prediction squared error is expected around "γ = 1 -σ 2 = 0.9025", for a sample with a reasonable number n of observations.

and Q2

B,1 versus λ, for each sample size. The "optimal" λ1's (for which R2 B,1 -Q2 B,1 are minimum) are also located in these graphs. One might notice that selecting model considering only Q2

B,1 would lead to bad conclusions, where its maximum (see the colored star on the curve of Q2 B,1 ) is around λ = 0.8 (for n = 50 or 200) or 0.3 (for n = 100) which is an over-fitting region according to the difference R2

B,1 -Q2 B,1 . The graph at the bottom left of Figure 1shows the number of selected variables for each value of λ and each of the 3 samples (n = 50, 100 and 200). These curves are decreasing until a first loss of the 950 noninformative variables of x which occurs naturally at lower λ values when n is large. Then one can observe a large plate, and finally, the curves become decreasing again, resulting in the loss of the 50 relevant variables of x in the model. Indeed, as n increases, the minimum correlation to remove the 950 non-informative variables of x is small, and the maximum correlation keeping the 50 informative variables is high. For the 3 selected models (with the optimal λ1 values), only the 50 relevant variables of x are selected. The right bottom graph of Figure 1shows boxplots of the corresponding first 50 regression coefficients for the 3 sample sizes. As the variables share the same distribution and are independent, the regression coefficient associated with the first 50 variables of x should naturally be close to 1/50, while the last 950 coefficients regression coefficient associated with the last 950 variables of x are equal to 0. One might notice from the three boxplots, the coefficients of the first 50 relevant variables get closer to this expected value 1/50 as n increases.

In the following, let's define the True Positive Rate (TPR) and the False Positive Rate (FPR) of variable (of x) selection as TPR = TP/P and FPR = FP/N where TP is the number of selected variables among those that should be selected, P is the total number of variables that should be selected, FP is the number of selected variables among those that should not be selected, and N is the total number of variables that should not be selected. Naturally, the optimal target is to have TPR = 1 and FPR = 0. In this simulation study, the TPR and FPR of variables of x selection are optimal for the 3 sample sizes. Note that the TPR and FPR are also optimal for the variable of y selection since the one-dimensional dependant variable is always retained in the final selected models. Finally, according to Equation ( 6), the metrics ||A B -D|| 2 /||D|| 2 can be calculated in this simulation framework: its value is approximately equal to 0.002 < < 1 and so shows reliable results. This "toy" example illustrated the principle of the selection criterion for the ddsPLS method. In Section 5, others simulations will lead to estimate its good numerical behavior in more complex contexts.

Simulation study

Different synthetic structures have been used to evaluate the numerical behavior of the proposed ddsPLS methodology. Since the multi-block framework can be reduced to a mono-block framework, only mono-block structure in the x part are considered in this simulation study via three specific designs.

The first design corresponds to two groups of variables in the x part, which are not correlated, and a 3-dimensional response variable y. The second design simulates spectroscopic data and will be detailed below.

For each design, 100 replications will be considered.

In this simulation study, various methodologies are compared with the proposed "ddsPLS" approach. The different methodologies are based on what has already been proposed by other authors and adapted to the bootstrap versions of γ. More precisely,

• "sPLS 1" , "sPLS 2" and "sPLS 3" build sparse PLS models as detailed in Lê [START_REF] Lê Cao | A sparse PLS for variable selection when integrating omics data[END_REF]. The additional notation "1" refers to the solution minimizing the difference R2 B,r -Q2 B,r for each component r. The notation "2" refers to the solution maximizing Q2 B,r for each component r, which is the closest solution to what is done currently using this model, corresponding to minimizing the cross-validation error (here replaced with maximizing the out-of-bag normalized error). Finally the notation "3" refers to model which builds variables while there exists (keepX , keepY ) such as Q2 R,(cum) > 0.0975 (based on CV results) where the couple (keepX , keepY )r is then choosing maximizing Q2 R,(cum) . This is what the most popular methodology for this method, see for example [START_REF] Tenenhaus | La régression PLS: théorie et pratique[END_REF].

• "NIPALS-PLS" corresponds to classical PLS algorithm where the number of components is chosen minimizing R2 B,r -Q2 B,r for each component. • "SPLS" builds sparse PLS models as detailed in [START_REF] Chun | Sparse partial least squares regression for simultaneous dimension reduction and variable selection[END_REF]. The parameter κ is left to the value "κ = 0.5" (knowing that κ ∈ ]]0, 1]]) which provides a balanced compromise between sparsity (low κ) and prediction results close to the classical PLS problem (high κ). Using a grid of values corresponding to η = {0.1, 0.2, . . . , 0.9} and R = {1, 2, . . . , 10}, the considered best model corresponds to the minimum of the mean squared prediction error.

For "SPLS" and "sPLS 3" , the maximum number of possible components is fixed to 10. Also, for those two methodologies, the metric Q 2 , based on the F -folds cross-validation, is used, with a number F of folds equal to F = (25, 50, 100, 100, 100) for n = (25,50,100,200,400) respectively. Note that, when n = (25, 50, 100), leave-one-out cross-validation is performed.

For the other approaches ("sPLS 1" , "sPLS 2" , "NIPALS-PLS" and "ddsPLS" ), Bootstrap versions of Q 2 or R 2 metrics are used with a number B of bootstrap samples adapted to the sample size n: more specifically B = (1000, 500, 300, 100, 100) for n = (25,50,100,200,400) respectively.

Design 1: simulated data

The first design considers The corresponding latent variable model is then:

A = √ 1 -σ 2 α313,
xj =    √ 1 -σ 2 (φ1 + φ2 + φ3)/ √ 3 + j for j = 1 . . . 50 √ 1 -σ 2 (φ4 + φ5)/ √ 2 + j for j = 51 . . . 100 j for j = 101 . . . 1000 and    y1 = √ 1 -σ 2 (φ1 + φ2 + φ3)/ √ 3 + ξ1 y2 = √ 1 -σ 2 (φ4 + φ5)/ √ 2 + ξ2 y3 = ξ3
(30) An analysis of matrices A and D described in (29) shows that R = 5, but there are multi-colinearities in A and the real objective for a relevant number of components might be R = 2, which is proved according to [START_REF] Helland | Comparison of prediction methods when only a few components are relevant[END_REF], as noticed in Remark 1. Also, only 100 variables over the 1000 variables of the variable x should be selected, moreover 2 out of the 3 variables of the dependent variable y should be selected. The theoretical γ is equal to γ * = 2(1 -σ 2 )/3 = 0.6534 according to (19).

According to Figure 2and for large n (=200 or 400), all the methodologies almost always build the correct number of components, with the exception of the "SPLS" approach. For low n (=25, 50 or 100), the "NIPALS-PLS" and "sPLS 3" methodologies mostly build far too many components, which is also the case for the "sPLS 2" and "SPLS" methodologies but in a mitigate manner with median values equal to 3. The "ddsPLS" methodology and "sPLS 1" almost always build the correct number of components with a clear advantage for the first-mentioned method.

According to Figure 3 , when n = 400 (large sample size), all the methodologies correctly provide "Q 2 " values close to γ * . For all the cases, "ddsPLS" underestimates the value of Q 2 . When n = 25, the bootstrap based versions of Q 2 seem to under-estimate γ * while the "sPLS 3" over-estimates it and "SPLS" shows high variability around this value. For "NIPALS-PLS" , "sPLS 3" and "SPLS" methodologies, the bad performances are due to over-fitting because of the high values retained for R for these approaches. However, the proposed Q 2 metric, Q2

B , used for "NIPALS-PLS" (and for "sPLS 1" and "sPLS 2" ), detects overfitting providing low Q 2 values, while "sPLS 3" provides high values of Q 2 which are often larger than the true expected value, this being essentially due to how Q 2 is calculated and not representative of the actual quality of the model. Finally, the "sPLS 1" and "sPLS 2" methodologies show slightly better results than "ddsPLS" in that low n context. Even if this metric seems worse for "ddsPLS" , this might be subject to over-fitting and other descriptors (such as ||A B -D|| 2 /||D|| 2 , or a prediction error on an independent data-set) should be used. This results have been observed and are accessible on the right part of the Figure 3 and in the Figure 4 . Those figures are commented below. In Figure 3 , from the boxplots of ||A B -D|| 2 /||D|| 2 , for each methodology and each sample size, the methodologies ranked in descending order of performance are: "ddsPLS" , "sPLS 3" , "SPLS" , "sPLS 1" , "sPLS 2" and "NIPALS-PLS" for low n. As n increases, this order is maintained and all the methodologies get better results for this metric.

We define the scaled RMSE (Root Mean-Squared Error) such as

∀j ∈ [[1, q]], scRMSEj = n test i=1 (yi,j -ŷi,j) 2 n test i=1 (yi,j -ȳj) 2 ,
estimated on a independent sample of size ntest where ȳj is estimated on the independent test sample and ŷi,j is estimated according to the current prediction model. The Figure 4shows the scaled RMSE for the different methodologies in all cases, for each response variable separately. Methodologies "ddsPLS" and "sPLS 3" show the best performances for low y1 and y2 while "ddsPLS" and "NIPALS-PLS" are the best for y3. More interestingly, "sPLS 3" and "SPLS" are especially not precise for y3, this is associated with the remark formulated in the previous paragraph associated with the potential over-fitting of those methodologies.

In terms of variable selection, the Figure 5insures that the True Positive Rate (TPR) of variable of x selection is mostly equal to 1 for all methodologies, except for "ddsPLS" . On the contrary, the False Positive Rate (FPR) is mostly equal to 0 only for "ddsPLS" and is closely equal to 0 for "SPLS" . Other methodologies based on "sPLS" decrease their FPR with n until n reaches 100. More precisely, the FPR of "sPLS 1" indeed slowly decreases with n, "sPLS 2" and "sPLS 3" increase their FPR for n ≥ 100. This shows that Q 2 maximization based methods (both for cross-validation and bootstrap versions) seem to build sparser models as the number n of individuals grows. On the contrary, Q 2 -R 2 minimization based methods decrease their FPR as n grows.

In the y part, regardless of n, all methodologies select always y1 and y2, and so the FPR equal 1 in all cases.

According to Figure 5 (b), all methodologies naturally decreases the selection rate of y3 as n increases, except for "sPLS 3" . Furthermore, among the bootstrapped methodologies, "ddsPLS" outperforms the others.

In terms of variable selection, this can lead to the conclusion that bootstrap-based versions provide better results than those based on cross-validation, and more precisely the one associated with the minimization of R2 B,r -Q2 B,r which perform better than those based on the maximization of Q2 B,r , for both x and y parts. Design 1 with varying q As to evaluate the performance of the proposed methodology if q varies, it has been considered to add extra response variables, as to reach q = {10, 50, 100}. Each of these extra response variables is not associated nor with the other response variables nor with the other response variables. As previously, N = 100 different data sets (with n = 50) have been generated to build different prediction models. Different information can be highlighted.

• All the models are built on 2 components.

• All the models select, in the response part, only the two relevant variables (the two first ones).

• The first row of Figure 6gives the scaled RMSE of the N models, for each scenario and for each response variable. An independent data set of size ntest = 1000 has been used to estimate those distributions.

It is clear that all the models predict well the two selected variables (y1 and y2) and do not "predict" for others but returns the empirical mean values.

• The Second row of Figure 6shows the frequency of selection in the covariate part for each scenario.

Only the 100 first covariates, associated with y1 and y2, are selected.

This experiment shows that the proposed "ddsPLS" methodology is robust to an increasing number q of response variables.

Design 2: simulated data inspired by real data

In this design, datasets with structures similar to spectroscopic data have been generated. Using the routine already used by [START_REF] Cloarec | Can we beat over-fitting[END_REF], the x part contains spectra observations and is divided into two blocks. The first one (X1) corresponds to spectroscopy with large overlapping peaks and a low number of variables, such as Ultraviolet, Visible and Near-Infrared spectrophotometry. The second block (X2) has narrow peaks and a high number of variables but with much less overlap between peaks. The blocks Xt, t = 1, 2 are built using the following procedure:

Xt = 1 -σ 2 t CtSt + σtEt (31)
with St = (st,r,j) (r,j) ∈ R R×p t , Ct = (ct,i,r) (i,r) ∈ R n×R and Et = e t,i i ∈ R n×p t where

∀(r, j) ∈ [[1, R]] × [[1, pt]], st,r,j = Kr k=1 h r,k exp -(j -π r,k ) 2 / 2ω r,k 2
, et,i ∼ N (0p t , Ip t ) and σt = 0.05.

The score matrix is here denoted Ct for "Concentration", and the matrix of the weights is denoted St for "Spectrum". Let pt be the number of variables in each block t and Rt the number of latent variables (Rt = 10 for both blocks). For each latent variable, r ∈ [ [START_REF] Liquet | Group and sparse group partial least square approaches applied in genomics context[END_REF]Rt]], the number Kr of peaks is selected randomly from a uniform discrete random distribution U[ [START_REF] Liquet | Group and sparse group partial least square approaches applied in genomics context[END_REF]5]] in order to generate, ∀k ∈ [ [START_REF] Liquet | Group and sparse group partial least square approaches applied in genomics context[END_REF]Kr]], a spectrum with the following characteristics:

• position of the peaks is sampled from a discrete uniform distribution π r,k ∼ U [[2, pt]];

• relative peak heights are sampled from a continuous uniform distribution h r,k ∼ U (0, 1);

• peak widths are sampled from a continuous uniform distribution ω r,k ∼ U (2, 3) and are proportional to the number of variables.

The previous parameters allow to build the matrix St. To build the matrix Ct such as the intensity of each latent variable in each observation, the following Gaussian distribution is used: ct,i,r ∼ N (µt,r, σ 2 t,r ) where • the expected intensity of each latent variable is randomly sampled from an exponential distribution: log (µt,r) ∼ U (0, 1), (b) Variable selection performance in the y part measured by the number of selected y variables where the gray areas correspond to theoretical objectives, that is two selected variables: y 1 and y 2 (which are always selected, providing a TPR equal to 1). • the variance of the intensity for each latent variable is randomly sampled from a uniform distribution σ 2 t,r ∼ U(0, σµµt,r) where σµ, when low enough, allows to generate positive values for ct,i,r with high probability.

Number of select variables

Moreover, in order to get the correlation between the two blocks X1 and X2, the first five columns of C2 have been replaced by the first five columns of C1. Figure 7provides a visualization of an extract from the x part for a dataset simulated according to the procedure.

The 5-dimensional y part is generated from the model Y = 1 -σ 2 y Cy + σyEy, where each row of Ey follows N (0, I), σy = 0.2 corresponds to the standard deviation of this observation noise, and the matrix Cy is used to manage the link with the x part. More precisely, the first 3 columns of Cy are the standardized first 3 columns of C1. The last two columns of Y are generated from a standardized normal distribution. Hence only the first 3 variables of the y part are related to the x part.

According to the way the spectra are generated, the objective number R of components of the underlying model is equal to 3 with a high probability. Figures 8 (a) and 8 (b) respectively show, for 100 replications, the number R of components retained in each model and the associated number of selected variables in the x part. Since the correct dimension is equal to R = 3, all the methodologies (apart from "NIPALS-PLS" , which builds far too many components) succeed to get most of times models with three components, and only "ddsPLS" has a clear majority of models built on three components. Moreover, according to Figure 8 (b), "sPLS 1" and "sPLS 2" build often non-sparse models with very large numbers of selected variables. "sPLS 3" also builds models with various numbers of selected variables and also builds non sparse models. The methodologies "ddsPLS" and "SPLS" build the sparsest models, with a clear advantage for the first methodology even if "ddsPLS" sometimes builds non sparse models. Note that these graphs do not make it possible to appreciate the quality of the selection procedure and the corresponding estimated model different methodologies; this point is discussed below. Since the simulated spectra are not sparse, it is not possible to directly evaluate the relevance of the variable selections. However, the 3 spectra, which are common to X1, X2 and Y are known in this simulation framework. It is reasonable to assume that the relevance of a variable is associated with the intensity of this variable in these three spectra of interest. Then, the selected variables are described in terms of TPR and FPR by varying a threshold α over the maximum values of the spectra for each variable along with the three spectra of interest as to build ROC (Receiver Operating Characteristic) curves and so create a common score of selectivity corresponding to the Area Under the Curve (AUC). For a relevant analysis, a common distribution of α for all methodologies has been used.

The top right graph of the Figure 9shows the AUC metrics computed for the models built over the 100 replications of datasets. It allows evaluating the performance of the methodologies in terms of variable selection. It seems that "sPLS 1" and "sPLS 2" are the least efficient methodologies to select the correct variables efficiently. The methodologies "sPLS 3" and "SPLS" get better AUC but keep high variability. Finally, the proposed "ddsPLS" approach provides the best AUC.

Top left graphs (density plots) of the Figure 9detail more precisely the quality of the selection procedure of the five competing methodologies (in columns), via the TPR versus α (first row), the FPR versus α (second row) and the FPR versus the TPR (third row).

• From the first row, one might notice that "sPLS 1" and "sPLS 2" present a huge loss of TPR for very low α and do not succeed to keep TPR close to 1 for low values of α. "sPLS 3" also suffers from a fall although less important and even succeeds to keep TPR close to 1 for low values of α. "SPLS" keeps this tendency, and "ddsPLS" performs even better. • According to the second row, it seems that "sPLS 1" and "sPLS 2" quickly (in terms of α) fall in FPR, which is encouraging in a certain way, while the three other methodologies keep it high for low values of α. This actually, quite paradoxically, allows considering that "ddsPLS" , "SPLS" and "sPLS 3" are, in fact, better methodologies because they do not tend to remove too quickly variables; feature that is often criticized in sparse methods. • The third row shows the ROC curves in density and once again confirms the previous observations:

"ddsPLS" , "SPLS" and "sPLS 3" , in decreasing order of performance, outperform the other methodologies since the high densities are concentrated at the top of their respective plot, corresponding to high TPR. Unfortunately, the precision is not sufficient to allow a correct differentiation of each distribution for high TPR's.

The five histograms in Figure 9 (bottom right graphs) show the conditional distribution of the FPR given TPR > 0.9. "ddsPLS" clearly provides the heaviest distribution in that TPR area and gets the heaviest lower distribution tail, which are the two reasonable objectives of a selection method. Figure 9: Variable selection performances for design 2. Top left graphs: density plots for the five competing methodologies (in colums) of the TPR versus α (first row), of the FPR versus α (second row) and of FPR versus TPR (third row), red areas correspond to densities higher than 4000 for the first two rows and 3000 otherwise. Top right graph: boxplots of AUC for the five competing methodologies. Bottom right graphs: histograms of FPR given TPR¿0.9 for the five competing methodologies.

In terms of Q 2 results (not provided here), "ddsPLS" , "SPLS" , "sPLS 1" and "sPLS 2" are the most efficient methodologies (respectively equal to 0.47, 0.46, 0.47 and 0.43 in median for the all simulated datasets) for R = 3 components built models. Note that those four methodologies get lower Q 2 than expected (γ = 3 2 y /5 = 0.597 in the considered simulation framework). Only "sPLS 3" gets higher median for Q 2 , equal to 0.55. For more complex models built with R > 3 components, the results also show that "sPLS 3" give a median Q 2 equal to 0.64 for R = 4, 0.75 for R = 5, 0.93 for R = 6, 0.91 for R = 8 and 0.99 for R = 10, while the four other methodologies get decreasing values for median Q 2 . Finally, in the context of that simulation design, it seems again that "sPLS 3" does not succeed to beat over-fitting.

Concerning the selection in the y part (which is only available for "ddsPLS" , "sPLS 1" , "sPLS 2" and "sPLS 3" methodologies), the following results were obtained. Below, the notation (yj, Nj) means that variable yj has been selected Nj times over the 100 replications:

• "ddsPLS" : (y1, 87), (y2, 87), (y3, 81), (y4, 18), (y5, 19),

• "sPLS 1" : (y1, 83), (y2, 86), (y3, 74), (y4, 37), (y5, 36),

• "sPLS 2" : (y1, 88), (y2, 84), (y3, 82), (y4, 53), (y5, 51),

• "sPLS 3" : (y1, 89), (y2, 91), (y3, 81), (y4, 32), (y5, 33).

It is clear that "ddsPLS" outperforms other methodologies regarding the FPR's in y part selection, while the TPR's appear to be equivalent for the four methodologies.

Conclusion

The proposed ddsPLS methodology is a sparse PLS method, both in x and y, based on soft-thresholding using bootstrap operations for model selection. This solution fixes sparsity in both x and y parts using a single parameter, which is also interpretable in terms of minimal empirical covariance allowed to be selected in the model. Moreover this solution shows very good results both in prediction and in selection, through the various simulation models analysed. While most of the other sparse PLS are based on cross-validation approaches, our solution uses bootstrap. This allows building smooth curves of the chosen criterion even in the ill-conditioned contexts of low n and large p data. Based on the proposed bootstrap version of sparse PLS, a new relevant definition of the Q 2 and the R 2 statistics denoted as Q2

B and R2 B where B is the number of bootstrap samples, has been introduced. Those new versions of Q 2 and R 2 do not mix the errors of the different folds but aggregate the fold estimations through the bootstrap mean estimator. Different aspects must be highlighted regarding those new statistics. The newly introduced γ is a theoretical limit for the Q 2 and R 2 statistics that our criterion respects while the tested competing methods violate it in the n << p context even for simple simulation schemes. The chosen criterion R2

B -Q2 B allows to efficiently select sparse models. Finally, the ddsPLS methodology can be applied to different data contexts such as biological data, omics data and chemometrics, where the sample size n is indeed much lower than the number p of variables.

An implementation of the proposed ddsPLS methodology is available at https://github.com/hlorenzo/ ddsPLS2. A vignette is also available to facilitate the use of the corresponding R package and the associated functions for any user not necessarily specialised in R or statistical modelling. Note that this package uses the foreach (v.1.5.0) package to parallelize the computations associated with the different bootstrap samples. These codes are written in C ++ thanks to the packages Rcpp (v.1.0.5) and RcppEigen (v.0.3.3.7.0) in order to reduce the computation time and to ensure a portability towards other languages, Python in particular.

Whatsmore, all simulation codes (regarding Toy Example, Design 1 and Design 2) are accessible via the repository https://github.com/hlorenzo/simulation_ddspls, where two files (README.md and simulation ddspls.html) help to use the different simulation functions.

Figure

  Figure 1gathers results of 3 simulations based on 3 sample sizes, respectively, n = 50, 100, 200. For the bootstrap-based selection criterion, B = 50 bootstrap samples are used. For each of these samples, the proposed algorithm keeps only one component. The 3 graphs at the top of Figure 1provide the plot of R2 B,1

Figure 1 :

 1 Figure 1: Simulation results for the toy example.

Figure 2 :

 2 Figure 2: Number of components built for design 1, gray areas correspond to the objective, that is R = 2 components.

Figure 3 :

 3 Figure 3: Main performance metrics for design 1. The horizontal line represents the value of Q 2 = γ = 2ε 2 /3 ≈ 0.6534 which corresponds to the objective, where ε = √ 1 -σ 2 .

Figure 4 :

 4 Figure 4: Scaled RMSE (Root Mean Squared Errors) for design 1, for each response variable, on an independent test sample of size n test = 1000. Horizontal lines represent the objective values ≈ (0.141, 0.141, 1).

Figure 5 :Figure 6 :

 56 Figure 5: Variable selection performances in the x and y parts for design 1

Figure 7 :

 7 Figure 7: First five observations of simulated datasets X 1 and X 2 from design 2.

  Number R of components built for each methodology for design 2. Gray area corresponds to the correct value R = 3. Number of selected variables in the x part for each methodology for design 2.

Figure 8 :

 8 Figure 8: Estimated number of components R (top) and number of selected variables (bottom) for design 2.

Acknowledgments

We thank the editor and the two anonymous referees for their comments that led to an improved version of this paper. Hadrien Lorenzo received the support of "Sartorius, Corporate Research Advanced Data Analytics" for this project.