Example for msma

2024-02-14 @Atsushi Kawaguchi

The msma package provides functions for a matrix decomposition method incorporating sparse and supervised modeling for a multiblock multivariable data analysis.

Preparation

Install package (as necessary)

if(!require("msma")) install.packages("msma")

Load package

library(msma)

Getting started

Simulated multiblock data (list data) by using the function simdata.

Sample size is 50. The correlation coeficient is 0.8. The numbers of columns for response and predictor can be specified by the argument Yps and Xps, respectively. The length of vecor represents the number of blocks. That is, response has three blocks with the numbers of columns being 3, 4, and 5 and predictor has one block with the number of columns being 3.

dataset0 = simdata(n = 50, rho = 0.8, Yps = c(3, 4, 5), Xps = 3, seed=1)
X0 = dataset0$X; Y0 = dataset0$Y 

The data generated here is applied to the msma function.

One Component

The argument comp can specify the number of components. The arguments lambdaX and lambdaY can specify the regularization parameters for X and Y, respectively.

First, we set comp=1, which will perform an analysis with 1 component.

fit01 = msma(X0, Y0, comp=1, lambdaX=0.05, lambdaY=1:3)
fit01

## Call:
## msma.default(X = X0, Y = Y0, comp = 1, lambdaX = 0.05, lambdaY = 1:3)
## 
## Numbers of non-zeros for X block: 
##        comp1
## block1     3
## 
## Numbers of non-zeros for X super: 
##         comp1
## comp1-1     1
## 
## Numbers of non-zeros for Y block: 
##        comp1
## block1     1
## block2     1
## block3     1
## 
## Numbers of non-zeros for Y super: 
##         comp1
## comp1-1     3

The plot function is available. In default setting, the block weights are displayed as a barplot.

plot(fit01)

Two Component

Next, we set comp=2, which will perform an analysis with 2 components.

fit02 = msma(X0, Y0, comp=2, lambdaX=0.03, lambdaY=0.01*(1:3))
fit02

## Call:
## msma.default(X = X0, Y = Y0, comp = 2, lambdaX = 0.03, lambdaY = 0.01 * 
##     (1:3))
## 
## Numbers of non-zeros for X block: 
##        comp1 comp2
## block1     3     3
## 
## Numbers of non-zeros for X super: 
##         comp1 comp2
## comp1-1     1     1
## 
## Numbers of non-zeros for Y block: 
##        comp1 comp2
## block1     3     3
## block2     4     4
## block3     5     5
## 
## Numbers of non-zeros for Y super: 
##         comp1 comp2
## comp1-1     3     3

Single Block

Two matrics are prepared by specifying arguments Yps and Xps.

dataset1 = simdata(n = 50, rho = 0.8, Yps = 5, Xps = 5, seed=1)
X1 = dataset1$X[[1]]; Y1 = dataset1$Y 

Principal Component Analysis (PCA)

If input is a matrix, a principal component analysis is implemented.

(fit111 = msma(X1, comp=5))

## Call:
## msma.default(X = X1, comp = 5)
## 
## Numbers of non-zeros for X block: 
##        comp1 comp2 comp3 comp4 comp5
## block1     5     5     5     5     5
## 
## Numbers of non-zeros for X super: 
##         comp1 comp2 comp3 comp4 comp5
## comp1-1     1     1     1     1     1

The weight (loading) vectors can be obtained as follows.

fit111$wbX

## $block1
##           comp1       comp2       comp3         comp4       comp5
## X.1.1 0.4309622 -0.74170223 -0.03672379  0.1325580413 -0.49520613
## X.1.2 0.4483196  0.31188303  0.63228246  0.5490205405  0.02310504
## X.1.3 0.4601597 -0.19547078 -0.38567734  0.1474129336  0.76129277
## X.1.4 0.4392794  0.55811865 -0.57117598 -0.0006449093 -0.41145448
## X.1.5 0.4566923  0.05386584  0.35196769 -0.8119567864  0.07331836

The bar plots of weight vectors are provided by the function plot. The component number is specified by the argument axes. The plot type is selected by the argument plottype. Furthermore, since this function uses the barplot function originally built into R, its arguments are also available. In the following example, on the horizontal axis, the magnification of the variable names is set to 0.7 by setting cex.names=0.7, and the variable names are oriented as las=2.

par(mfrow=c(1,2))
plot(fit111, axes = 1, plottype="bar", cex.names=0.7, las=2)
plot(fit111, axes = 2, plottype="bar", cex.names=0.7, las=2)

The score vectors for first six subjects.

lapply(fit111$sbX, head)

## $block1
##            [,1]          [,2]        [,3]        [,4]        [,5]
## [1,]  0.7097369  0.0487564120  0.10746733 -0.02462727 -0.00598565
## [2,] -0.6976955 -0.5423072581 -0.98211121 -0.23652145 -0.16120137
## [3,]  2.4367362 -0.0238218850 -0.32403419 -0.44206969 -0.47004393
## [4,] -2.4460385 -0.0007036966  0.08112764  0.14263545 -0.45584684
## [5,]  1.7708133  0.9741849574 -0.64716134  0.09377875 -0.78085072
## [6,] -0.8043943 -0.9469205017 -0.34705994 -0.62641753  0.26617649

The scatter plots for the score vectors specified by the argument v. The argument axes is specified by the two length vector represents which components are displayed.

par(mfrow=c(1,2))
plot(fit111, v="score", axes = 1:2, plottype="scatter")
plot(fit111, v="score", axes = 2:3, plottype="scatter")

When the argument v was specified as “cpev”, the cummulative eigenvalues are plotted.

par(mfrow=c(1,2))
plot(fit111, v="cpev", ylim=c(0.7, 1))

Other functions in R for PCA

There is the R function prcomp to implement PCA.

(fit1112 = prcomp(X1, scale=TRUE))

## Standard deviations (1, .., p=5):
## [1] 2.0446732 0.5899513 0.4458638 0.3926788 0.3439156
## 
## Rotation (n x k) = (5 x 5):
##             PC1         PC2         PC3           PC4         PC5
## [1,] -0.4309746  0.74172462 -0.03722419 -1.351296e-01  0.49442882
## [2,] -0.4483141 -0.31171881  0.63044575 -5.510830e-01 -0.02629751
## [3,] -0.4601629  0.19547669 -0.38616901 -1.416349e-01 -0.76213651
## [4,] -0.4392701 -0.55816074 -0.57114566  4.296727e-05  0.41144993
## [5,] -0.4566918 -0.05405032  0.35470924  8.111640e-01 -0.06859643

summary(fit1112)

## Importance of components:
##                           PC1     PC2     PC3     PC4     PC5
## Standard deviation     2.0447 0.58995 0.44586 0.39268 0.34392
## Proportion of Variance 0.8361 0.06961 0.03976 0.03084 0.02366
## Cumulative Proportion  0.8361 0.90575 0.94551 0.97634 1.00000

This Rotation is almost the same as the output of msma, but it can be made closer by setting the argument ceps as follows.

fit1113 = msma(X1, comp=5, ceps=0.0000001)
fit1113$wbX

## $block1
##           comp1       comp2       comp3       comp4       comp5
## X.1.1 0.4309745 -0.74172365 -0.03694568  0.13514153 -0.49444798
## X.1.2 0.4483141  0.31172789  0.63155896  0.54980542  0.02621947
## X.1.3 0.4601628 -0.19547781 -0.38588003  0.14252697  0.76211630
## X.1.4 0.4392701  0.55815708 -0.57114810  0.00105051 -0.41145010
## X.1.5 0.4566918  0.05404487  0.35306479 -0.81187175  0.06871165

Plotting the scores with the signs turned over, we see that similar scores are calculated.

par(mfrow=c(1,2))
biplot(fit1112)
plot(-fit1113$sbX[[1]][,1:2],xlab="Component 1",ylab="Component 2")

The ggfortify package is also available for the PCA plot.

Sparse PCA

If lambdaX (>0) is specified, a sparse principal component analysis is implemented.

(fit112 = msma(X1, comp=5, lambdaX=0.1))

## Call:
## msma.default(X = X1, comp = 5, lambdaX = 0.1)
## 
## Numbers of non-zeros for X block: 
##        comp1 comp2 comp3 comp4 comp5
## block1     5     4     4     5     4
## 
## Numbers of non-zeros for X super: 
##         comp1 comp2 comp3 comp4 comp5
## comp1-1     1     1     1     1     1

par(mfrow=c(1,2))
plot(fit112, axes = 1, plottype="bar", las=2)
plot(fit112, axes = 2, plottype="bar", las=2)

Supervised Sparse PCA

The outcome Z is generated.

set.seed(1); Z = rbinom(50, 1, 0.5)

If the outcome Z is specified, a supervised sparse principal component analysis is implemented.

(fit113 = msma(X1, Z=Z, comp=5, lambdaX=0.02))

## Call:
## msma.default(X = X1, Z = Z, comp = 5, lambdaX = 0.02)
## 
## Numbers of non-zeros for X block: 
##        comp1 comp2 comp3 comp4 comp5
## block1     5     5     5     4     5
## 
## Numbers of non-zeros for X super: 
##         comp1 comp2 comp3 comp4 comp5
## comp1-1     1     1     1     1     1

