Data-Driven Sparse PLS (ddsPLS)

Construction and prediction errors

The ddsPLS methodology uses the principle of construction and prediction errors. The first one, denoted as \(R^2\), evaluates the precision of the model on the training data-set. The second one, denoted as \(Q^2\), evaluates the precision of the model on a test data-set, independent from the train data-set. The \(R^2\) is well known to be sensible to over-fitting and an accepted rule of thumb, among PLS users, is to select model for which the difference \(R^2-Q^2\) is minimum.

The ddsPLS methodology is based on bootstrap versions of the \(R^2\) and the \(Q^2\) defined below.

Bootstrap versions of the \(R^2\) and the \(Q^2\)

More precisely, for the bootstrap sample of index \(b\), among a total of \(B\) bootstrap samples, the \(R^2_b\) and the \(Q^2_b\) are defined as \[ \begin{array}{cc} \begin{array}{rccc} &R^{2}_b & = & 1-\dfrac{ \left|\left|\mathbf{y}_{\text{IN}(b)} -\hat{\mathbf{y}}_{\text{IN}(b)}\right|\right|^2 }{ \left|\left|\mathbf{y}_{\text{IN}(b)} -\bar{\mathbf{y}}_{\text{IN}(b)}\right|\right|^2 }, \end{array} & \begin{array}{rccc} &Q^{2}_b & = & 1-\dfrac{ \left|\left|\mathbf{y}_{\text{OOB}(b)} -\hat{\mathbf{y}}_{\text{OOB}(b)}\right|\right|^2 }{ \left|\left|\mathbf{y}_{\text{OOB}(b)} -\bar{\mathbf{y}}_{\text{IN}(b)}\right|\right|^2 }, \end{array} \end{array} \] where the IN(b) (IN for “In Bag”) is the list on indices of the observations selected in the bootstrap sample \(b\) and OOB(b) (OOB for “Out-Of-Bag”) is the list on indices not selected in the bootstrap sample. Then, the subscripts IN(b) correspond to values of the object taken for in-bag indices (respectively for subscript notation OOB(b) and out-of-bag indices). Also, the “\(\hat{\mathbf{y}}\)” notation corresponds to the estimation of “\(\mathbf{y}\)” by the current model (based on the in-bag sample) and “\(\bar{\mathbf{y}}\)” stands for the mean estimator of “\(\mathbf{y}\)”. The ddsPLS methodology aggregates the \(B\) descriptors \(R^2_b\) and \(Q^2_b\) as follows: \[ \bar{R}^{2}_{B}=\frac{1}{B}\sum_{b=1}^{B}R^2_b ~~~\mbox{and}~~~ \bar{Q}^{2}_{B}=\frac{1}{B}\sum_{b=1}^{B}Q^2_b. \]

Even if the notations of the metrics \(\bar{R}^{2}_{B}\) and \(\bar{Q}^{2}_{B}\) show a square, they can be negative. Indeed, for large samples, they actually compare the quality of the built model to the mean prediction model:

• if the metrics is \(>0\) then the model works better than the mean estimator,
• if the metrics is \(<0\) then the model works worse than the mean estimator.

As the objective of a linear prediction model is to build models better than the mean prediction model, a rule of thumb is to select models for which \[ \begin{array}{ccc} \mbox{(A)} &:& \bar{Q}^{2}_{B}>0. \end{array} \]

A second rule of thumb is necessary to determine if a new component is relevant to improve the overall prediction power of the model. If the \(r^{th}\)-component is tested, then the condition writes

\[ \begin{array}{ccc} \mbox{(Ar)} &:& \bar{Q}^{2}_{B,r}>0, \end{array} \]

where \(\bar{Q}^{2}_{B,r}\) is the “component-version” of \(\bar{Q}^{2}_{B}\)” defined in Appendix A.1.

Remark 1 The metrics \(\bar{R}^{2}_{B}\) and \(\bar{Q}^{2}_{B}\) are in fact estimators of a same statistic \(\gamma\), associated with a prediction model \(\mathcal{P}\): \[ \begin{array}{ccc} \gamma(\mathcal{P}) & = & 1-\dfrac{ \sum_{j=1}^q \mbox{var} ({y}_{j} -y_{j}^{(\mathcal{P})}) }{ \sum_{j=1}^q \mbox{var} ({y}_{j})}, \end{array} \] where \(y_{j}^{(\mathcal{P})})\) is the prediction of the \(j^{th}\) component of \(\mathbf{y}\) by the model \(\mathcal{P}\). The closer this statistic is to 1, the more accurate the model \(\mathcal{P}\) is.

