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², evaluates the precision of the model on the training data-set. The second one, denoted as Q², evaluates the precision of the model on a test data-set, independent from the train data-set. The R² 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² − Q² is minimum.

The ddsPLS methodology is based on bootstrap versions of the R² and the Q² defined below.

Bootstrap versions of the R² and the Q²

More precisely, for the bootstrap sample of index b, among a total of B bootstrap samples, the R_b² and the Q_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 “y” by the current model (based on the in-bag sample) and “$\bar{\mathbf{y}}$” stands for the mean estimator of “y”. The ddsPLS methodology aggregates the B descriptors R_b² and Q_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 R̄_B² and Q̄_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 Q̄_B, r² is the “component-version” of Q̄_B²” defined in Appendix A.1.

Remark 1 The metrics R̄_B² and Q̄_B² are in fact estimators of a same statistic γ, associated with a prediction model 𝒫: $$ \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^(𝒫)) is the prediction of the j^th component of y by the model 𝒫. The closer this statistic is to 1, the more accurate the model 𝒫 is.

Automatic tunning of the parameters

As said before, a ddsPLS model of r components is based on r parameters (λ₁, …, λ_r), denoted 𝒫_{λ1, …, λr}. The metrics R̄_B², Q̄_B² or Q̄_B, r² are functions of the values taken by the estimated ddsPLS models $\widehat{\mathcal{P}}$, which depend on the r estimated λ̂₁, …, λ̂_r parameter. The ddsPLS methodology is based on a set of to be tested values for estimating each λ_s, ∀s ∈ [ [1, r] ], which is denoted as Λ 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 λ̂₁, …, λ̂_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 λ̂_r = ∅ 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 (λ̂₁, …, λ̂_r − 1).

A high-dimensional structure

For the sack of the proposed simulations, we use R = 2 (associated to the dimension of ϕ), p = 1000 (associated to the dimension of A′ϕ and ϵ) and q = 3 (associated to the dimension of D′ϕ and ξ).

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 ψ = (ϕ′, ϵ_{1, …, 100}′/σ, ϵ_{101, …, 1000}′, ξ_1, 2′, ξ₃)′ ∼ 𝒩(0_{2 + 1000 + 3}, 𝕀_{2 + 1000 + 3}),

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

∀(i, j) ∈ [ [1, p] ] × [ [1, q] ], var(x_i) = var(y_j) = 1.

Remark 2 Taking into account Remark 1 and the current statistical model, we can define a theoretical (understanding if n → +∞) value for γ(𝒫) 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 Q̄_B² helps interpretability, for theorethical work only. If Q̄_B²<  < γ^⋆ then the corresponding model is not enough efficient. If Q̄_B²>  > γ^⋆ then the corresponding model overfits the data.

Correlation structure

More precisely we propose to study the following matrices A and 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 α ∈ [0, 1] can be easily interpreted. Indeed, α controls the correlation between the components x_1…50 and x_51…75. Also β₀ = 0.1. It indirectly controls the association between x and y₂. The effects of α on those two associations are detailed in the following table.

If α > 1/2, the second variable response component, y₂, is hard to predict and the variables x_51…75 are hard to select. The objective is to build models predicting y₁ and y₂ but not y₃. Also components 1 to 75 of 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 Λ = {0} and look at its prediction performances, through the Q̄_B² 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.63 0.24 0.56 0.31
#> 
#> 
#> The explained variance (in %):
#> -----------------------
#> 
#> In total: 62.83
#> -  -  -  
#> 
#> Per component or cumulated:
#> -  -  -  -  -  -  -  -  -  
#>                Comp. 1 Comp. 2
#> Per component     38.9   23.93
#> Cumulative        38.9   62.83
#> 
#> Per response variable:
#> -  -  -  -  -  -  -  -
#>            Y1    Y2   Y3
#> Comp. 1 78.74 34.30 3.65
#> Comp. 2 13.90 57.78 0.11
#> 
#> Per response variable per component:
#> -  -  -  -  -  -  -  -  -  -  -  -  
#>            Y1    Y2   Y3
#> Comp. 1 78.74 34.30 3.65
#> Comp. 2 13.90 57.78 0.11
#> 
#> ...and cumulated to:
#> -  -  -  -  -  -  - 
#>            Y1    Y2   Y3
#> Comp. 1 78.74 34.30 3.65
#> Comp. 2 92.64 92.09 3.75
#> 
#> For the X block:
#> -  -  -  -  -  -  -  -  -  
#>                Comp. 1 Comp. 2
#> Per component     5.32    2.45
#> Cumulative        5.32    7.76
#> =====================                =====================
#>                      ================

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.3575314
#> Comp. 2      0 0.5582657
print(sum_up$R2Q2[,c(1,4)])
#>            lambda        Q2
#> Comp. 1 0.5172414 0.3314618
#> Comp. 2 0.5172414 0.5982472

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 λ in abscissa:

• the • points correspond to λ values for which the constraint of the optimization problem detailed in the section Automatic tunning of the parameters are active.
• the large ∘ 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 R̄_B² − Q̄_B². This is the optimized metrics.
• type="Q2" to draw the Q̄_B² 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², this ∀b ∈ [ [1, B] ]. Since this is tricky to interpret negative values for Q_b, r² (apart from describing models which perform worse than the mean prediction model) negative values for Q̄_B, r² is necessarily hardly interpretable. However, the proportion of models with positive Q_b, r² can be interpreted as a probability to finally build a model 𝒫 with a positive γ(𝒫). 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 λ 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₃ 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 R̄_B² − Q̄_B²

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 Q̄_B² 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 R̄_B² − Q̄_B²) represents a model for which the Q̄_B² 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 Q̄_B²” denoted Q̄_B, r²

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 y_OOB(b) for the model built on r components. The interpretation of this metrics is that if Q_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: 0x5640407b3a40>

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

A.3 Definitions of the explained variances…