par(mfrow=c(1,2))
plot(fit113, axes = 1, plottype="bar", las=2)
plot(fit113, axes = 2, plottype="bar", las=2)

Partial Least Squres (PLS)

If the another input Y1 is specified, a partial least squres is implemented.

(fit121 = msma(X1, Y1, comp=2))

## Call:
## msma.default(X = X1, Y = Y1, comp = 2)
## 
## Numbers of non-zeros for X block: 
##        comp1 comp2
## block1     5     5
## 
## Numbers of non-zeros for X super: 
##         comp1 comp2
## comp1-1     1     1
## 
## Numbers of non-zeros for Y block: 
##        comp1 comp2
## block1     5     5
## 
## Numbers of non-zeros for Y super: 
##         comp1 comp2
## comp1-1     1     1

The component number is specified by the argument axes. When the argument XY was specified as “XY”, the scatter plots for Y score against X score are plotted.

par(mfrow=c(1,2))
plot(fit121, axes = 1, XY="XY")
plot(fit121, axes = 2, XY="XY")

Sparse PLS

If lambdaX and lambdaY are specified, a sparse PLS is implemented.

(fit122 = msma(X1, Y1, comp=2, lambdaX=0.5, lambdaY=0.5))

## Call:
## msma.default(X = X1, Y = Y1, comp = 2, lambdaX = 0.5, lambdaY = 0.5)
## 
## Numbers of non-zeros for X block: 
##        comp1 comp2
## block1     2     2
## 
## Numbers of non-zeros for X super: 
##         comp1 comp2
## comp1-1     1     1
## 
## Numbers of non-zeros for Y block: 
##        comp1 comp2
## block1     2     2
## 
## Numbers of non-zeros for Y super: 
##         comp1 comp2
## comp1-1     1     1

par(mfrow=c(1,2))
plot(fit122, axes = 1, XY="XY")
plot(fit122, axes = 2, XY="XY")

Supervised Sparse PLS

If the outcome Z is specified, a supervised sparse PLS is implemented.

(fit123 = msma(X1, Y1, Z, comp=2, lambdaX=0.5, lambdaY=0.5))

## Call:
## msma.default(X = X1, Y = Y1, Z = Z, comp = 2, lambdaX = 0.5, 
##     lambdaY = 0.5)
## 
## Numbers of non-zeros for X block: 
##        comp1 comp2
## block1     2     2
## 
## Numbers of non-zeros for X super: 
##         comp1 comp2
## comp1-1     1     1
## 
## Numbers of non-zeros for Y block: 
##        comp1 comp2
## block1     2     2
## 
## Numbers of non-zeros for Y super: 
##         comp1 comp2
## comp1-1     1     1

par(mfrow=c(1,2))
plot(fit123, axes = 1, XY="XY")
plot(fit123, axes = 2, XY="XY")

Multi Block

Multiblock data is a list of data matrix.

dataset2 = simdata(n = 50, rho = 0.8, Yps = c(2, 3), Xps = c(3, 4), seed=1)
X2 = dataset2$X; Y2 = dataset2$Y 

PCA

The input class is list.

class(X2)

## [1] "list"

The list length is 2 for 2 blocks.

length(X2)

## [1] 2

list of data matrix structure.

lapply(X2, dim)

## [[1]]
## [1] 50  3
## 
## [[2]]
## [1] 50  4

The function msma is applied to this list X2 as follows.

(fit211 = msma(X2, comp=1))

## Call:
## msma.default(X = X2, comp = 1)
## 
## Numbers of non-zeros for X block: 
##        comp1
## block1     3
## block2     4
## 
## Numbers of non-zeros for X super: 
##         comp1
## comp1-1     2

The bar plots for the block and super weights (loadings) specified the argument block.

par(mfrow=c(1,2))
plot(fit211, axes = 1, plottype="bar", block="block", las=2)
plot(fit211, axes = 1, plottype="bar", block="super")

Sparse PCA

If lambdaX with the length of 2 (same as the length of blocks) are specified, a multiblock sparse PCA is implemented.

(fit212 = msma(X2, comp=1, lambdaX=c(0.5, 0.5)))

## Call:
## msma.default(X = X2, comp = 1, lambdaX = c(0.5, 0.5))
## 
## Numbers of non-zeros for X block: 
##        comp1
## block1     3
## block2     2
## 
## Numbers of non-zeros for X super: 
##         comp1
## comp1-1     2

The bar plots for the block and super weights (loadings).

par(mfrow=c(1,2))
plot(fit212, axes = 1, plottype="bar", block="block", las=2)
plot(fit212, axes = 1, plottype="bar", block="super")

Supervised Sparse PCA

If the outcome Z is specified, a supervised analysis is implemented.

(fit213 = msma(X2, Z=Z, comp=1, lambdaX=c(0.5, 0.5)))

## Call:
## msma.default(X = X2, Z = Z, comp = 1, lambdaX = c(0.5, 0.5))
## 
## Numbers of non-zeros for X block: 
##        comp1
## block1     3
## block2     2
## 
## Numbers of non-zeros for X super: 
##         comp1
## comp1-1     2

par(mfrow=c(1,2))
plot(fit213, axes = 1, plottype="bar", block="block", las=2)
plot(fit213, axes = 1, plottype="bar", block="super")

Nested Component

A vector of length 2 can be given to the comp argument to perform the nested component analysis, which is a method to consider multiple components even in the super component. The first element of the vector corresponds to the number of block components and the second element corresponds to the number of (nested) super components.

(fit214 = msma(X2, comp=c(2,3)))

## Call:
## msma.default(X = X2, comp = c(2, 3))
## 
## Numbers of non-zeros for X block: 
##        comp1 comp2
## block1     3     3
## block2     4     4
## 
## Numbers of non-zeros for X super: 
## $comp1
## comp1-1 comp1-2 comp1-3 
##       2       2       2 
## 
## $comp2
## comp2-1 comp2-2 comp2-3 
##       2       2       2

In this example, there are 2 block components and 3 super components.

fit214$wbX

## $block1
##            comp1       comp2
## X.1.1 -0.5307011 -0.75618454
## X.1.2 -0.6006433  0.01688668
## X.1.3 -0.5979833  0.65414049
## 
## $block2
##           comp1         comp2
## X.2.1 0.4841242 -6.914840e-01
## X.2.2 0.5172474  5.175554e-05
## X.2.3 0.5241242  7.171312e-01
## X.2.4 0.4726233 -8.702155e-02

For the block weights, the number of blocks is 2 since there are two data matrices as shown as follows, and the number of rows is 3 and 4, the number of variables in each.

The number of components is 2 for the first element of the vector specified by the comp argument, which is the number of columns in each matrix.

par(mfrow=c(1,2))
plot(fit214, axes = 1, axes2 = 1, plottype="bar", block="block", las=2)
plot(fit214, axes = 2, axes2 = 1, plottype="bar", block="block", las=2)

fit214$wsX

## $comp1
##             comp1      comp2      comp3
## block1 -0.5198025 -0.8542864 -0.5242706
## block2  0.8542864 -0.5198025  0.8515517
## 
## $comp2
##             comp1      comp2      comp3
## block1 -0.3887274 -0.9213528  0.9899107
## block2  0.9213528 -0.3887274 -0.1416924

par(mfrow=c(2,3))
for(j in 1:2) for(i in 1:3) plot(fit214, axes = j, axes2 = i, plottype="bar", block="super")

PLS

If the another input (list) Y2 is specified, the partial least squared is implemented.

(fit221 = msma(X2, Y2, comp=1))

## Call:
## msma.default(X = X2, Y = Y2, comp = 1)
## 
## Numbers of non-zeros for X block: 
##        comp1
## block1     3
## block2     4
## 
## Numbers of non-zeros for X super: 
##         comp1
## comp1-1     2
## 
## Numbers of non-zeros for Y block: 
##        comp1
## block1     2
## block2     3
## 
## Numbers of non-zeros for Y super: 
##         comp1
## comp1-1     2

par(mfrow=c(1,2))
plot(fit221, axes = 1, plottype="bar", block="block", XY="X", las=2)
plot(fit221, axes = 1, plottype="bar", block="super", XY="X")

par(mfrow=c(1,2))
plot(fit221, axes = 1, plottype="bar", block="block", XY="Y", las=2)
plot(fit221, axes = 1, plottype="bar", block="super", XY="Y")

Sparse PLS

The regularized parameters lambdaX and lambdaY are specified vectors with same length with the length of lists X2 and Y2, respectively.

(fit222 = msma(X2, Y2, comp=1, lambdaX=c(0.5, 0.5), lambdaY=c(0.5, 0.5)))

## Call:
## msma.default(X = X2, Y = Y2, comp = 1, lambdaX = c(0.5, 0.5), 
##     lambdaY = c(0.5, 0.5))
## 
## Numbers of non-zeros for X block: 
##        comp1
## block1     2
## block2     2
## 
## Numbers of non-zeros for X super: 
##         comp1
## comp1-1     2
## 
## Numbers of non-zeros for Y block: 
##        comp1
## block1     1
## block2     3
## 
## Numbers of non-zeros for Y super: 
##         comp1
## comp1-1     2