Automatic tunning of the parameters

As said before, a ddsPLS model of \(r\) components is based on \(r\) parameters \((\lambda_1,\dots, \lambda_r)\), denoted \(\mathcal{P}_{{\lambda}_1,\dots,{\lambda}_{r}}\). The metrics \(\bar{R}^{2}_{B}\), \(\bar{Q}^{2}_{B}\) or \(\bar{Q}^{2}_{B,r}\) are functions of the values taken by the estimated ddsPLS models \(\widehat{\mathcal{P}}\), which depend on the \(r\) estimated \(\widehat{\lambda}_1,\dots,\widehat{\lambda}_{r}\) parameter. The ddsPLS methodology is based on a set of to be tested values for estimating each \(\lambda_s\), \(\forall s\in[\![1,r]\!]\), which is denoted as \(\Lambda\) in the following. Their values are pick in \([0,1]\).

We denote by \[ \bar{R}^{2}_{B}(\widehat{\mathcal{P}}_{\widehat{\lambda}_1,\dots,\widehat{\lambda}_{r}})~~~\text{and}~~~\bar{Q}^{2}_{B}(\widehat{\mathcal{P}}_{\widehat{\lambda}_1,\dots,\widehat{\lambda}_{r}}) \] the values of the two metrics for the estimated ddsPLS model \(\widehat{\mathcal{P}}\) of \(r\) components, based on the \(r\) estimated regularization coefficients \(\widehat{\lambda}_1,\dots,\widehat{\lambda}_{r}\). For each component, the ddsPLS methodology seeks the model which minimizes the difference between those two metrics, more precisely:

\[ \begin{array}{cccc} \widehat{\lambda}_r &=& \mbox{arg min}_{\lambda\in \Lambda} & \bar{R}^{2}_{B}( \widehat{\mathcal{P}}_{\widehat{\lambda}_1,\dots,\widehat{\lambda}_{r-1},\lambda} ) - \bar{Q}^{2}_{B}( \widehat{\mathcal{P}}_{\widehat{\lambda}_1,\dots,\widehat{\lambda}_{r-1},\lambda} ),\\ && \mbox{s.t} & \left\{\begin{array}{l} \bar{Q}^{2}_{B}( \widehat{\mathcal{P}}_{\widehat{\lambda}_1,\dots,\widehat{\lambda}_{r-1},\lambda} )>\bar{Q}^{2}_{B}( \widehat{\mathcal{P}}_{\widehat{\lambda}_1,\dots,\widehat{\lambda}_{r-1}} ),\\ \bar{Q}^{2}_{B,r}( \widehat{\mathcal{P}}_{\widehat{\lambda}_1,\dots,\widehat{\lambda}_{r-1},\lambda} )>0. \end{array}\right. \end{array} \]

Nota bene

The \(r^{th}\)-component is not built if \(\widehat{\lambda}_r =\emptyset\) and so the selected model is a \((r-1)\)-component (if \(r-1>0\) or is the mean estimator model) with estimated values for the regularization coefficients \((\widehat{\lambda}_1,\dots,\widehat{\lambda}_{r-1})\).

A high-dimensional structure

For the sack of the proposed simulations, we use \(R=2\) (associated to the dimension of \(\boldsymbol{\phi}\)), \(p=1000\) (associated to the dimension of \(\mathbf{A}'\boldsymbol{\phi}\) and \(\boldsymbol{\epsilon}\)) and \(q=3\) (associated to the dimension of \(\mathbf{D}'\boldsymbol{\phi}\) and \(\boldsymbol{\xi}\)).

The total number of such data-sets is rarely higher than \(p=1000\) and the constraint \(n<\!<p\) holds most of times. This is an high-dimensionnal data-set.

Each observation of those random vectors are generated following a multivariate normal distribution such as \[ \boldsymbol{\psi}=(\boldsymbol{\phi}',\boldsymbol{\epsilon}_{1,\dots,100}'/\sigma,\boldsymbol{\epsilon}_{101,\dots,1000}',\boldsymbol{\xi}_{1,2}',\xi_3)'\sim \mathcal{N}\left(\mathbf{0}_{2+1000+3},\mathbb{I}_{2+1000+3}\right), \]

where \(\sigma\) is the standard deviation of the additive noise. A response variable \(\mathbf{y}\) of \(q=3\) components is generated as a linear combination of the latent variable \(\boldsymbol{\phi}\) to which is added a Gaussian noise \(\boldsymbol{\xi}\). The equivalent process generates a variable \(\mathbf{x}\) of \(p=1000\) components, from the matrix \(\mathbf{A}\) and Gaussian additive noise \({\sigma}=\sqrt{1-0.95^2}\approx0.312\). The columns of \(\mathbf{A}\) and \(\mathbf{D}\) are normalized such as

\[ \forall(i,j)\in[\![1,p]\!]\times[\![1,q]\!],\ \mbox{var}({x}_i)=\mbox{var}({y}_j)=1. \]

Remark 2 Taking into account Remark 1 and the current statistical model, we can define a theoretical (understanding if \(n\rightarrow +\infty\)) value for \(\gamma(\mathcal{P})\) which is \[ \begin{array}{ccc} \gamma^\star & = & 1-\dfrac{ \sum_{j=1}^q \mbox{var} (\epsilon_j) }{ \sum_{j=1}^q \mbox{var} ({y}_{j})}= 2(1-\sigma^2)/3\approx0.602. \end{array} \] Comparing this theorethical value to the values of \(\bar{Q}^2_B\) helps interpretability, for theorethical work only. If \(\bar{Q}^2_B<\!<\gamma^\star\) then the corresponding model is not enough efficient. If \(\bar{Q}^2_B>\!>\gamma^\star\) then the corresponding model overfits the data.

Correlation structure

More precisely we propose to study the following matrices \(\mathbf{A}\) and \(\mathbf{D}\) \[ \begin{array}{c l c c} & \mathbf{A} =\sqrt{1-\sigma^2} \left(\begin{array}{ccc} \boldsymbol{1}_{50}' & \sqrt{\alpha}\boldsymbol{1}_{25}' & \boldsymbol{0}_{25}' & \boldsymbol{0}_{900}'\\ \boldsymbol{0}_{50}' & \sqrt{1-\alpha}\boldsymbol{1}_{25}' & \boldsymbol{0}_{25}' & \boldsymbol{0}_{900}'\\ \boldsymbol{0}_{50}' & \boldsymbol{0}_{25}' &\boldsymbol{1}_{25}' & \boldsymbol{0}_{900}'\\ \end{array} \right)_{(3,1000)} &\text{ and }& \mathbf{D}=\sqrt{1-\sigma^2} \left(\begin{array}{ccc} 1 & \sqrt{\beta_0} & 0 \\ 0 & \sqrt{1-\beta_0} & 0 \\ 0 & 0 & 0 \\ \end{array} \right)_{(3,3)}, \end{array} \] where \(\alpha\in [0,1]\) can be easily interpreted. Indeed, \(\alpha\) controls the correlation between the components \(\mathbf{x}_{1\dots 50}\) and \(\mathbf{x}_{51\dots 75}\). Also \(\beta_0=0.1\). It indirectly controls the association between \(\mathbf{x}\) and \(y_2\). The effects of \(\alpha\) on those two associations are detailed in the following table.

If \(\alpha> 1/2\), the second variable response component, \(y_2\), is hard to predict and the variables \(\mathbf{x}_{51\dots 75}\) are hard to select. The objective is to build models predicting \(y_1\) and \(y_2\) but not \(y_3\). Also components 1 to 75 of \(\mathbf{x}\) should be selected.

The function getData, printed in Appendix A.2 allows to simulate a couple (X,Y) according to this structure.

Compare with “non selection” ddsPLS model

As a popint of comparison, we can build the ddsPLS model for which \(\Lambda=\{0\}\) and look at its prediction performances, through the \(\bar{Q}_B^2\) statistics.

mo0 <- ddsPLS( datas$X, datas$Y,lambdas = 0,
               n_B=n_B,NCORES=NCORES,verbose = FALSE)
sum0 <- summary(mo0,return = TRUE)
print(sum0$R2Q2[,c(1,4)])

mo0 <- ddsPLS( datas$X, datas$Y,lambdas = 0,
               n_B=n_B,NCORES=1,verbose = FALSE)
sum0 <- summary(mo0,return = TRUE)
#>                       ______________
#>                      |    ddsPLS    |
#> =====================----------------=====================
#> The optimal ddsPLS model is built on 2 component(s)
#> 
#> The bootstrap quality statistics:
#> ---------------------------------
#>         lambda   R2 R2_r   Q2 Q2_r
#> Comp. 1      0 0.39 0.39 0.36 0.36
#> Comp. 2      0 0.62 0.23 0.55 0.29
#> 
#> 
#> The explained variance (in %):
#> -----------------------
#> 
#> In total: 61.75
#> -  -  -  
#> 
#> Per component or cumulated:
#> -  -  -  -  -  -  -  -  -  
#>                Comp. 1 Comp. 2
#> Per component    39.01   22.74
#> Cumulative       39.01   61.75
#> 
#> Per response variable:
#> -  -  -  -  -  -  -  -
#>            Y1    Y2   Y3
#> Comp. 1 83.36 33.33 0.35
#> Comp. 2  9.48 58.37 0.37
#> 
#> Per response variable per component:
#> -  -  -  -  -  -  -  -  -  -  -  -  
#>            Y1    Y2   Y3
#> Comp. 1 83.36 33.33 0.35
#> Comp. 2  9.48 58.37 0.37
#> 
#> ...and cumulated to:
#> -  -  -  -  -  -  - 
#>            Y1    Y2   Y3
#> Comp. 1 83.36 33.33 0.35
#> Comp. 2 92.84 91.70 0.72
#> 
#> For the X block:
#> -  -  -  -  -  -  -  -  -  
#>                Comp. 1 Comp. 2
#> Per component     5.22    2.41
#> Cumulative        5.22    7.63
#> =====================                =====================
#>                      ================

It is possible to compare the prediction qualities of the two models using

print(sum0$R2Q2[,c(1,4)])
#>         lambda        Q2
#> Comp. 1      0 0.3575840
#> Comp. 2      0 0.5457179
print(sum_up$R2Q2[,c(1,4)])
#>            lambda        Q2
#> Comp. 1 0.5172414 0.3246724
#> Comp. 2 0.4482759 0.5977299

In that context, the sparse ddsPLS approach allows to get better prediction rate than the “non selection” ddsPLS model.

The S3-method plot

It is also possible to plot different things thanks to the plot S3-method. In the representation with \(\lambda\) in abscissa:

• the \(\bullet\) points correspond to \(\lambda\) values for which the constraint of the optimization problem detailed in the section Automatic tunning of the parameters are active.
• the large \(\circ\) points corresponds to the selected value.

The different values given to the argument type would give representation that helps the analyst concluding on the final quality of the model. The different values are

• type="predict" to draw the predicted values of y against the observed. This can be useful to locate potential outliers (observations away from the distribution…) that would drive the model (… but close to the bisector).
• type="criterion" to draw the values of the metrics \(\bar{R}_{B}^2-\bar{Q}_{B}^2\). This is the optimized metrics.
• type="Q2" to draw the \(\bar{Q}_{B}^2\) metrics which represents the overall prediction quality of the build model, one component after another.
• type="prop" to draw the proportion of bootstrap models with a positive \(Q_{b,r}^2\), this \(\forall b \in [\![1,B]\!]\). Since this is tricky to interpret negative values for \(Q_{b,r}^2\) (apart from describing models which perform worse than the mean prediction model) negative values for \(\bar{Q}_{B,r}^2\) is necessarily hardly interpretable. However, the proportion of models with positive \(Q_{b,r}^2\) can be interpreted as a probability to finally build a model \(\mathcal{P}\) with a positive \(\gamma(\mathcal{P})\). This can be interesting to look at this metrics to interpret
• type="weightX" or type="weightY" to draw the values of the weights for each component for the X block of for the Y block. If there is an a priori, such as a functional one, on the variables of X (resp. Y), this a priori (which is not currently taken into account in the ddsPLS model) must certainly have an impact on the values of the weights parameters. This can be characterized by a structure of the weights on each component. In the opposite case, if the model is not enough sparse or too sparse, for example, the analyst is invited to modify the parameterization of the model, by limiting the grid of accessible \(\lambda\) for example.

Predicted versus observed values of the response

Simply specifying type="predict".

plot(mo,type="predict",legend.position = "topleft")

Since the the variable \(y_3\) has not been selected, its predicted values are constantly equal to the mean estimation. The two other columns of Y are described with more than 90% accuracy and no observation seems to guide the model at the expense of other observations.

The criterion \(\bar{R}_{B}^2-\bar{Q}_{B}^2\)

To plot the criterion, the plot argument type must be set to criterion. It is possible to move the legend with the legend.position argument.

plot(mo,type="criterion",legend.position = "top")

The legend title gives the total explained variance by the model built on the two components while the legend itself gives the explained variance by each of the considered component.

The \(\bar{Q}_{B}^2\) metrics

the previous figure does not provide information on the prediction quality of the model. This information can be found using the same S3-method fixing the parameter type to type="Q2".

plot(mo,type="Q2",legend.position = "bottomleft")

According to that figure it is clear that the chosen optimal model (minimizing \(\bar{R}_{B}^2-\bar{Q}_{B}^2\)) represents a model for which the \(\bar{Q}_{B}^2\) on each of its component is not far from being maximum.

The proportion of informational bootstrap models

plot(mo,type="prop",legend.position = "bottomleft")

plot(mo,type="weightY",mar=c(4,7,2,1))

plot(mo,type="weightX",cex.names = 0.5 )

A.1 Definition of the “component-\(r\) \(\bar{Q}_B^2\)” denoted \(\bar{Q}_{B,r}^2\)

For a given bootstrap sample indexed by \(b\), we define \[ \begin{array}{rccc} &Q^{2}_{b,r} & = & 1-\dfrac{ \left|\left|\mathbf{y}_{\text{OOB}(b)} -\hat{\mathbf{y}}_{\text{OOB}(b)}^{(r)}\right|\right|^2 }{ \left|\left|\mathbf{y}_{\text{OOB}(b)} -\bar{\mathbf{y}}_{\text{IN}(b)}^{(r-1)}\right|\right|^2 }, \end{array} \]

where, \(\hat{\mathbf{y}}_{\text{OOB}(b)}^{(r)}\) is the prediction of \({\mathbf{y}}_{\text{OOB}(b)}\) for the model built on \(r\) components. The interpretation of this metrics is that if \(Q^{2}_{b,r}>0\) then the \(r\)-components based model predicts better than the \(r-1\)-components based model. Naturally, the aggregated version of this metrics is

\[ \bar{Q}^{2}_{B,r}=\frac{1}{B}\sum_{b=1}^{B}Q^2_{b,r} \]

A.2 Function to simulate the high-dimensional data sets

We propose the following function

getData
#> function(n=100,alpha=0.4,beta_0=0.2,sigma=0.3,
#>                     p1=50,p2=25,p3=25,p=1000){
#>   R1 = R2 = R3 <- 1
#>   d <- R1+R2+R3
#>   A0 <- matrix(c(
#>     c(rep(1/sqrt(R1),p1),rep(sqrt(alpha),p2),rep(0,p3),rep(0,p-p1-p2-p3)),
#>     c(rep(0,p1),rep(sqrt(1-alpha),p2),rep(0,p3),rep(0,p-p1-p2-p3)),
#>     c(rep(0,p1),rep(0,p2),rep(1,p3),rep(0,p-p1-p2-p3))
#>   ),nrow = d,byrow = TRUE)
#>   A <- eps*A0
#>   D0 <- matrix(c(1,0,0,
#>                  sqrt(beta_0),sqrt(1-beta_0),0,
#>                  0,0,0),nrow = d,byrow = FALSE)
#>   D <- eps*D0
#>   q <- ncol(D)
#>   L_total <- q+p
#>   psi <- MASS::mvrnorm(n,mu = rep(0,d+L_total),Sigma = diag(d+L_total))
#>   phi <- psi[,1:d,drop=F]
#>   errorX <- matrix(rep(sqrt(1-apply(A^2,2,sum)),n),n,byrow = TRUE)
#>   errorY <- matrix(rep(sqrt(1-apply(D^2,2,sum)),n),n,byrow = TRUE)
#>   X <- phi%*%A + errorX*psi[,d+1:p,drop=F]
#>   Y <- phi%*%D + errorY*psi[,d+p+1:q,drop=F]
#>   list(X=X,Y=Y)
#> }
#> <bytecode: 0x110f55f58>

where the output is a list of two matrices X and Y.

A.3 Definitions of the explained variances…