Nothing
Data-Driven Sparse PLS (ddsPLS)
In ddsPLS: Data-Driven Sparse Partial Least Squares
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
library(ddsPLS)
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 = T)
A <- eps*A0
D0 <- matrix(c(1,0,0,
sqrt(beta_0),sqrt(1-beta_0),0,
0,0,0),nrow = d,byrow = F)
D <- eps*D0
q <- ncol(D)
L_total <- q+p
psi <- matrix(rnorm((d+L_total)*n,mean = 0,1),n)
phi <- psi[,1:d,drop=F]
errorX <- matrix(rep(sqrt(1-apply(A^2,2,sum)),n),n,byrow = T)
errorY <- matrix(rep(sqrt(1-apply(D^2,2,sum)),n),n,byrow = T)
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)
}
R2mean_diff_Q2mean <- cbind(
c(0.0390018049710666,0.0341751037010943,0.0300566345943933,0.0270372843280886,0.0251730242086019,0.0242227185400858,0.0234891342998435,0.0193133989423563,0.0169097166961483,0.0171263425383856,0.0170385680914704,0.0163279588077112,0.0152494718863702,0.0144263961042266,0.0139846253273588,0.0138793727811233,0.0138875727882883,0.0139142088056178,0.013947769332828,0.0151373262113938,0.0151964897734,0.0152786585690076,0.0184760931873834,0.0245523530157278,0.0329589020372058,0.0414793872444706,0.0487733189333913,0.0513953941897703,0.0513953941897703,0.0513953941897703)
,
c(0.0511147118390687,0.0366502939988904,0.0244192885247774,0.0157596067948765,0.0103693740213229,0.00774690072471196,0.00654474100829194,0.00602821715754076,0.00456584615102629,0.000647908840480049,-0.0012117730626553,-0.00324779272702536,-0.00339597202321107,-0.00349873504246012,-0.00347050350187528,-0.00343289518836143,-0.0033818961424037,-0.00330807936244493,-0.00318820602484982,-0.00296287672757123,-0.00264380624349336,-0.00161555612407471,-0.00511044138700312,-0.00517941187733331,-0.0027506954114751,-0.000122694727787143,0,0,0,0)
)
Q2mean <- cbind(
c(
0.377432184277843,0.380140639849587,0.381440573733071,0.381211078689433,0.379661476366436,0.377071669682569,0.373718997129306,0.371737203943991,0.366583285363937,0.357192628930804,0.346720293574788,0.337697081091153,0.332215307215219,0.330014817204357,0.329463268875761,0.329394099789224,0.329391348298594,0.32939602992879,0.329401599141234,0.32892978442128,0.328934064910596,0.328947551675465,0.327555460720922,0.325057907988492,0.321266212859659,0.318924669092771,0.315393736421526,0.295485665498064,0.295485665498064,0.295485665498064
),
c(
0.550769148647778,0.564051120615041,0.574049282463525,0.579753590218078,0.582232070714579,0.582828526524413,0.582822316455892,0.582729371555363,0.583602822903569,0.586270800202963,0.587558962483909,0.588999196398028,0.589131477829779,0.589208302519195,0.589207789772327,0.58920691960863,0.589205165067086,0.589201245713807,0.589190361787023,0.589150058389219,0.581193572554757,0.540787461471241,0.287353467160313,-0.174809002088594,-0.730811514935191,-0.992123542028864,-1,-1,-1,-1
)
)
PropQ2hPos <- cbind(
c(
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1
),
c(
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0.995,0.97,0.81,0.52,0.17,0.005,0,0,0,0
)
)
lambda_optim <- list(
c(TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE),
c(TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE),
c(FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE)
)
lambdaSparse <- c(0.5172414,0.4482759)
Overview
The objective is to explain a block-matrix $\mathbf{Y}$ thanks to a block-matrix $\mathbf{X}$. Each block describes $n$ observations through $q$ numerical variables for the block-matrix $\mathbf{Y}$ and $p$ numerical variables for the block-matrix $\mathbf{X}$. The links are assumed to be linear such as the objective is to estimated a linear matrix transformation $\mathbf{B}$ such as
$$
\mathbf{Y} = \mathbf{X}\mathbf{B} + \mathbf{E},
$$
where $\mathbf{E}$ is an additive random noise.
In the general case, the number of observations $n$ can be lower than the number of descriptors $p$ and most of regression methods cannot handle the estimation of the matrix $\mathbf{B}$, often denoted $n<!<p$. The ddsPLS methodology deals with this framework.
The objective of the proposed method is to estimate the matrix $\mathbf{B}$ and to simultaneously select the relevant variables in $\mathbf{Y}$ and in $\mathbf{X}$.
The methodology gives quality descriptors of the optimal built model, in terms of explained variance and prediction error, and provides a ranking of variable importance.
Tuning parameters
Only two kinds of parameters need to be tuned in the ddsPLS methodology. This tuning is automatically data driven. The proposed criterion uses the power of bootstrap, a classical statistical tool that generates different samples from an initial one, particularly interesting when the sample size $n$ is small.
The tuning parameters are:
- the number $R$ of independent links (latent components) between the matrix $\mathbf{Y}$ and the matrix $\mathbf{X}$.
- For each latent component $r\in [![1,R]!]$, the regularization coefficient $\lambda_r\in [0,1]$, which is easily interpretable as the minimum level of absolute correlation between the variables of $\mathbf{Y}$ and $\mathbf{X}$ allowed in the building of the current component. For instance, the $\mathbf{X}$ selected variables of a built component with $\lambda_r=0.2$ are at least correlated to 0.2 with the $\mathbf{Y}$ selected variables.
More precisely, the ddsPLS methodology detailed in this vignette is based on $R$ soft-thresholded estimations of the covariance matrices between a $q$-dimensional response matrix $\mathbf{Y}$ and a $p$-dimensional covariable matrix $\mathbf{X}$. Each soft-threshold parameter is the $\lambda_r$ parameter introduced above. It relies on the following latent variable model.
Quality descriptors
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})$.
Download and load the package
The package depends on a low number of low level packages. They are of two types:
-
Parallelization packages: foreach, parallel and doParallel (which depends on both of the previous packages).
-
$\mathbf{C}^{++}$-developer packages: Rcpp (for basic $\mathbf{C}^{++}$ development) and RcppEigen (for inner mathematical operations).
devtools::install_github("hlorenzo/ddsPLS")
library(ddsPLS)
library(ddsPLS)
A Latent Variable Model
In the following, we study a synthetic structure defined as, in general
$$
\left{
\begin{array}{l}
\mathbf{x} = \mathbf{A}'\boldsymbol{\phi} + \boldsymbol{\epsilon},\
\mathbf{A}\in\mathbb{R}^{R\times p},
\
\mathbf{y} = \mathbf{D}'\boldsymbol{\phi}+ \boldsymbol{\xi},\
\mathbf{D}\in\mathbb{R}^{R\times q},
\end{array}
\right.
$$
where $R$ is the total number of eigenvectors of $\mathbf{A}'\mathbf{A}=\mbox{var}\left(\mathbf{x}\right)$ with non-null projections on $\mathbf{A}'\mathbf{D}=\mbox{cov}\left(\mathbf{x},\mathbf{y}\right)$. In the PLS context this is the theoretical number of components.
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.
| Value of $\alpha$ | $\alpha\approx 0$ |$\alpha\approx 1$ |
|:--|:--:|:--:|
|$\mbox{cor}(\mathbf{x}{1\dots 50},\mathbf{x}{51\dots 75})$|Strong|Low|
|$\mbox{cor}(\mathbf{x},y_2)$|Strong|Low|
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.
A low correlated data-set
We assume that $\alpha=0.1$ and we simulate the data thanks to the function given in Appendix A.2. For this simulation, a sample of $n=200$ observations is chosen
eps <- 0.95
n <- 100
alpha <- 0.1
datas <- getData(n=n,alpha=alpha,beta_0=0.1,sigma=sqrt(1-eps^2),
p1=50,p2=25,p3=25,p=300)
It is also a necessity to fix the grid of $\lambda$, denoted as $\Lambda$ in the previous section and lambdas in the following. Also the number n_B of bootstrap samples is chosen, fixed to 20.
lambdas <- seq(0,1,length.out = 30)
n_B <- 20
Finally we have chosen to build the model working on NCORES=4 CPUs.
NCORES <- 4
mo <- ddsPLS( datas$X, datas$Y,lambdas = lambdas,
n_B=n_B,NCORES=NCORES,verbose = F)
mo <- ddsPLS( datas$X, datas$Y,lambdas = lambdaSparse,
n_B=n_B,NCORES=1,verbose = F)
Once the model is built, it is easy to check its performance thanks to the summary S3-method.
sum_up <- summary(mo,return = T)
A $2$-components model has been built and different metrics are detailed. We can discuss different points:
- the estimated values for $\lambda_1$ and $\lambda_2$ are 0.41 and 0.55 which means that no correlation between any X or Y variables lower than 0.41 will be considered for building the components and so for building the regression matrix.
- the variable $\mathbf{Y}_3$ is not selected.
- the metrics of prediction error is equal to $\bar{Q}_{B}^2=0.58$ which is below the theoretical value $\gamma^\star\approx0.602$.
Also, this S3-method can be used to get the list of the selected variables (and the figure of the Yes/No selction variable, as visible). The high dimension of the data does not allow to conclude directly, the following command allows to say that the variables selected in the X data set are exactly the one desired (since the two lists are exactly the same).
setdiff(1:75,mo$Selection$X)
Evaluate the quality of the model
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 = F)
sum0 <- summary(mo0,return = T)
print(sum0$R2Q2[,c(1,4)])
mo0 <- ddsPLS( datas$X, datas$Y,lambdas = 0,
n_B=n_B,NCORES=1,verbose = F)
sum0 <- summary(mo0,return = T)
It is possible to compare the prediction qualities of the two models using
print(sum0$R2Q2[,c(1,4)])
print(sum_up$R2Q2[,c(1,4)])
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 interprettype="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")
h_opt <- 2
matplot(lambdas,R2mean_diff_Q2mean,type = "l",ylab="",xlab=expression(lambda),
main=bquote(bar(R)[B]^2-bar(Q)[B]^2))
for(s in 1:h_opt){
points(lambdas,R2mean_diff_Q2mean[,s],type = "p",pch=16,cex=lambda_optim[[s]],col=s)
points(lambdaSparse[s],R2mean_diff_Q2mean[which.min(abs(lambdas-lambdaSparse[s])),s],pch=1,cex=2,col=s)
}
legend("top",paste("Comp.",1:h_opt," (",round(mo$varExplained$Comp),"%)",sep=""),
col = 1:h_opt,pch=16,bty = "n",
title = paste("Total explained variance ",round(mo$varExplained$Cumu)[h_opt],"%",sep=""))
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")
h_opt <- 2
matplot(lambdas,Q2mean,type = "l",ylab="",xlab=expression(lambda),
main=bquote(bar(R)[B]^2-bar(Q)[B]^2))
for(s in 1:h_opt){
points(lambdas,Q2mean[,s],type = "p",pch=16,cex=lambda_optim[[s]],col=s)
points(lambdaSparse[s],Q2mean[which.min(abs(lambdas-lambdaSparse[s])),s],pch=1,cex=2,col=s)
}
abline(h=0,lwd=2,lty=2)
legend("bottomleft",paste("Comp.",1:h_opt," (",round(mo$varExplained$Comp),"%)",sep=""),
col = 1:h_opt,pch=16,bty = "n",
title = paste("Total explained variance ",round(mo$varExplained$Cumu)[h_opt],"%",sep=""))
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")
h_opt <- 2
# Plot of Prop of positive Q2h
matplot(lambdas,PropQ2hPos,type = "l",ylab="",xlab=expression(lambda),
main=bquote("Proportion of models with positive"~Q["b,r"]^2))
abline(h=((1:10)/10)[-5],col="gray80",lwd=0.5,lty=3)
abline(h=5/10,col="gray60",lwd=0.7,lty=1)
text(min(lambdas),1/2,labels = "1/2",pos = 4,col="gray40")
for(s in 1:h_opt){
points(lambdas,PropQ2hPos[,s],type = "p",pch=16,cex=mo$lambda_optim[[s]],col=s)
points(lambdaSparse[s],PropQ2hPos[which.min(abs(lambdas-lambdaSparse[s])),s],pch=1,cex=2,col=s)
}
legend("bottomleft",paste("Comp.",1:h_opt," (",round(mo$varExplained$Comp),"%)",sep=""),
col = 1:h_opt,pch=16,bty = "n",
title = paste("Total explained variance ",round(mo$varExplained$Cumu)[h_opt],"%",sep=""))
Plot of the weights
It is also possible to directly plot the values of the weights for each component. For the Y block and for the X block, such as
plot(mo,type="weightY",mar=c(4,7,2,1))
where it is clear that the first component explains $y_1$ while the second one explains $y_2$. Looking at the weights on the block X:
plot(mo,type="weightX",cex.names = 0.5 )
the variables $\mathbf{x}{1\dots 50}$ (resp. $\mathbf{x}{51\dots 75}$) are selected on the first (resp. second) component, which is associated only with $y_1$ (resp. $y_2$).
Appendix
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
where the output is a list of two matrices X and Y.
A.3 Definitions of the explained variances...
A.3.1 Cumulated
For a model built on $R$ components
$$
\begin{array}{rccc}
&\text{ExpVar}{1:R} & = &
\left(1-\frac{1}{n}\sum{i=1}^n\dfrac{
\left|\left|\mathbf{y}i
-\hat{\mathbf{y}}^{(1:R)}_i\right|\right|^2
}{
\left|\left|\mathbf{y}_i
-\boldsymbol{\mu}{\mathbf{y}}\right|\right|^2
}
\right)*100,
\end{array}
$$
where
$$
\hat{\mathbf{y}}^{(1:R)}
=
\left(
\mathbf{X}-\boldsymbol{\mu}{\mathbf{x}}
\right)
\mathbf{U}{(1:R)}\left(\mathbf{P}{(1:R)}'\mathbf{U}{(1:R)}\right)^{-1}\mathbf{C}{(1:R)}
+
\boldsymbol{\mu}{\mathbf{y}},
$$
and $\mathbf{U}{(1:R)}$ is the concatenation of the $R$ weights, $\mathbf{P}{(1:R)}$ is the concatenation of the $R$ scores $\mathbf{p}r=\mathbf{X}^{(r)'}\mathbf{t}_r/\mathbf{t}_r'\mathbf{t}_r$ and $\mathbf{C}{(1:R)}$ is the concatenation of the $R$ scores $\mathbf{c}r=\boldsymbol{\Pi}_r\mathbf{Y}^{(r)'}\mathbf{t}_r/\mathbf{t}_r'\mathbf{t}_r$. The deflated matrices are defined such as
$$
\mathbf{X}{^{(r+1)}} =\mathbf{X}^{(r)}-\mathbf{t}_r\mathbf{p}_r',\ \mathbf{Y}{^{(r+1)}} = \mathbf{Y}{^{(r)}} - \mathbf{t}_r\mathbf{c}_r',
$$
and $\boldsymbol{\Pi}_r$ is the diagonal matrix with 1 element if the associated response variable is selected and 0 elsewhere. Also, $\boldsymbol{\mu}{\mathbf{y}}$ and $\boldsymbol{\mu}_{\mathbf{x}}$ are the estimated mean matrices.
A.3.2 Per component
$$
\begin{array}{rccc}
&\text{ExpVar}{R} & = &
\left(1-\frac{1}{n}\sum{i=1}^n\dfrac{
\left|\left|\mathbf{y}i
-\hat{\mathbf{y}}^{(R)}_i\right|\right|^2
}{
\left|\left|\mathbf{y}_i
-\boldsymbol{\mu}{\mathbf{y}}\right|\right|^2
}
\right)*100,
\end{array}
$$
A.3.3 Per response variable
$$
\begin{array}{rccc}
&\text{ExpVar}{1:R}^{(j)} & = &
\left( 1-\frac{1}{n}\sum{i=1}^n\dfrac{
\left|\left|y_{i,j}
-\hat{{y}}^{(1:R)}{i,j}\right|\right|^2
}{
\left|\left|{y}{i,j}
-{\mu}_{{y}_j}\right|\right|^2
}
\right)*100,
\end{array}
$$
A.3.4 Per response variable per component
$$
\begin{array}{rccc}
&\text{ExpVar}{R}^{(j)} & = &
\left(
1-\frac{1}{n}\sum{i=1}^n\dfrac{
\left|\left|y_{i,j}
-\hat{{y}}^{(R)}{i,j}\right|\right|^2
}{
\left|\left|{y}{i,j}
-{\mu}_{{y}_j}\right|\right|^2
}
\right)*100,
\end{array}
$$
Try the ddsPLS package in your browser
Any scripts or data that you put into this service are public.
ddsPLS documentation built on May 31, 2023, 7:50 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
library(ddsPLS)
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 = T) A <- eps*A0 D0 <- matrix(c(1,0,0, sqrt(beta_0),sqrt(1-beta_0),0, 0,0,0),nrow = d,byrow = F) D <- eps*D0 q <- ncol(D) L_total <- q+p psi <- matrix(rnorm((d+L_total)*n,mean = 0,1),n) phi <- psi[,1:d,drop=F] errorX <- matrix(rep(sqrt(1-apply(A^2,2,sum)),n),n,byrow = T) errorY <- matrix(rep(sqrt(1-apply(D^2,2,sum)),n),n,byrow = T) 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) } R2mean_diff_Q2mean <- cbind( c(0.0390018049710666,0.0341751037010943,0.0300566345943933,0.0270372843280886,0.0251730242086019,0.0242227185400858,0.0234891342998435,0.0193133989423563,0.0169097166961483,0.0171263425383856,0.0170385680914704,0.0163279588077112,0.0152494718863702,0.0144263961042266,0.0139846253273588,0.0138793727811233,0.0138875727882883,0.0139142088056178,0.013947769332828,0.0151373262113938,0.0151964897734,0.0152786585690076,0.0184760931873834,0.0245523530157278,0.0329589020372058,0.0414793872444706,0.0487733189333913,0.0513953941897703,0.0513953941897703,0.0513953941897703) , c(0.0511147118390687,0.0366502939988904,0.0244192885247774,0.0157596067948765,0.0103693740213229,0.00774690072471196,0.00654474100829194,0.00602821715754076,0.00456584615102629,0.000647908840480049,-0.0012117730626553,-0.00324779272702536,-0.00339597202321107,-0.00349873504246012,-0.00347050350187528,-0.00343289518836143,-0.0033818961424037,-0.00330807936244493,-0.00318820602484982,-0.00296287672757123,-0.00264380624349336,-0.00161555612407471,-0.00511044138700312,-0.00517941187733331,-0.0027506954114751,-0.000122694727787143,0,0,0,0) ) Q2mean <- cbind( c( 0.377432184277843,0.380140639849587,0.381440573733071,0.381211078689433,0.379661476366436,0.377071669682569,0.373718997129306,0.371737203943991,0.366583285363937,0.357192628930804,0.346720293574788,0.337697081091153,0.332215307215219,0.330014817204357,0.329463268875761,0.329394099789224,0.329391348298594,0.32939602992879,0.329401599141234,0.32892978442128,0.328934064910596,0.328947551675465,0.327555460720922,0.325057907988492,0.321266212859659,0.318924669092771,0.315393736421526,0.295485665498064,0.295485665498064,0.295485665498064 ), c( 0.550769148647778,0.564051120615041,0.574049282463525,0.579753590218078,0.582232070714579,0.582828526524413,0.582822316455892,0.582729371555363,0.583602822903569,0.586270800202963,0.587558962483909,0.588999196398028,0.589131477829779,0.589208302519195,0.589207789772327,0.58920691960863,0.589205165067086,0.589201245713807,0.589190361787023,0.589150058389219,0.581193572554757,0.540787461471241,0.287353467160313,-0.174809002088594,-0.730811514935191,-0.992123542028864,-1,-1,-1,-1 ) ) PropQ2hPos <- cbind( c( 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1 ), c( 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0.995,0.97,0.81,0.52,0.17,0.005,0,0,0,0 ) ) lambda_optim <- list( c(TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE), c(TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE), c(FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE) ) lambdaSparse <- c(0.5172414,0.4482759)
Overview
The objective is to explain a block-matrix $\mathbf{Y}$ thanks to a block-matrix $\mathbf{X}$. Each block describes $n$ observations through $q$ numerical variables for the block-matrix $\mathbf{Y}$ and $p$ numerical variables for the block-matrix $\mathbf{X}$. The links are assumed to be linear such as the objective is to estimated a linear matrix transformation $\mathbf{B}$ such as
$$ \mathbf{Y} = \mathbf{X}\mathbf{B} + \mathbf{E}, $$ where $\mathbf{E}$ is an additive random noise.
In the general case, the number of observations $n$ can be lower than the number of descriptors $p$ and most of regression methods cannot handle the estimation of the matrix $\mathbf{B}$, often denoted $n<!<p$. The ddsPLS methodology deals with this framework.
The objective of the proposed method is to estimate the matrix $\mathbf{B}$ and to simultaneously select the relevant variables in $\mathbf{Y}$ and in $\mathbf{X}$.
The methodology gives quality descriptors of the optimal built model, in terms of explained variance and prediction error, and provides a ranking of variable importance.
Tuning parameters
Only two kinds of parameters need to be tuned in the ddsPLS methodology. This tuning is automatically data driven. The proposed criterion uses the power of bootstrap, a classical statistical tool that generates different samples from an initial one, particularly interesting when the sample size $n$ is small.
The tuning parameters are:
- the number $R$ of independent links (latent components) between the matrix $\mathbf{Y}$ and the matrix $\mathbf{X}$.
- For each latent component $r\in [![1,R]!]$, the regularization coefficient $\lambda_r\in [0,1]$, which is easily interpretable as the minimum level of absolute correlation between the variables of $\mathbf{Y}$ and $\mathbf{X}$ allowed in the building of the current component. For instance, the $\mathbf{X}$ selected variables of a built component with $\lambda_r=0.2$ are at least correlated to 0.2 with the $\mathbf{Y}$ selected variables.
More precisely, the ddsPLS methodology detailed in this vignette is based on $R$ soft-thresholded estimations of the covariance matrices between a $q$-dimensional response matrix $\mathbf{Y}$ and a $p$-dimensional covariable matrix $\mathbf{X}$. Each soft-threshold parameter is the $\lambda_r$ parameter introduced above. It relies on the following latent variable model.
Quality descriptors
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})$.
Download and load the package
The package depends on a low number of low level packages. They are of two types:
-
Parallelization packages:
foreach,parallelanddoParallel(which depends on both of the previous packages). -
$\mathbf{C}^{++}$-developer packages:
Rcpp(for basic $\mathbf{C}^{++}$ development) andRcppEigen(for inner mathematical operations).
devtools::install_github("hlorenzo/ddsPLS") library(ddsPLS)
library(ddsPLS)
A Latent Variable Model
In the following, we study a synthetic structure defined as, in general
$$
\left{
\begin{array}{l}
\mathbf{x} = \mathbf{A}'\boldsymbol{\phi} + \boldsymbol{\epsilon},\
\mathbf{A}\in\mathbb{R}^{R\times p},
\
\mathbf{y} = \mathbf{D}'\boldsymbol{\phi}+ \boldsymbol{\xi},\
\mathbf{D}\in\mathbb{R}^{R\times q},
\end{array}
\right.
$$
where $R$ is the total number of eigenvectors of $\mathbf{A}'\mathbf{A}=\mbox{var}\left(\mathbf{x}\right)$ with non-null projections on $\mathbf{A}'\mathbf{D}=\mbox{cov}\left(\mathbf{x},\mathbf{y}\right)$. In the PLS context this is the theoretical number of components.
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.
| Value of $\alpha$ | $\alpha\approx 0$ |$\alpha\approx 1$ | |:--|:--:|:--:| |$\mbox{cor}(\mathbf{x}{1\dots 50},\mathbf{x}{51\dots 75})$|Strong|Low| |$\mbox{cor}(\mathbf{x},y_2)$|Strong|Low|
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.
A low correlated data-set
We assume that $\alpha=0.1$ and we simulate the data thanks to the function given in Appendix A.2. For this simulation, a sample of $n=200$ observations is chosen
eps <- 0.95 n <- 100 alpha <- 0.1 datas <- getData(n=n,alpha=alpha,beta_0=0.1,sigma=sqrt(1-eps^2), p1=50,p2=25,p3=25,p=300)
It is also a necessity to fix the grid of $\lambda$, denoted as $\Lambda$ in the previous section and lambdas in the following. Also the number n_B of bootstrap samples is chosen, fixed to 20.
lambdas <- seq(0,1,length.out = 30) n_B <- 20
Finally we have chosen to build the model working on NCORES=4 CPUs.
NCORES <- 4 mo <- ddsPLS( datas$X, datas$Y,lambdas = lambdas, n_B=n_B,NCORES=NCORES,verbose = F)
mo <- ddsPLS( datas$X, datas$Y,lambdas = lambdaSparse, n_B=n_B,NCORES=1,verbose = F)
Once the model is built, it is easy to check its performance thanks to the summary S3-method.
sum_up <- summary(mo,return = T)
A $2$-components model has been built and different metrics are detailed. We can discuss different points:
- the estimated values for $\lambda_1$ and $\lambda_2$ are 0.41 and 0.55 which means that no correlation between any X or Y variables lower than 0.41 will be considered for building the components and so for building the regression matrix.
- the variable $\mathbf{Y}_3$ is not selected.
- the metrics of prediction error is equal to $\bar{Q}_{B}^2=0.58$ which is below the theoretical value $\gamma^\star\approx0.602$.
Also, this S3-method can be used to get the list of the selected variables (and the figure of the Yes/No selction variable, as visible). The high dimension of the data does not allow to conclude directly, the following command allows to say that the variables selected in the X data set are exactly the one desired (since the two lists are exactly the same).
setdiff(1:75,mo$Selection$X)
Evaluate the quality of the model
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 = F) sum0 <- summary(mo0,return = T) print(sum0$R2Q2[,c(1,4)])
mo0 <- ddsPLS( datas$X, datas$Y,lambdas = 0, n_B=n_B,NCORES=1,verbose = F) sum0 <- summary(mo0,return = T)
It is possible to compare the prediction qualities of the two models using
print(sum0$R2Q2[,c(1,4)]) print(sum_up$R2Q2[,c(1,4)])
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 interprettype="weightX"ortype="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")
h_opt <- 2 matplot(lambdas,R2mean_diff_Q2mean,type = "l",ylab="",xlab=expression(lambda), main=bquote(bar(R)[B]^2-bar(Q)[B]^2)) for(s in 1:h_opt){ points(lambdas,R2mean_diff_Q2mean[,s],type = "p",pch=16,cex=lambda_optim[[s]],col=s) points(lambdaSparse[s],R2mean_diff_Q2mean[which.min(abs(lambdas-lambdaSparse[s])),s],pch=1,cex=2,col=s) } legend("top",paste("Comp.",1:h_opt," (",round(mo$varExplained$Comp),"%)",sep=""), col = 1:h_opt,pch=16,bty = "n", title = paste("Total explained variance ",round(mo$varExplained$Cumu)[h_opt],"%",sep=""))
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")
h_opt <- 2 matplot(lambdas,Q2mean,type = "l",ylab="",xlab=expression(lambda), main=bquote(bar(R)[B]^2-bar(Q)[B]^2)) for(s in 1:h_opt){ points(lambdas,Q2mean[,s],type = "p",pch=16,cex=lambda_optim[[s]],col=s) points(lambdaSparse[s],Q2mean[which.min(abs(lambdas-lambdaSparse[s])),s],pch=1,cex=2,col=s) } abline(h=0,lwd=2,lty=2) legend("bottomleft",paste("Comp.",1:h_opt," (",round(mo$varExplained$Comp),"%)",sep=""), col = 1:h_opt,pch=16,bty = "n", title = paste("Total explained variance ",round(mo$varExplained$Cumu)[h_opt],"%",sep=""))
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")
h_opt <- 2 # Plot of Prop of positive Q2h matplot(lambdas,PropQ2hPos,type = "l",ylab="",xlab=expression(lambda), main=bquote("Proportion of models with positive"~Q["b,r"]^2)) abline(h=((1:10)/10)[-5],col="gray80",lwd=0.5,lty=3) abline(h=5/10,col="gray60",lwd=0.7,lty=1) text(min(lambdas),1/2,labels = "1/2",pos = 4,col="gray40") for(s in 1:h_opt){ points(lambdas,PropQ2hPos[,s],type = "p",pch=16,cex=mo$lambda_optim[[s]],col=s) points(lambdaSparse[s],PropQ2hPos[which.min(abs(lambdas-lambdaSparse[s])),s],pch=1,cex=2,col=s) } legend("bottomleft",paste("Comp.",1:h_opt," (",round(mo$varExplained$Comp),"%)",sep=""), col = 1:h_opt,pch=16,bty = "n", title = paste("Total explained variance ",round(mo$varExplained$Cumu)[h_opt],"%",sep=""))
Plot of the weights
It is also possible to directly plot the values of the weights for each component. For the Y block and for the X block, such as
plot(mo,type="weightY",mar=c(4,7,2,1))
where it is clear that the first component explains $y_1$ while the second one explains $y_2$. Looking at the weights on the block X:
plot(mo,type="weightX",cex.names = 0.5 )
the variables $\mathbf{x}{1\dots 50}$ (resp. $\mathbf{x}{51\dots 75}$) are selected on the first (resp. second) component, which is associated only with $y_1$ (resp. $y_2$).
Appendix
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
where the output is a list of two matrices X and Y.
A.3 Definitions of the explained variances...
A.3.1 Cumulated
For a model built on $R$ components
$$ \begin{array}{rccc} &\text{ExpVar}{1:R} & = & \left(1-\frac{1}{n}\sum{i=1}^n\dfrac{ \left|\left|\mathbf{y}i -\hat{\mathbf{y}}^{(1:R)}_i\right|\right|^2 }{ \left|\left|\mathbf{y}_i -\boldsymbol{\mu}{\mathbf{y}}\right|\right|^2 } \right)*100, \end{array} $$
where
$$ \hat{\mathbf{y}}^{(1:R)} = \left( \mathbf{X}-\boldsymbol{\mu}{\mathbf{x}} \right) \mathbf{U}{(1:R)}\left(\mathbf{P}{(1:R)}'\mathbf{U}{(1:R)}\right)^{-1}\mathbf{C}{(1:R)} + \boldsymbol{\mu}{\mathbf{y}}, $$ and $\mathbf{U}{(1:R)}$ is the concatenation of the $R$ weights, $\mathbf{P}{(1:R)}$ is the concatenation of the $R$ scores $\mathbf{p}r=\mathbf{X}^{(r)'}\mathbf{t}_r/\mathbf{t}_r'\mathbf{t}_r$ and $\mathbf{C}{(1:R)}$ is the concatenation of the $R$ scores $\mathbf{c}r=\boldsymbol{\Pi}_r\mathbf{Y}^{(r)'}\mathbf{t}_r/\mathbf{t}_r'\mathbf{t}_r$. The deflated matrices are defined such as $$ \mathbf{X}{^{(r+1)}} =\mathbf{X}^{(r)}-\mathbf{t}_r\mathbf{p}_r',\ \mathbf{Y}{^{(r+1)}} = \mathbf{Y}{^{(r)}} - \mathbf{t}_r\mathbf{c}_r', $$ and $\boldsymbol{\Pi}_r$ is the diagonal matrix with 1 element if the associated response variable is selected and 0 elsewhere. Also, $\boldsymbol{\mu}{\mathbf{y}}$ and $\boldsymbol{\mu}_{\mathbf{x}}$ are the estimated mean matrices.
A.3.2 Per component
$$ \begin{array}{rccc} &\text{ExpVar}{R} & = & \left(1-\frac{1}{n}\sum{i=1}^n\dfrac{ \left|\left|\mathbf{y}i -\hat{\mathbf{y}}^{(R)}_i\right|\right|^2 }{ \left|\left|\mathbf{y}_i -\boldsymbol{\mu}{\mathbf{y}}\right|\right|^2 } \right)*100, \end{array} $$
A.3.3 Per response variable
$$ \begin{array}{rccc} &\text{ExpVar}{1:R}^{(j)} & = & \left( 1-\frac{1}{n}\sum{i=1}^n\dfrac{ \left|\left|y_{i,j} -\hat{{y}}^{(1:R)}{i,j}\right|\right|^2 }{ \left|\left|{y}{i,j} -{\mu}_{{y}_j}\right|\right|^2 } \right)*100, \end{array} $$
A.3.4 Per response variable per component
$$ \begin{array}{rccc} &\text{ExpVar}{R}^{(j)} & = & \left( 1-\frac{1}{n}\sum{i=1}^n\dfrac{ \left|\left|y_{i,j} -\hat{{y}}^{(R)}{i,j}\right|\right|^2 }{ \left|\left|{y}{i,j} -{\mu}_{{y}_j}\right|\right|^2 } \right)*100, \end{array} $$
Try the ddsPLS package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.