par(mfrow=c(1,2))
plot(fit222, axes = 1, plottype="bar", block="block", XY="X", las=2)
plot(fit222, axes = 1, plottype="bar", block="super", XY="X")

par(mfrow=c(1,2))
plot(fit222, axes = 1, plottype="bar", block="block", XY="Y", las=2)
plot(fit222, axes = 1, plottype="bar", block="super", XY="Y")

Supervised Sparse PLS

(fit223 = msma(X2, Y2, Z, comp=1, lambdaX=c(0.5, 0.5), lambdaY=c(0.5, 0.5)))

## Call:
## msma.default(X = X2, Y = Y2, Z = Z, comp = 1, lambdaX = c(0.5, 
##     0.5), lambdaY = c(0.5, 0.5))
## 
## Numbers of non-zeros for X block: 
##        comp1
## block1     2
## block2     2
## 
## Numbers of non-zeros for X super: 
##         comp1
## comp1-1     2
## 
## Numbers of non-zeros for Y block: 
##        comp1
## block1     1
## block2     3
## 
## Numbers of non-zeros for Y super: 
##         comp1
## comp1-1     2

par(mfrow=c(1,2))
plot(fit223, axes = 1, plottype="bar", block="block", XY="X", las=2)
plot(fit223, axes = 1, plottype="bar", block="super", XY="X")

par(mfrow=c(1,2))
plot(fit223, axes = 1, plottype="bar", block="block", XY="Y", las=2)
plot(fit223, axes = 1, plottype="bar", block="super", XY="Y")

Parameter Selection

Number of Components

number of components search

Single Block

PCA

(ncomp11 = ncompsearch(X1, comps = c(1, 5, 10*(1:2)), nfold=5))

## $criterion
## [1] "CV"
## 
## $comps
## $comps[[1]]
## [1]  1  5 10 20
## 
## $comps[[2]]
## [1] 1
## 
## 
## $mincriterion
## [1] 1.044353e-31 1.000000e+00
## 
## $criterions
## $criterions[[1]]
## [1] 1.583224e-01 1.044353e-31 3.891936e-05 3.149335e-04
## 
## $criterions[[2]]
## [1] 1 1 1 1
## 
## 
## $optncomp
## [1] 5 1
## 
## $optlambdaX
## NULL
## 
## $optlambdaY
## NULL
## 
## $optlambdaXsup
## NULL
## 
## $optlambdaYsup
## NULL
## 
## attr(,"class")
## [1] "ncompsearch"

plot(ncomp11)

(ncomp12 = ncompsearch(X1, comps = 20, criterion="BIC"))

## $criterion
## [1] "BIC"
## 
## $comps
## $comps[[1]]
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20
## 
## $comps[[2]]
## [1] 1
## 
## 
## $mincriterion
##       comp5 comp1.comp1 
##   -70.84933   -70.49840 
## 
## $criterions
## $criterions$bic
##      comp1      comp2      comp3      comp4      comp5      comp6      comp7 
##  -1.718502  -2.161108  -2.598561  -3.322634 -70.849334 -65.470146 -65.368445 
##      comp8      comp9     comp10     comp11     comp12     comp13     comp14 
## -65.256004 -65.144626 -65.041593 -64.921986 -64.798407 -64.678425 -64.593260 
##     comp15     comp16     comp17     comp18     comp19     comp20 
## -64.468167 -64.374664 -64.275334 -64.138925 -64.029054 -63.917875 
## 
## $criterions$bic2
##  comp1.comp1  comp2.comp1  comp3.comp1  comp4.comp1  comp5.comp1  comp6.comp1 
##    -70.49840    -73.11129    -74.17235         -Inf    -73.81548    -69.52601 
##  comp7.comp1  comp8.comp1  comp9.comp1 comp10.comp1 comp11.comp1 comp12.comp1 
##    -65.20757    -65.12301    -65.06962    -64.96089    -64.88026    -64.81762 
## comp13.comp1 comp14.comp1 comp15.comp1 comp16.comp1 comp17.comp1 comp18.comp1 
##    -64.71127    -64.61512    -64.56091    -64.48897    -64.44835    -64.36345 
## comp19.comp1 comp20.comp1 
##    -64.26312    -64.19979 
## 
## 
## $optncomp
## [1] 5 1
## 
## $optlambdaX
## NULL
## 
## $optlambdaY
## NULL
## 
## $optlambdaXsup
## NULL
## 
## $optlambdaYsup
## NULL
## 
## attr(,"class")
## [1] "ncompsearch"

plot(ncomp12)

PLS

ncomp21 = ncompsearch(X2, Y2, comps = c(1, 5, 10*(1:2)), nfold=5)
plot(ncomp21)

Multi Block

The multi block structure has

dataset3 = simdata(n = 50, rho = 0.8, Yps = rep(4, 5), Xps = rep(4, 5), seed=1)
X3 = dataset3$X; Y3 = dataset3$Y 

PCA

ncomp31 = ncompsearch(X3, comps = 20, criterion="BIC")
plot(ncomp31)

Nested

(ncomp32 = ncompsearch(X3, comps = list(20, 20), criterion="BIC"))

## $criterion
## [1] "BIC"
## 
## $comps
## $comps[[1]]
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20
## 
## $comps[[2]]
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20
## 
## 
## $mincriterion
## $mincriterion$comp4
## [1] -71.35651
## 
## $mincriterion$comp5
##     comp1     comp2     comp3     comp4     comp5     comp6     comp7     comp8 
## -62.09665 -61.98622 -61.87579 -61.76537 -61.65494 -61.54451 -61.43408 -61.32365 
##     comp9    comp10    comp11    comp12    comp13    comp14    comp15    comp16 
## -61.21322 -61.10279 -60.99236 -60.88193 -60.77150 -60.66107 -60.55064 -60.44021 
##    comp17    comp18    comp19    comp20 
## -60.32979 -60.21936 -60.10893 -59.99850 
## 
## 
## $criterions
## $criterions$bic
##      comp1      comp2      comp3      comp4      comp5      comp6      comp7 
##  -1.825223  -2.241347  -2.958985 -71.356506 -65.458124 -65.308952 -65.171428 
##      comp8      comp9     comp10     comp11     comp12     comp13     comp14 
## -65.037531 -64.893623 -64.749608 -64.606181 -64.467683 -64.315116 -64.184393 
##     comp15     comp16     comp17     comp18     comp19     comp20 
## -64.028380 -63.896878 -63.762334 -63.624459 -63.495309 -63.365134 
## 
## $criterions$bic2
## $criterions$bic2$comp1
##       comp1       comp2       comp3       comp4       comp5       comp6 
##   0.9714158   0.7048689   0.3285174  -0.4043394 -70.1748081 -67.7655048 
##       comp7       comp8       comp9      comp10      comp11      comp12 
## -67.5712460 -67.7718902 -67.5677851 -67.6442733 -67.4184431 -67.3674295 
##      comp13      comp14      comp15      comp16      comp17      comp18 
## -67.1229968 -66.8735021 -66.7556268 -66.6656846 -66.6941339 -66.5614065 
##      comp19      comp20 
## -66.3343175 -66.1795437 
## 
## $criterions$bic2$comp2
##       comp1       comp2       comp3       comp4       comp5       comp6 
##   0.3896461   0.0920536  -0.2677686  -0.9695676 -69.9271549 -67.4334486 
##       comp7       comp8       comp9      comp10      comp11      comp12 
## -67.3191452 -67.1180722 -67.1080342 -67.0136077 -66.8514965 -66.6046869 
##      comp13      comp14      comp15      comp16      comp17      comp18 
## -66.4873417 -66.4897346 -66.3233707 -66.2485745 -66.1583394 -65.8363530 
##      comp19      comp20 
## -65.7566891 -65.6339959 
## 
## $criterions$bic2$comp3
##       comp1       comp2       comp3       comp4       comp5       comp6 
##   2.4279355   2.1068041   1.6379884   0.8785787 -67.6880762 -65.1878125 
##       comp7       comp8       comp9      comp10      comp11      comp12 
## -65.0691299 -64.8096419 -64.8433985 -64.6429184 -64.5593624 -64.3852060 
##      comp13      comp14      comp15      comp16      comp17      comp18 
## -64.2765664 -64.1790566 -64.1007365 -63.9825735 -63.8279595 -63.7506078 
##      comp19      comp20 
## -63.4982199 -63.2869695 
## 
## $criterions$bic2$comp4
##      comp1      comp2      comp3      comp4      comp5      comp6      comp7 
##   4.356252   4.011048   3.618331   2.814239 -66.116523 -62.885107 -62.764886 
##      comp8      comp9     comp10     comp11     comp12     comp13     comp14 
## -62.563141 -62.514718 -62.424698 -62.242041 -61.990384 -61.887556 -61.767548 
##     comp15     comp16     comp17     comp18     comp19     comp20 
## -61.640053 -61.494525 -61.375280 -61.304241 -61.128873 -61.074902 
## 
## $criterions$bic2$comp5
##     comp1     comp2     comp3     comp4     comp5     comp6     comp7     comp8 
## -62.09665 -61.98622 -61.87579 -61.76537 -61.65494 -61.54451 -61.43408 -61.32365 
##     comp9    comp10    comp11    comp12    comp13    comp14    comp15    comp16 
## -61.21322 -61.10279 -60.99236 -60.88193 -60.77150 -60.66107 -60.55064 -60.44021 
##    comp17    comp18    comp19    comp20 
## -60.32979 -60.21936 -60.10893 -59.99850 
## 
## $criterions$bic2$comp6
##     comp1     comp2     comp3     comp4     comp5     comp6     comp7     comp8 
## -54.50030 -54.38987 -54.27944 -54.16902 -54.05859 -53.94816 -53.83773 -53.72730 
##     comp9    comp10    comp11    comp12    comp13    comp14    comp15    comp16 
## -53.61687 -53.50644 -53.39601 -53.28558 -53.17515 -53.06472 -52.95429 -52.84386 
##    comp17    comp18    comp19    comp20 
## -52.73344 -52.62301 -52.51258 -52.40215 
## 
## $criterions$bic2$comp7
##     comp1     comp2     comp3     comp4     comp5     comp6     comp7     comp8 
## -52.26402 -52.15359 -52.04316 -51.93274 -51.82231 -51.71188 -51.60145 -51.49102 
##     comp9    comp10    comp11    comp12    comp13    comp14    comp15    comp16 
## -51.38059 -51.27016 -51.15973 -51.04930 -50.93887 -50.82844 -50.71801 -50.60759 
##    comp17    comp18    comp19    comp20 
## -50.49716 -50.38673 -50.27630 -50.16587 
## 
## $criterions$bic2$comp8
##     comp1     comp2     comp3     comp4     comp5     comp6     comp7     comp8 
## -50.05248 -49.94205 -49.83162 -49.72119 -49.61076 -49.50033 -49.38990 -49.27948 
##     comp9    comp10    comp11    comp12    comp13    comp14    comp15    comp16 
## -49.16905 -49.05862 -48.94819 -48.83776 -48.72733 -48.61690 -48.50647 -48.39604 
##    comp17    comp18    comp19    comp20 
## -48.28561 -48.17518 -48.06475 -47.95432 
## 
## $criterions$bic2$comp9
##     comp1     comp2     comp3     comp4     comp5     comp6     comp7     comp8 
## -47.85680 -47.74637 -47.63594 -47.52551 -47.41508 -47.30465 -47.19422 -47.08379 
##     comp9    comp10    comp11    comp12    comp13    comp14    comp15    comp16 
## -46.97337 -46.86294 -46.75251 -46.64208 -46.53165 -46.42122 -46.31079 -46.20036 
##    comp17    comp18    comp19    comp20 
## -46.08993 -45.97950 -45.86907 -45.75864 
## 
## $criterions$bic2$comp10
##     comp1     comp2     comp3     comp4     comp5     comp6     comp7     comp8 
## -45.66784 -45.55741 -45.44698 -45.33655 -45.22612 -45.11569 -45.00526 -44.89483 
##     comp9    comp10    comp11    comp12    comp13    comp14    comp15    comp16 
## -44.78440 -44.67398 -44.56355 -44.45312 -44.34269 -44.23226 -44.12183 -44.01140 
##    comp17    comp18    comp19    comp20 
## -43.90097 -43.79054 -43.68011 -43.56968 
## 
## $criterions$bic2$comp11
##     comp1     comp2     comp3     comp4     comp5     comp6     comp7     comp8 
## -43.44288 -43.33245 -43.22202 -43.11159 -43.00116 -42.89073 -42.78030 -42.66988 
##     comp9    comp10    comp11    comp12    comp13    comp14    comp15    comp16 
## -42.55945 -42.44902 -42.33859 -42.22816 -42.11773 -42.00730 -41.89687 -41.78644 
##    comp17    comp18    comp19    comp20 
## -41.67601 -41.56558 -41.45515 -41.34472 
## 
## $criterions$bic2$comp12
##     comp1     comp2     comp3     comp4     comp5     comp6     comp7     comp8 
## -41.21638 -41.10595 -40.99552 -40.88510 -40.77467 -40.66424 -40.55381 -40.44338 
##     comp9    comp10    comp11    comp12    comp13    comp14    comp15    comp16 
## -40.33295 -40.22252 -40.11209 -40.00166 -39.89123 -39.78080 -39.67037 -39.55994 
##    comp17    comp18    comp19    comp20 
## -39.44952 -39.33909 -39.22866 -39.11823 
## 
## $criterions$bic2$comp13
##     comp1     comp2     comp3     comp4     comp5     comp6     comp7     comp8 
## -39.01504 -38.90461 -38.79418 -38.68375 -38.57332 -38.46289 -38.35246 -38.24203 
##     comp9    comp10    comp11    comp12    comp13    comp14    comp15    comp16 
## -38.13161 -38.02118 -37.91075 -37.80032 -37.68989 -37.57946 -37.46903 -37.35860 
##    comp17    comp18    comp19    comp20 
## -37.24817 -37.13774 -37.02731 -36.91688 
## 
## $criterions$bic2$comp14
##     comp1     comp2     comp3     comp4     comp5     comp6     comp7     comp8 
## -36.80469 -36.69426 -36.58383 -36.47340 -36.36297 -36.25254 -36.14212 -36.03169 
##     comp9    comp10    comp11    comp12    comp13    comp14    comp15    comp16 
## -35.92126 -35.81083 -35.70040 -35.58997 -35.47954 -35.36911 -35.25868 -35.14825 
##    comp17    comp18    comp19    comp20 
## -35.03782 -34.92739 -34.81696 -34.70654 
## 
## $criterions$bic2$comp15
##     comp1     comp2     comp3     comp4     comp5     comp6     comp7     comp8 
## -34.61125 -34.50082 -34.39039 -34.27996 -34.16953 -34.05910 -33.94867 -33.83824 
##     comp9    comp10    comp11    comp12    comp13    comp14    comp15    comp16 
## -33.72781 -33.61739 -33.50696 -33.39653 -33.28610 -33.17567 -33.06524 -32.95481 
##    comp17    comp18    comp19    comp20 
## -32.84438 -32.73395 -32.62352 -32.51309 
## 
## $criterions$bic2$comp16
##     comp1     comp2     comp3     comp4     comp5     comp6     comp7     comp8 
## -32.39008 -32.27965 -32.16923 -32.05880 -31.94837 -31.83794 -31.72751 -31.61708 
##     comp9    comp10    comp11    comp12    comp13    comp14    comp15    comp16 
## -31.50665 -31.39622 -31.28579 -31.17536 -31.06493 -30.95450 -30.84407 -30.73365 
##    comp17    comp18    comp19    comp20 
## -30.62322 -30.51279 -30.40236 -30.29193 
## 
## $criterions$bic2$comp17
##     comp1     comp2     comp3     comp4     comp5     comp6     comp7     comp8 
## -30.19018 -30.07975 -29.96932 -29.85889 -29.74846 -29.63803 -29.52760 -29.41717 
##     comp9    comp10    comp11    comp12    comp13    comp14    comp15    comp16 
## -29.30675 -29.19632 -29.08589 -28.97546 -28.86503 -28.75460 -28.64417 -28.53374 
##    comp17    comp18    comp19    comp20 
## -28.42331 -28.31288 -28.20245 -28.09202 
## 
## $criterions$bic2$comp18
##     comp1     comp2     comp3     comp4     comp5     comp6     comp7     comp8 
## -27.98137 -27.87094 -27.76051 -27.65008 -27.53965 -27.42922 -27.31879 -27.20836 
##     comp9    comp10    comp11    comp12    comp13    comp14    comp15    comp16 
## -27.09794 -26.98751 -26.87708 -26.76665 -26.65622 -26.54579 -26.43536 -26.32493 
##    comp17    comp18    comp19    comp20 
## -26.21450 -26.10407 -25.99364 -25.88321 
## 
## $criterions$bic2$comp19
##     comp1     comp2     comp3     comp4     comp5     comp6     comp7     comp8 
## -25.81292 -25.70249 -25.59206 -25.48163 -25.37120 -25.26077 -25.15034 -25.03991 
##     comp9    comp10    comp11    comp12    comp13    comp14    comp15    comp16 
## -24.92948 -24.81906 -24.70863 -24.59820 -24.48777 -24.37734 -24.26691 -24.15648 
##    comp17    comp18    comp19    comp20 
## -24.04605 -23.93562 -23.82519 -23.71476 
## 
## $criterions$bic2$comp20
##     comp1     comp2     comp3     comp4     comp5     comp6     comp7     comp8 
## -23.56844 -23.45801 -23.34758 -23.23715 -23.12672 -23.01629 -22.90586 -22.79543 
##     comp9    comp10    comp11    comp12    comp13    comp14    comp15    comp16 
## -22.68500 -22.57457 -22.46415 -22.35372 -22.24329 -22.13286 -22.02243 -21.91200 
##    comp17    comp18    comp19    comp20 
## -21.80157 -21.69114 -21.58071 -21.47028 
## 
## 
## 
## $optncomp
## [1] 4 5
## 
## $optlambdaX
## NULL
## 
## $optlambdaY
## NULL
## 
## $optlambdaXsup
## NULL
## 
## $optlambdaYsup
## NULL
## 
## attr(,"class")
## [1] "ncompsearch"

par(mfrow=c(1,2))
plot(ncomp32,1)
plot(ncomp32,2)

PLS

ncomp41 = ncompsearch(X3, Y3, comps = c(1, 5, 10*(1:2)), criterion="BIC")
plot(ncomp41)

Combined Search

The number of components and regularized parameters can be selected by the function optparasearch. The following options are available.

criteria = c("BIC", "CV")
search.methods = c("regparaonly", "regpara1st", "ncomp1st", "simultaneous")

Single Block

PCA

• The regparaonly method searches for the regularized parameters with a fixed number of components.

(opt11 = optparasearch(X1, search.method = "regparaonly", criterion="BIC", comp=ncomp11$optncomp))

## $optncomp
## [1] 5 1
## 
## $optlambdaX
## [1] 0.4011227
## 
## $search.method
## [1] "regparaonly"
## 
## $criterion
## [1] "BIC"
## 
## $criterion4ncomp
## [1] "BIC"
## 
## attr(,"class")
## [1] "optparasearch"

(fit311 = msma(X1, comp=opt11$optncomp, lambdaX=opt11$optlambdaX))

## Call:
## msma.default(X = X1, comp = opt11$optncomp, lambdaX = opt11$optlambdaX)
## 
## Numbers of non-zeros for X block: 
##        comp1 comp2 comp3 comp4 comp5
## block1     5     2     1     2     2
## 
## Numbers of non-zeros for X super: 
##         comp1 comp2 comp3 comp4 comp5
## comp1-1     1     1     1     1     1

• The ncomp1st method identifies the number of components with a regularized parameter of 0, then searches for the regularized parameters with the selected number of components.

(opt12 = optparasearch(X1, search.method = "ncomp1st", criterion="BIC"))

## $criterion
## [1] "BIC"
## 
## $comps
## $comps[[1]]
##  [1]  1  2  3  4  5  6  7  8  9 10
## 
## $comps[[2]]
## [1] 1
## 
## 
## $mincriterion
##       comp5 comp1.comp1 
##   -70.84933   -70.49840 
## 
## $criterions
## $criterions$bic
##      comp1      comp2      comp3      comp4      comp5      comp6      comp7 
##  -1.718502  -2.161108  -2.598561  -3.322634 -70.849334 -65.470146 -65.368445 
##      comp8      comp9     comp10 
## -65.256004 -65.144626 -65.041593 
## 
## $criterions$bic2
##  comp1.comp1  comp2.comp1  comp3.comp1  comp4.comp1  comp5.comp1  comp6.comp1 
##    -70.49840    -73.11129    -74.17235         -Inf    -73.81548    -69.52601 
##  comp7.comp1  comp8.comp1  comp9.comp1 comp10.comp1 
##    -65.20757    -65.12301    -65.06962    -64.96089 
## 
## 
## $optncomp
## [1] 5 1
## 
## $optlambdaX
## [1] 0.4011227
## 
## $search.method
## [1] "ncomp1st"
## 
## $criterion4ncomp
## [1] "BIC"
## 
## attr(,"class")
## [1] "optparasearch"

(fit312 = msma(X1, comp=opt12$optncomp, lambdaX=opt12$optlambdaX))

## Call:
## msma.default(X = X1, comp = opt12$optncomp, lambdaX = opt12$optlambdaX)
## 
## Numbers of non-zeros for X block: 
##        comp1 comp2 comp3 comp4 comp5
## block1     5     2     1     2     2
## 
## Numbers of non-zeros for X super: 
##         comp1 comp2 comp3 comp4 comp5
## comp1-1     1     1     1     1     1

• The regpara1st identifies the regularized parameters by fixing the number of components, then searching for the number of components with the selected regularized parameters.

(opt13 = optparasearch(X1, search.method = "regpara1st", criterion="BIC"))

## $criterion
## [1] "BIC"
## 
## $comps
## $comps[[1]]
##  [1]  1  2  3  4  5  6  7  8  9 10
## 
## $comps[[2]]
## [1] 1
## 
## 
## $mincriterion
##       comp5 comp1.comp1 
##   -71.23845   -69.27306 
## 
## $criterions
## $criterions$bic
##      comp1      comp2      comp3      comp4      comp5      comp6      comp7 
##  -1.718325  -2.202584  -2.657239  -3.466906 -71.238453 -65.620115 -65.559015 
##      comp8      comp9     comp10 
## -65.493668 -65.444509 -65.403246 
## 
## $criterions$bic2
##  comp1.comp1  comp2.comp1  comp3.comp1  comp4.comp1  comp5.comp1  comp6.comp1 
##    -69.27306    -73.16344    -74.32756    -73.79695    -73.78251    -69.90172 
##  comp7.comp1  comp8.comp1  comp9.comp1 comp10.comp1 
##    -65.27116    -65.20281    -65.10658    -65.02588 
## 
## 
## $optncomp
## [1] 5 1
## 
## $optlambdaX
## [1] 0.2578646
## 
## $optlambdaY
## NULL
## 
## $optlambdaXsup
## NULL
## 
## $optlambdaYsup
## NULL
## 
## $search.method
## [1] "regpara1st"
## 
## $criterion4ncomp
## [1] "BIC"
## 
## attr(,"class")
## [1] "optparasearch"

(fit313 = msma(X1, comp=opt13$optncomp, lambdaX=opt13$optlambdaX))

## Call:
## msma.default(X = X1, comp = opt13$optncomp, lambdaX = opt13$optlambdaX)
## 
## Numbers of non-zeros for X block: 
##        comp1 comp2 comp3 comp4 comp5
## block1     5     2     3     2     3
## 
## Numbers of non-zeros for X super: 
##         comp1 comp2 comp3 comp4 comp5
## comp1-1     1     1     1     1     1

• The simultaneous method identifies the number of components by searching the regularized parameters in each component.

(opt14 = optparasearch(X1, search.method = "simultaneous", criterion="BIC"))

## $criterion
## [1] "BIC"
## 
## $comps
## $comps[[1]]
##  [1]  1  2  3  4  5  6  7  8  9 10
## 
## $comps[[2]]
## [1] 1
## 
## 
## $mincriterion
## cve.criterion       lambdaX 
##  -72.08886650    0.02865162 
## 
## $criterions
## $criterions[[1]]
## cve.criterion cve.criterion cve.criterion cve.criterion cve.criterion 
##     -1.718502     -2.202584     -2.677000     -3.497959    -72.088866 
## cve.criterion cve.criterion cve.criterion cve.criterion cve.criterion 
##    -65.571063    -65.463953    -65.493668    -65.444509    -65.403246 
## 
## $criterions[[2]]
##    lambdaX    lambdaX    lambdaX    lambdaX    lambdaX    lambdaX    lambdaX 
## 0.02865162 0.25786460 0.34381947 0.34381947 0.40112271 0.08595487 0.08595487 
##    lambdaX    lambdaX    lambdaX 
## 0.25786460 0.25786460 0.25786460 
## 
## 
## $optncomp
## [1] 5 1
## 
## $optlambdaX
## [1] 0.4011227
## 
## $optlambdaY
## NULL
## 
## $optlambdaXsup
## NULL
## 
## $optlambdaYsup
## NULL
## 
## $search.method
## [1] "simultaneous"
## 
## $criterion4ncomp
## [1] "BIC"
## 
## attr(,"class")
## [1] "optparasearch"

(fit314 = msma(X1, comp=opt14$optncomp, lambdaX=opt14$optlambdaX))

## Call:
## msma.default(X = X1, comp = opt14$optncomp, lambdaX = opt14$optlambdaX)
## 
## Numbers of non-zeros for X block: 
##        comp1 comp2 comp3 comp4 comp5
## block1     5     2     1     2     2
## 
## Numbers of non-zeros for X super: 
##         comp1 comp2 comp3 comp4 comp5
## comp1-1     1     1     1     1     1

The argument maxpct4ncomp=0.5 means that 0.5\(\lambda\) is used as the regularized parameter when the number of components is searched and where \(\lambda\) is the maximum of the regularized parameters among the possible candidates.

(opt132 = optparasearch(X1, search.method = "ncomp1st", criterion="BIC", maxpct4ncomp=0.5))

## $criterion
## [1] "BIC"
## 
## $comps
## $comps[[1]]
##  [1]  1  2  3  4  5  6  7  8  9 10
## 
## $comps[[2]]
## [1] 1
## 
## 
## $mincriterion
## cve.criterion       lambdaX 
##  -72.40872188    0.01432581 
## 
## $criterions
## $criterions[[1]]
## cve.criterion cve.criterion cve.criterion cve.criterion cve.criterion 
##     -1.718502     -2.190527     -2.649988     -3.450188    -72.408722 
## cve.criterion cve.criterion cve.criterion cve.criterion cve.criterion 
##    -65.571063    -65.463953    -65.426285    -65.326271    -65.239644 
## 
## $criterions[[2]]
##    lambdaX    lambdaX    lambdaX    lambdaX    lambdaX    lambdaX    lambdaX 
## 0.01432581 0.17190973 0.17190973 0.18623555 0.17190973 0.08595487 0.08595487 
##    lambdaX    lambdaX    lambdaX 
## 0.05730324 0.05730324 0.05730324 
## 
## 
## $optncomp
## [1] 5 1
## 
## $optlambdaX
## [1] 0.4011227
## 
## $search.method
## [1] "ncomp1st"
## 
## $criterion4ncomp
## [1] "BIC"
## 
## attr(,"class")
## [1] "optparasearch"

(fit3132 = msma(X1, comp=opt132$optncomp, lambdaX=opt132$optlambdaX))

## Call:
## msma.default(X = X1, comp = opt132$optncomp, lambdaX = opt132$optlambdaX)
## 
## Numbers of non-zeros for X block: 
##        comp1 comp2 comp3 comp4 comp5
## block1     5     2     1     2     2
## 
## Numbers of non-zeros for X super: 
##         comp1 comp2 comp3 comp4 comp5
## comp1-1     1     1     1     1     1

The result with the argument regpara1st depends on the number of components and the default value is 10. The number of components is set as follows.

(opt133 = optparasearch(X1, search.method = "regpara1st", criterion="BIC", comp=5))

## $criterion
## [1] "BIC"
## 
## $comps
## $comps[[1]]
## [1] 1 2 3 4 5
## 
## $comps[[2]]
## [1] 1
## 
## 
## $mincriterion
##       comp5 comp1.comp1 
##   -72.08887   -70.73420 
## 
## $criterions
## $criterions$bic
##      comp1      comp2      comp3      comp4      comp5 
##  -1.711131  -2.183698  -2.634114  -3.400602 -72.088866 
## 
## $criterions$bic2
## comp1.comp1 comp2.comp1 comp3.comp1 comp4.comp1 comp5.comp1 
##   -70.73420   -73.22687   -72.43315   -73.85812   -73.84430 
## 
## 
## $optncomp
## [1] 5 1
## 
## $optlambdaX
## [1] 0.4011227
## 
## $optlambdaY
## NULL
## 
## $optlambdaXsup
## NULL
## 
## $optlambdaYsup
## NULL
## 
## $search.method
## [1] "regpara1st"
## 
## $criterion4ncomp
## [1] "BIC"
## 
## attr(,"class")
## [1] "optparasearch"

(fit3133 = msma(X1, comp=opt133$optncomp, lambdaX=opt133$optlambdaX))

## Call:
## msma.default(X = X1, comp = opt133$optncomp, lambdaX = opt133$optlambdaX)
## 
## Numbers of non-zeros for X block: 
##        comp1 comp2 comp3 comp4 comp5
## block1     5     2     1     2     2
## 
## Numbers of non-zeros for X super: 
##         comp1 comp2 comp3 comp4 comp5
## comp1-1     1     1     1     1     1

PLS

For PLS, two parameters \(\lambda_X\) and \(\lambda_Y\) are used in arguments lambdaX and lambdaY to control sparseness for data X and Y, respectively.

(opt21 = optparasearch(X2, Y2, search.method = "regparaonly", criterion="BIC"))

## $optncomp
## [1] 10  1
## 
## $optlambdaX
##   lambdaX1   lambdaX2 
## 0.08245431 0.48828296 
## 
## $optlambdaY
##  lambdaY1  lambdaY2 
## 0.4318262 0.2256272 
## 
## $optlambdaXsup
## lambdaXsup1 
##   0.4273084 
## 
## $optlambdaYsup
## lambdaYsup1 
##   0.6187171 
## 
## $search.method
## [1] "regparaonly"
## 
## $criterion
## [1] "BIC"
## 
## $criterion4ncomp
## [1] "BIC"
## 
## attr(,"class")
## [1] "optparasearch"

(fit321 = msma(X2, Y2, comp=opt21$optncomp, lambdaX=opt21$optlambdaX, lambdaY=opt21$optlambdaY))

## Call:
## msma.default(X = X2, Y = Y2, comp = opt21$optncomp, lambdaX = opt21$optlambdaX, 
##     lambdaY = opt21$optlambdaY)
## 
## Numbers of non-zeros for X block: 
##        comp1 comp2 comp3 comp4 comp5 comp6 comp7 comp8 comp9 comp10
## block1     3     3     3     2     3     3     2     3     3      3
## block2     2     2     1     2     2     1     1     1     2      1
## 
## Numbers of non-zeros for X super: 
##         comp1 comp2 comp3 comp4 comp5 comp6 comp7 comp8 comp9 comp10
## comp1-1     2     2     2     2     2     2     2     2     2      2
## 
## Numbers of non-zeros for Y block: 
##        comp1 comp2 comp3 comp4 comp5 comp6 comp7 comp8 comp9 comp10
## block1     1     1     2     1     1     1     1     1     1      1
## block2     3     2     3     2     2     2     3     3     2      2
## 
## Numbers of non-zeros for Y super: 
##         comp1 comp2 comp3 comp4 comp5 comp6 comp7 comp8 comp9 comp10
## comp1-1     2     2     2     2     2     2     2     2     2      2

Multi Block

PCA

(opt31 = optparasearch(X3, search.method = "regparaonly", criterion="BIC"))

## $optncomp
## [1] 10  1
## 
## $optlambdaX
##  lambdaX1  lambdaX2  lambdaX3  lambdaX4  lambdaX5 
## 0.1310995 0.1513608 0.1590729 0.1300171 0.3984990 
## 
## $optlambdaXsup
## lambdaXsup1 
##   0.1493833 
## 
## $search.method
## [1] "regparaonly"
## 
## $criterion
## [1] "BIC"
## 
## $criterion4ncomp
## [1] "BIC"
## 
## attr(,"class")
## [1] "optparasearch"

(fit331 = msma(X3, comp=opt31$optncomp, lambdaX=opt31$optlambdaX, lambdaXsup=opt31$optlambdaXsup))

## Call:
## msma.default(X = X3, comp = opt31$optncomp, lambdaX = opt31$optlambdaX, 
##     lambdaXsup = opt31$optlambdaXsup)
## 
## Numbers of non-zeros for X block: 
##        comp1 comp2 comp3 comp4 comp5 comp6 comp7 comp8 comp9 comp10
## block1     0     4     0     2     3     0     0     3     0      3
## block2     0     0     2     2     0     0     0     2     2      0
## block3     4     0     2     0     0     4     0     0     4      4
## block4     0     2     2     3     0     2     0     0     0      0
## block5     0     0     2     2     0     2     1     0     0      1
## 
## Numbers of non-zeros for X super: 
##         comp1 comp2 comp3 comp4 comp5 comp6 comp7 comp8 comp9 comp10
## comp1-1     1     2     4     4     1     3     1     2     2      3

(opt32 = optparasearch(X3, search.method = "regparaonly", criterion="BIC", whichselect="X"))

## $optncomp
## [1] 10  1
## 
## $optlambdaX
##   lambdaX1   lambdaX2   lambdaX3   lambdaX4   lambdaX5 
## 0.42607353 0.56760297 0.11930469 0.48756415 0.09962475 
## 
## $search.method
## [1] "regparaonly"
## 
## $criterion
## [1] "BIC"
## 
## $criterion4ncomp
## [1] "BIC"
## 
## attr(,"class")
## [1] "optparasearch"

(fit332 = msma(X3, comp=opt32$optncomp, lambdaX=opt32$optlambdaX, lambdaXsup=opt32$optlambdaXsup))

## Call:
## msma.default(X = X3, comp = opt32$optncomp, lambdaX = opt32$optlambdaX, 
##     lambdaXsup = opt32$optlambdaXsup)
## 
## Numbers of non-zeros for X block: 
##        comp1 comp2 comp3 comp4 comp5 comp6 comp7 comp8 comp9 comp10
## block1     2     2     1     1     2     1     1     1     1      1
## block2     2     1     1     1     1     1     1     1     1      1
## block3     4     4     3     3     3     4     4     4     4      4
## block4     1     2     1     2     1     1     1     1     1      1
## block5     4     3     3     4     2     3     2     2     3      3
## 
## Numbers of non-zeros for X super: 
##         comp1 comp2 comp3 comp4 comp5 comp6 comp7 comp8 comp9 comp10
## comp1-1     5     5     5     5     5     5     5     5     5      5

(opt33 = optparasearch(X3, search.method = "regparaonly", criterion="BIC", whichselect="Xsup"))

## $optncomp
## [1] 10  1
## 
## $optlambdaXsup
## [1] 0.1493833
## 
## $search.method
## [1] "regparaonly"
## 
## $criterion
## [1] "BIC"
## 
## $criterion4ncomp
## [1] "BIC"
## 
## attr(,"class")
## [1] "optparasearch"

(fit333 = msma(X3, comp=opt33$optncomp, lambdaX=opt33$optlambdaX, lambdaXsup=opt33$optlambdaXsup))

## Call:
## msma.default(X = X3, comp = opt33$optncomp, lambdaX = opt33$optlambdaX, 
##     lambdaXsup = opt33$optlambdaXsup)
## 
## Numbers of non-zeros for X block: 
##        comp1 comp2 comp3 comp4 comp5 comp6 comp7 comp8 comp9 comp10
## block1     0     0     3     0     2     0     0     1     2      0
## block2     0     0     4     3     0     4     0     0     0      2
## block3     4     0     0     4     3     0     2     0     0      0
## block4     0     4     3     3     4     0     0     0     0      0
## block5     4     0     4     2     0     0     2     0     0      0
## 
## Numbers of non-zeros for X super: 
##         comp1 comp2 comp3 comp4 comp5 comp6 comp7 comp8 comp9 comp10
## comp1-1     2     1     4     4     3     1     2     1     1      1

Nested

ncomp1st

(opt341 = optparasearch(X3, search.method = "ncomp1st", criterion="BIC", comp=c(8, 8)))

## $criterion
## [1] "BIC"
## 
## $comps
## $comps[[1]]
## [1] 1 2 3 4 5 6 7 8
## 
## $comps[[2]]
## [1] 1 2 3 4 5 6 7 8
## 
## 
## $mincriterion
## $mincriterion$comp4
## [1] -71.35651
## 
## $mincriterion$comp5
##     comp1     comp2     comp3     comp4     comp5     comp6     comp7     comp8 
## -67.39726 -67.28683 -67.17640 -67.06597 -66.95554 -66.84511 -66.73468 -66.62425 
## 
## 
## $criterions
## $criterions$bic
##      comp1      comp2      comp3      comp4      comp5      comp6      comp7 
##  -1.825223  -2.241347  -2.958985 -71.356506 -65.458124 -65.308952 -65.171428 
##      comp8 
## -65.037531 
## 
## $criterions$bic2
## $criterions$bic2$comp1
##       comp1       comp2       comp3       comp4       comp5       comp6 
##   0.9714158   0.7048689   0.3285174  -0.4043394 -70.1748081 -67.7655048 
##       comp7       comp8 
## -67.5712460 -67.7718902 
## 
## $criterions$bic2$comp2
##       comp1       comp2       comp3       comp4       comp5       comp6 
##  -0.9355045  -1.2330970  -1.5929192  -2.2947183 -71.2523055 -68.7585992 
##       comp7       comp8 
## -68.6442958 -68.4432228 
## 
## $criterions$bic2$comp3
##       comp1       comp2       comp3       comp4       comp5       comp6 
##  -0.2223657  -0.5434971  -1.0123128  -1.7717225 -70.3383775 -67.8381137 
##       comp7       comp8 
## -67.7194311 -67.4599431 
## 
## $criterions$bic2$comp4
##        comp1        comp2        comp3        comp4        comp5        comp6 
##   0.38080030   0.03559579  -0.35712052  -1.16121248 -70.09197475 -66.86055880 
##        comp7        comp8 
## -66.74033788 -66.53859245 
## 
## $criterions$bic2$comp5
##     comp1     comp2     comp3     comp4     comp5     comp6     comp7     comp8 
## -67.39726 -67.28683 -67.17640 -67.06597 -66.95554 -66.84511 -66.73468 -66.62425 
## 
## $criterions$bic2$comp6
##     comp1     comp2     comp3     comp4     comp5     comp6     comp7     comp8 
## -61.12606 -61.01563 -60.90520 -60.79477 -60.68434 -60.57391 -60.46348 -60.35305 
## 
## $criterions$bic2$comp7
##     comp1     comp2     comp3     comp4     comp5     comp6     comp7     comp8 
## -60.21493 -60.10450 -59.99407 -59.88364 -59.77321 -59.66278 -59.55235 -59.44192 
## 
## $criterions$bic2$comp8
##     comp1     comp2     comp3     comp4     comp5     comp6     comp7     comp8 
## -59.32853 -59.21810 -59.10768 -58.99725 -58.88682 -58.77639 -58.66596 -58.55553 
## 
## 
## 
## $optncomp
## [1] 4 5
## 
## $optlambdaX
##  lambdaX1  lambdaX2  lambdaX3  lambdaX4  lambdaX5 
## 0.4916233 0.0378402 0.1193047 0.4875641 0.2324578 
## 
## $optlambdaXsup
## lambdaXsup1 
##  0.03734583 
## 
## $search.method
## [1] "ncomp1st"
## 
## $criterion4ncomp
## [1] "BIC"
## 
## attr(,"class")
## [1] "optparasearch"

(fit341 = msma(X3, comp=opt341$optncomp, lambdaX=opt341$optlambdaX, lambdaXsup=opt341$optlambdaXsup))

## Call:
## msma.default(X = X3, comp = opt341$optncomp, lambdaX = opt341$optlambdaX, 
##     lambdaXsup = opt341$optlambdaXsup)
## 
## Numbers of non-zeros for X block: 
##        comp1 comp2 comp3 comp4
## block1     2     2     1     1
## block2     4     4     4     2
## block3     4     3     3     4
## block4     1     2     1     2
## block5     4     3     2     3
## 
## Numbers of non-zeros for X super: 
## $comp1
## comp1-1 comp1-2 comp1-3 comp1-4 comp1-5 
##       5       3       4       3       4 
## 
## $comp2
## comp2-1 comp2-2 comp2-3 comp2-4 comp2-5 
##       5       5       5       4       5 
## 
## $comp3
## comp3-1 comp3-2 comp3-3 comp3-4 comp3-5 
##       5       5       4       5       5 
## 
## $comp4
## comp4-1 comp4-2 comp4-3 comp4-4 comp4-5 
##       5       4       4       5       5

regparaonly

(opt342 = optparasearch(X3, search.method = "regparaonly", criterion="BIC", comp=c(4, 5)))

## $optncomp
## [1] 4 5
## 
## $optlambdaX
##  lambdaX1  lambdaX2  lambdaX3  lambdaX4  lambdaX5 
## 0.4916233 0.0378402 0.1193047 0.4875641 0.2324578 
## 
## $optlambdaXsup
## lambdaXsup1 
##  0.03734583 
## 
## $search.method
## [1] "regparaonly"
## 
## $criterion
## [1] "BIC"
## 
## $criterion4ncomp
## [1] "BIC"
## 
## attr(,"class")
## [1] "optparasearch"

(fit342 = msma(X3, comp=opt342$optncomp, lambdaX=opt342$optlambdaX, lambdaXsup=opt342$optlambdaXsup))

## Call:
## msma.default(X = X3, comp = opt342$optncomp, lambdaX = opt342$optlambdaX, 
##     lambdaXsup = opt342$optlambdaXsup)
## 
## Numbers of non-zeros for X block: 
##        comp1 comp2 comp3 comp4
## block1     2     2     1     1
## block2     4     4     4     2
## block3     4     3     3     4
## block4     1     2     1     2
## block5     4     3     2     3
## 
## Numbers of non-zeros for X super: 
## $comp1
## comp1-1 comp1-2 comp1-3 comp1-4 comp1-5 
##       5       3       4       3       4 
## 
## $comp2
## comp2-1 comp2-2 comp2-3 comp2-4 comp2-5 
##       5       5       5       4       5 
## 
## $comp3
## comp3-1 comp3-2 comp3-3 comp3-4 comp3-5 
##       5       5       4       5       5 
## 
## $comp4
## comp4-1 comp4-2 comp4-3 comp4-4 comp4-5 
##       5       4       4       5       5

regpara1st

(opt344 = optparasearch(X3, search.method = "regpara1st", criterion="BIC", comp=c(8, 8)))

## $criterion
## [1] "BIC"
## 
## $comps
## $comps[[1]]
## [1] 1 2 3 4 5 6 7 8
## 
## $comps[[2]]
## [1] 1 2 3 4 5 6 7 8
## 
## 
## $mincriterion
## $mincriterion$comp8
## [1] -65.88419
## 
## $mincriterion$comp1
##       comp1       comp2       comp3       comp4       comp5       comp6 
##  -0.3319289  -0.7590936  -1.5204997 -69.4643431 -67.5110041 -67.3244010 
##       comp7       comp8 
## -67.2743347 -67.2668227 
## 
## 
## $criterions
## $criterions$bic
##      comp1      comp2      comp3      comp4      comp5      comp6      comp7 
##  -0.935017  -1.849394  -2.266516  -2.889812  -3.629910  -3.931964  -4.557453 
##      comp8 
## -65.884185 
## 
## $criterions$bic2
## $criterions$bic2$comp1
##       comp1       comp2       comp3       comp4       comp5       comp6 
##  -0.3319289  -0.7590936  -1.5204997 -69.4643431 -67.5110041 -67.3244010 
##       comp7       comp8 
## -67.2743347 -67.2668227 
## 
## $criterions$bic2$comp2
##       comp1       comp2       comp3       comp4       comp5       comp6 
##  -0.6233029  -1.0594666  -1.1643906 -71.2072857 -65.6936200 -65.5903894 
##       comp7       comp8 
## -65.5242167 -65.5096375 
## 
## $criterions$bic2$comp3
##      comp1      comp2      comp3      comp4      comp5      comp6      comp7 
##  -1.366974  -2.177044 -71.771974 -69.316501 -69.323144 -69.229313 -69.168400 
##      comp8 
## -69.035486 
## 
## $criterions$bic2$comp4
##       comp1       comp2       comp3       comp4       comp5       comp6 
##  -0.8947938  -1.6824952 -71.6163400 -69.8465635 -69.9392749 -70.1294851 
##       comp7       comp8 
## -70.0144783 -69.8638346 
## 
## $criterions$bic2$comp5
##       comp1       comp2       comp3       comp4       comp5       comp6 
##  -0.8635976  -1.7316051 -71.7560115 -68.2340788 -68.2091618 -68.0457070 
##       comp7       comp8 
## -67.9708741 -67.9389925 
## 
## $criterions$bic2$comp6
##     comp1     comp2     comp3     comp4     comp5     comp6     comp7     comp8 
## -71.97580 -69.41968 -69.39759 -69.37551 -69.35342 -69.33134 -69.30925 -69.28716 
## 
## $criterions$bic2$comp7
##     comp1     comp2     comp3     comp4     comp5     comp6     comp7     comp8 
## -72.88314 -69.29750 -69.27542 -69.07266 -69.23125 -69.02849 -69.00640 -68.84040 
## 
## $criterions$bic2$comp8
##     comp1     comp2     comp3     comp4     comp5     comp6     comp7     comp8 
##      -Inf -70.08564 -69.88699 -69.86491 -69.84282 -69.55459 -69.79865 -69.77656 
## 
## 
## 
## $optncomp
## [1] 8 1
## 
## $optlambdaX
##  lambdaX1  lambdaX2  lambdaX3  lambdaX4  lambdaX5 
## 0.1966493 0.5297628 0.3976823 0.4550599 0.0664165 
## 
## $optlambdaY
## NULL
## 
## $optlambdaXsup
## lambdaXsup1 
##  0.07469167 
## 
## $optlambdaYsup
## NULL
## 
## $search.method
## [1] "regpara1st"
## 
## $criterion4ncomp
## [1] "BIC"
## 
## attr(,"class")
## [1] "optparasearch"

(fit344 = msma(X3, comp=opt344$optncomp, lambdaX=opt344$optlambdaX, lambdaXsup=opt344$optlambdaXsup))

## Call:
## msma.default(X = X3, comp = opt344$optncomp, lambdaX = opt344$optlambdaX, 
##     lambdaXsup = opt344$optlambdaXsup)
## 
## Numbers of non-zeros for X block: 
##        comp1 comp2 comp3 comp4 comp5 comp6 comp7 comp8
## block1     0     4     2     0     1     0     0     2
## block2     2     2     0     1     1     0     0     0
## block3     2     0     1     1     0     1     0     0
## block4     1     2     1     2     0     0     0     0
## block5     4     4     0     0     2     0     4     0
## 
## Numbers of non-zeros for X super: 
##         comp1 comp2 comp3 comp4 comp5 comp6 comp7 comp8
## comp1-1     4     4     3     3     3     1     1     1

This is computationally expensive and takes much longer to execute.

(opt345 = optparasearch(X3, search.method = "simultaneous", criterion="BIC", comp=c(8, 8)))
(fit345 = msma(X3, comp=opt345$optncomp, lambdaX=opt345$optlambdaX))

PLS

This is computationally expensive and takes much longer to execute due to the large number of blocks.

(opt41 = optparasearch(X3, Y3, search.method = "regparaonly", criterion="BIC"))
(fit341 = msma(X3, Y3, comp=opt41$optncomp, lambdaX=opt41$optlambdaX, lambdaY=opt41$optlambdaY, lambdaXsup=opt41$optlambdaXsup, lambdaYsup=opt41$optlambdaYsup))

In this example, it works by narrowing down the parameters as follows.

(opt42 = optparasearch(X3, Y3, search.method = "regparaonly", criterion="BIC", whichselect=c("Xsup","Ysup")))

## $optncomp
## [1] 10  1
## 
## $optlambdaXsup
## lambdaXsup1 
##   0.2459844 
## 
## $optlambdaYsup
## lambdaYsup1 
##   0.2072225 
## 
## $search.method
## [1] "regparaonly"
## 
## $criterion
## [1] "BIC"
## 
## $criterion4ncomp
## [1] "BIC"
## 
## attr(,"class")
## [1] "optparasearch"

(fit342 = msma(X3, Y3, comp=opt42$optncomp, lambdaX=opt42$optlambdaX, lambdaY=opt42$optlambdaY, lambdaXsup=opt42$optlambdaXsup, lambdaYsup=opt42$optlambdaYsup))

## Call:
## msma.default(X = X3, Y = Y3, comp = opt42$optncomp, lambdaX = opt42$optlambdaX, 
##     lambdaY = opt42$optlambdaY, lambdaXsup = opt42$optlambdaXsup, 
##     lambdaYsup = opt42$optlambdaYsup)
## 
## Numbers of non-zeros for X block: 
##        comp1 comp2 comp3 comp4 comp5 comp6 comp7 comp8 comp9 comp10
## block1     4     2     0     2     0     0     0     0     3      0
## block2     0     0     4     2     0     0     0     2     3      0
## block3     0     4     0     0     3     0     0     3     0      3
## block4     0     3     0     0     4     4     3     0     0      0
## block5     4     0     0     4     3     2     0     0     0      0
## 
## Numbers of non-zeros for X super: 
##         comp1 comp2 comp3 comp4 comp5 comp6 comp7 comp8 comp9 comp10
## comp1-1     2     3     1     3     3     1     1     2     2      1
## 
## Numbers of non-zeros for Y block: 
##        comp1 comp2 comp3 comp4 comp5 comp6 comp7 comp8 comp9 comp10
## block1     0     4     0     0     0     4     1     3     3      4
## block2     0     4     4     2     0     3     0     0     2      2
## block3     2     0     2     3     3     0     0     0     0      0
## block4     4     4     3     0     0     0     2     0     2      2
## block5     2     3     0     3     0     3     0     0     0      2
## 
## Numbers of non-zeros for Y super: 
##         comp1 comp2 comp3 comp4 comp5 comp6 comp7 comp8 comp9 comp10
## comp1-1     3     4     3     3     1     3     2     1     3      4

Another example dataset is generated.

dataset4 = simdata(n = 50, rho = 0.8, Yps = rep(4, 2), Xps = rep(4, 3), seed=1)
X4 = dataset4$X; Y4 = dataset4$Y 

With this number of blocks, the calculation can be performed in a relatively short time.

(opt43 = optparasearch(X4, Y4, search.method = "regparaonly", criterion="BIC"))

## $optncomp
## [1] 10  1
## 
## $optlambdaX
##  lambdaX1  lambdaX2  lambdaX3 
## 0.1659782 0.4483126 0.1461557 
## 
## $optlambdaY
##  lambdaY1  lambdaY2 
## 0.5138342 0.4073084 
## 
## $optlambdaXsup
## lambdaXsup1 
##   0.1571942 
## 
## $optlambdaYsup
## lambdaYsup1 
##   0.1728372 
## 
## $search.method
## [1] "regparaonly"
## 
## $criterion
## [1] "BIC"
## 
## $criterion4ncomp
## [1] "BIC"
## 
## attr(,"class")
## [1] "optparasearch"

(fit343 = msma(X4, Y4, comp=opt43$optncomp, lambdaX=opt43$optlambdaX, lambdaY=opt43$optlambdaY, lambdaXsup=opt43$optlambdaXsup, lambdaYsup=opt43$optlambdaYsup))

## Call:
## msma.default(X = X4, Y = Y4, comp = opt43$optncomp, lambdaX = opt43$optlambdaX, 
##     lambdaY = opt43$optlambdaY, lambdaXsup = opt43$optlambdaXsup, 
##     lambdaYsup = opt43$optlambdaYsup)
## 
## Numbers of non-zeros for X block: 
##        comp1 comp2 comp3 comp4 comp5 comp6 comp7 comp8 comp9 comp10
## block1     2     3     3     2     0     2     2     2     2      2
## block2     2     1     1     2     0     2     3     3     2      3
## block3     3     3     3     0     3     2     2     2     2      2
## 
## Numbers of non-zeros for X super: 
##         comp1 comp2 comp3 comp4 comp5 comp6 comp7 comp8 comp9 comp10
## comp1-1     3     3     3     2     1     3     3     3     3      3
## 
## Numbers of non-zeros for Y block: 
##        comp1 comp2 comp3 comp4 comp5 comp6 comp7 comp8 comp9 comp10
## block1     2     1     1     1     0     0     0     0     0      0
## block2     3     2     2     0     2     1     1     1     1      2
## 
## Numbers of non-zeros for Y super: 
##         comp1 comp2 comp3 comp4 comp5 comp6 comp7 comp8 comp9 comp10
## comp1-1     2     2     2     1     1     1     1     1     1      1