Variable Selection in the Regularized Simultaneous Component Analysis Method for Multi-Source Data Integration
Interdisciplinary research often involves analyzing data obtained from different data sources with respect to the same subjects, objects, or experimental units. for example, global positioning systems (GPS) data have been coupled with travel diary data, resulting in a better understanding of traveling behavior. The GPS data and the travel diary data are very different in nature, and, to analyze the two types of data jointly, one often uses data integration techniques, such as the regularized simultaneous component analysis (regularized ScA) method. Regularized ScA is an extension of the (sparse) principle component analysis model to the cases where at least two data blocks are jointly analyzed, which -in order to reveal the joint and unique sources of variation -heavily relies on proper selection of the set of variables (i.e., component loadings) in the components. Regularized ScA requires a proper variable selection method to either identify the optimal values for tuning parameters or stably select variables. By means of two simulation studies with various noise and sparseness levels in simulated data, we compare six variable selection methods, which are cross-validation (CV) with the "one-standard-error" rule, repeated double CV (rdCV), BIC, Bolasso with CV, stability selection, and index of sparseness (IS)a lesser known (compared to the first five methods) but computationally efficient method. Results show that iS is the best-performing variable selection method.
<!>Discussion<!>Regularized ScA. Let<!>Table 4.<!>C<!>Repeated double cross-validation (rdCV).
<p>As a result of recent technological developments, often data from varying types of sources with respect to the same investigation units are gathered and analyzed jointly, which is referred to as multi-source data integration (also known as multi-block data analysis, linked data analysis, and in a broader sense, data fusion 1 ). In health research, joint analysis combining global positioning systems (GPS) data and self-report travel diary data for the same subjects has been shown to be insightful for understanding people's traveling behavior, purpose, and immediate environment, providing critical information relevant to health research 2 . In metabolomics, to gain a comprehensive picture of the metabolism in a biological system, researchers have conducted joint analysis on the measures obtained from two different instrumental methods, which are Mass-spectrometry (MS) with gas chromatography (GC/MS) and MS with liquid chromatography (LC/MS) [3][4][5] , on the same samples. Multi-source data integration has also been found useful in epigenetics (e.g., joint analysis on genetic information and environmental factors) 6 , in epidemiology (e.g., joint analysis on behavioral data and genetic data) 7 , and in longitudinal and life course studies (e.g., joint analysis on longitudinal survey data and bio-measures) 8 , to name a few.</p><p>A popular multi-source data integration methodology often used in social and behavior research, bioinformatics, and analytical chemistry [9][10][11][12][13][14] is the simultaneous component based data integration method (SCA for short). In essence, SCA is an extension of the well-known principal component analysis (PCA) model 15 to the cases where more than one data block is analyzed. Here, a data block can be, for example, survey data, genetic data, and behavioral data. Under certain constraints imposed on all data blocks, information shared across all data blocks can be extracted and represented by a few components. Thus, by means of dimension reduction, SCA is used to explore and interpret the internal structure that binds all data blocks together. Recent extensions of SCA have greatly improved the flexibility and the usefulness of the method by incorporating regularization such as the Lasso 16 and the Group Lasso 17 , resulting in the regularized simultaneous component analysis method (regularized SCA for short) 13,[18][19][20] . Regularized SCA reveals not only the information shared across all data blocks, Department of Methodology and Statistics, Tilburg University, Tilburg, 5000, LE, The Netherlands. *email: z.gu@ tilburguniversity.edu open simulation study extended the first one by integrating four data blocks rather than two data blocks. Both simulation studies followed the same simulation design, and therefore, in the remainder of the section, we outline the design of the first simulation study in details and mention the second simulation study when necessary.</p><p>In the first simulation study, the data were generated in five steps.</p><p>Step 1: Two data matrices, denoted by X 1 and X 2 , were generated. Here we considered three situations:</p><p>20 120 2 20 30   and</p><p>where, for all three situations, ∼ . . .</p><p>The choice of how to generate initial structures in this step has little influence on the final results as it only contributes to the true model part; other choices could also have been made, for example using an autoregressive structure on the covariance matrices. Then, the concatenated data matrix with respect to rows, denoted by</p><p>2 , was of dimension 20 × 50, 20 × 150, and 80 × 50, respectively. In the following, we use the first situation (i.e., Eq. 1) as an example to explain the remaining steps.</p><p>Step 2: Using singular value decomposition (SVD), we decomposed ∼ X C into UΣV. We defined the "true" component score matrix, denoted by T true , as the matrix containing the three left singular vectors in U corresponding to the three largest singular values. Let Σ ∼ denote the diagonal matrix containing the three largest singular values, and let ∼ V denote the matrix containing the three right singular vectors corresponding to the three largest singular values. Then, the non-sparse component loading matrix, denoted by  . We assumed that the first component of</p><p>Step 4: We replaced some loadings in p</p><p>. Note that for notational convenience we used the same symbols for the sparsified loading vectors as previously.</p><p>Step 5: We computed = X T P ( )</p><p>true T , and added a noise matrix, denoted by E, to X C true to generate the final simulated dataset, denoted by X C generated , so that</p><p>, where the scalar α is a scaling factor. The cells in E were generated from N(0, 1). Note that an implicit assumption of PCA and also SCA is independent and identically distributed noise; other types of noise structure may affect the results. By adjusting α, we were able to control the proportion of noise variance in X C generated . We considered two noise levels: 0.5% and 30% of variance in X C generated were attributable to noise. In summary, the first simulation study included the following design factors:</p><p>• Three situations of X 1 and X 2 (i.e., Eqs. 1, 2 and 3).</p><p>• Two sparseness levels in p 1 1 , p 1 2 , p 2 2 , and p 3 1 : 30% and 50%. • Two noise levels: 0.5%, and 30%.</p><p>The design factors were fully crossed, resulting in × × = 3 2 2 12 design cells. In each design cell, we simulated 20 datasets following the above five steps, and therefore in total 240 datasets were simulated. Then, for each dataset, we conducted the regularized SCA analysis and compared the results generated by the model selection methods, which are CV with "one-standard-error" rule, rdCV, BIC, IS, Bolasso with CV, and stability selection. The design of the second simulation study also involved five steps similar to the first simulation, but we made the following changes. In Step 1 of the second simulation study, we considered only one situation: In summary, the second simulation study included the following two design factors:</p><p>, and p 1 4 : 30% and 50%. • Two noise levels: 0.5%, and 30%. www.nature.com/scientificreports www.nature.com/scientificreports/ The design factors were fully crossed, resulting in × = 2 2 4 design cells. In each design cell, we simulated 20 datasets following the above five steps, and therefore in total 80 datasets were simulated.</p><p>Performance measures. To compare the variable selection methods, we used two types of performance measures. The first type concerned the component loading matrix, and the second type concerned the component score matrix. The first type consisted of three performance measures. Let PC denote the estimated concatenated component loading matrix. The first performance measure, denoted by PL, was the proportion of non-zero and zero loadings correctly identified in PC compared to P C true :</p><p>P number of correctly selected non zero loadings number of correctly identified zero loadings total number of loadings in C true Notice that ∈ PL [0, 1]. Intuitively, for regularized SCA, the best model selection method should be the one that generating the highest PL among the methods. In addition to PL, we also used PL non-0 loadings , defined as  18 . Note that we are interested in the associations among items within a component, and the associations are indicated by the signs of the loadings. Take Component 2 for example. The three non-zero loadings have the same sign (in this case "−" sign), meaning that mother's "activities with children", father's "activities with children", and father's "communication with children" are positively associated with each other. Two loadings having opposite signs indicates a negative association between the two items. We remind the reader that, when interpreting the loadings and the associations among them, one should also take into account how the items are scored (see Table 1). For example, a higher score on "relationship with parters" indicates a more satisfied relationship. A higher score on "argue with partners" indicates a less violent relationship. www.nature.com/scientificreports www.nature.com/scientificreports/ We used PL non-0 loadings to evaluate how well a model selection method assisted correctly retaining non-zero loadings and used PL 0 loadings to evaluate how well a model selection method assisted correctly identifying zero loadings.</p><p>In this study, we focused on the component loading matrix, and we used the variable selection methods to help us identify non-zero and zero loadings, but the component score matrix was also important. Ideally, we would prefer an estimated component score matrix as close as possible to the true component score matrix. Therefore, the second type of performance measure evaluated the degree of similarity between T true and the estimated component score matrix T, quantified by Tucker congruence ϕ 42</p><p>Ideally, a good model selection method for regularized SCA is the one that makes ϕ close to 1.</p><p>Results. We used the R package RegularizedSCA (version 0.5.5) 20 to estimate the regularized SCA model; the R script for replicating the study is included in the supplementary material. All columns in the simulated datasets were mean-centered and scaled to norm one. We used the Group Lasso penalty to identify component structure (i.e., common/distinctive components) and used the Lasso penalty to impose sparseness within a component. For details, please see the Methods section.</p><p>Figures 2, 3, 4 and 5 summarize the results of the first simulation, where two data blocks were integrated. Specifically, Figs. 2, 4 and 5, by means of boxplots, present the performance measures PL (Eq. 6), PL non-0 loadings (Eq. 7), and PL 0 loadings (Eq. 8), respectively. Figure 3 presents the boxplots of Tucker congruence measures. For each figure, the upper, middel, and bottom panels correspond to the first, second, and third situations of X 1 and X 2 (i.e., Eqs. 1, 2 and 3), respectively. The reader may notice that most methods (except for BIC and Bolasso) did not differ much in Tucker congruence, and therefore, we focus on discussing PL, PL non-0 loadings , and PL 0 loadings and mention Tucker conguence only when necessary.</p><p>Based on the figures, we concluded the following. First, CV with "one-standard-error" rule and rdCV did not outperform the other methods in most cases in terms of correctly identifying non-zero and zero loadings (see Fig. 2). Figures 4 and 5 show that the two methods tended to retain more non-zero loadings than needed, resulting in high PL non-0 loadings but low PL 0 loadings , which is a known feature of CV-based methods 43 . Second, stability selection was the best-performing method in terms of PL. However, as we have explained in the Methods section, in order for the method to work in the simulation, we assumed that the correct number of non-zero loadings was known a priori, which is unrealistic in practice. Third, IS was the second best-performing method (Fig. 2), witnessed by a balanced, high PL non-0 loadings (Fig. 4) and high PL 0 loadings (Fig. 5). Fourth, BIC performed worse than the other methods (except for Bolasso) when the noise level was high (i.e., 30%). Figures 4 and 5 suggest that BIC consistently favored very sparse results, resulting in very high PL 0 loadings but low PL non-0 loadings , which in turn lead to low Tucker congruence values (Fig. 3). Finally, Bolasso performed the worst among all the methods in terms of PL and Tucker congruence. This is primarily because the algorithm is very strict: A loading was identified as a non-zero loading only if the loading was estimated to be different from zero in all 50 repetitions (see the Methods section). As a result, the algorithm generated an estimated loading matrix with too many zeros -that is, very high PL 0 loadings in Fig. 5 and very low PL non-0 loadings in Fig. 4. Figures 6, 7, 8 and 9 present the results of the second simulation study, where four data blocks were integrated. It may be noted that the four figures are very similar to the Figs. 2, 3, 4 and 5, and therefore, similar conclusions can be made for the second simulation study. For the sake of simplicity, we do not discuss the Figs. 6, 7, 8, and 9.</p><p>Based on the two simulation studies, we conclude that, in practice, IS is the best-performing variable selection method for regularized SCA. In addition, more research is needed to improve the stability selection algorithm for regularized SCA so that it will no longer rely on the unrealistic assumption that the correct number of total non-zero loading is known a priori. empirical examples. In this section, we present three empirical applications of regularized SCA combined with IS for variable selection. We used the first two empirical examples to explain to the reader how to interpret the estimated component loading matrix generated by regularized SCA together with IS in applied research. The third empirical example is the parent-child relationship data discussed in the Introduction section. We reanalyzed the data by using IS and compared the results with Table 2. We remind the reader that, to evaluate and to interpret the results generated by regularzed SCA, one typically resorts to both the estimated component loading matrix and the estimated component score matrix. In this article, because we focus on variable selection in the component loading matrix, we refrain from discussing the interpretation of the estimated component score matrix in the remainder of this section. Furthemore, for detailed explanation on the use of regularized SCA and the interpretation of the results, we refer to Gu and Van Deun 18 .</p><p>We used the following setup for IS: 50 Lasso tuning parameter values (equally spaced ranging from 0.0000001 to the smallest value making the entire estimated component loading matrix a zero matrix), and 50 Group Lasso tuning parameter values (equally spaced ranging from 0.0000001 to the smallest value making the entire estimated component loading matrix a zero matrix). All columns in the empirical datasets were mean-centered and scaled to norm one before the regularized SCA analysis was performed.</p><p>Joint analysis of the Herring data. In food science, researchers are often interested in the chemical/physical characteristics and the sensory characteristics of a certain food item and analyze the characteristics jointly. An example is the Herring data obtained from a ripening experiment 44,45 . In this article, we used part of the original Herring data 20 , consisting of two datablocks. The first block contained the physical and chemical changes, including pHB, ProteinM, ProteinB, Water, AshM, Fat, TCAIndexM, TCAIndexB, TCAM, and TCAB, of 21 salted herring samples. The meaning of the labels of the physical and chemical changes can be found at http://www. models.life.ku.dk/Ripening_of_Herring. The second block contained the sensory data, including features such as ripened, rawness, malt, stockfish smell, sweetness, salty, spice, softness, toughness, and watery, of the same 21 samples. An interesting research question is whether certain physical and chemical changes are associated with certain sensory characteristics of the herrings. It may be noted that, in this article, we do not discuss how to identify the number of components R (see the Methods section), and for this topic, we refer to Gu and Van Deun 18 . A previous study 18 suggested that, for the Herring data, the reasonable number of components R was 4. Therefore, we performed the regularized SCA analysis with IS and = R 4, and the estimated component loading matrix is presented in Table 3. The table suggests that, for each component, not all variables were important. For example, for Component 1, variables pHB, Water, and AshM from the block of "physical and chemical changes" and variables Ripened, Rawness, Stockfish smell, Sweetness, and Spice from the "sensory" block were important and therefore their loadings were different from zero. To interpret the associations among the variables of Component 1, we primarily look at the signs of the non-zero loadings. For example, for Component 1, variables pHB, Water, Rawness, Sweetness, and Spice were negatively associated with variables AshM, Ripened, Stockfish smell. The remaining three components can be interpreted in the same way. Joint analysis of metabolomics data. In metabolomics, researchers often use multiple instrumental methods to measure as many metabolites as possible and perform joint analyses by combining the measures on the same metabolites gathered from different instrumental methods 5 . The dataset used in this article contained measures of 28 samples of Escherichia coli (E. coli) obtained from using two measurement methods, which were mass spectrometry with gas chromatograph (GC/MS) and mass spectrometry with liquid chromatography (LC/MS) 3,4 . The dataset contained a block of GC/MS data with 144 metabolites and a block of LC/MS data with 44 metabolites. For a detailed description of the dataset, including the experimental design and conditions for obtaining the measures, we refer to Smilde, Van der Werf, Bijlsma, Van der Werff-van der Vat, and Jellema 5 . A previous study 19 suggested that the appropriate number of components R was five. We thus performed the regularized SCA analysis with IS and = R 5. It may be noted that, in this example, because of the large number of variables, a table of estimated component loading matrix such like Table 3 usually is not practical. Instead, researchers typically use a heatmap so as to get some impression about the sparseness of the loading matrix. Figure 10 presents such a heatmap for the estimated component loading matrix. We found that many loadings in Fig. 10 were very close or equal to zero. As a side note, for this study, researchers typically focus on interpreting the estimated component score matrix instead of the estimated component loading matrix (see, e.g., Van Deun, Wilderjans, van den Berg, Antoniadis, and Van Mechelen 46 ).</p><p>Re-analysis of the parent-child relationship survey data. Table 4 presents the estimated component loading matrix obtained by using IS. The orders of the components were adjusted by using Tucker congruence so that the components in Table 4 are comparable to the components in Table 2 which were generated by using CV 18 . The two estimated component loading matrices in Tables 4 and 2 are comparable, and the conclusions based on the two tables are almost the same. For example, for Component 1 of both tables, the last 7 variables from the "Mother"</p><!><p>In this article, we examined six variable selection methods suitable for regularized SCA. The popular CV-based variable selection methods, including CV with "one-standard-error" rule and rdCV, did not outperform other methods. This result may be surprising to many researchers, especially considering that CV seems to be the standard practice when it comes to variable selection. The poor recovery rate of component loadings by using the CV-based methods in the simulations showed that the CV-based methods retained more loadings than needed. Stability selection is a promising method, but at this moment we do not know how to identify an accurate lower bound for the expected non-zero loadings (i.e., Q), making it impossible to tune λ L . Thus, we advocate the use of IS. It is possible that a hybrid method combining IS and stability selection may perform better than IS. For  www.nature.com/scientificreports www.nature.com/scientificreports/ example, one first uses IS to decide the total number of non-zero loadings and then uses stability selection given the total number of non-zero loadings. Further examination on this idea is needed.</p><p>We focused on determining the status of the components (i.e., common/distinctive structure) and their level of sparseness. Another important issue that remains to be fully understood is the selection of the number of components R. Because the goal of this article is to understand variable selection methods for the component loading matrices, the selection of R is beyond the scope of this article. For interested readers, we refer to Bro, Kjeldahl, Smilde, and Kiers 47 , Gu and Van Deun 18 , and Måge, Smilde, and van der Kloet 48 . We believe that more studies are needed to evaluate the performance of model selection methods for determining R and the performance of variable selection. This may be done sequentially (i.e., first determining R and then, given R, performing variable selection) but also simultaneously (for example, using the index of sparseness to determine R and to perform variable selection at the same time). Finally, we call for studies on comparing the performance of variable selection methods in regularized models. The six variable selection methods studied in this article originated in sparse PCA literature. Therefore, we suspect that stability selection and IS would still outperform the other five methods in the sparse PCA settings. However, we are not aware of any study that compares variable selection methods in sparse PCA.</p><p>Admittedly, the six methods studied in this article do not constitute an exhaustive list of all possible variable selection methods for regularized SCA. Other variable selection methods exist, such as the method by Qi, Luo, and Zhao 49 , the information criterion by Chen and Chen 43 , and the numerical convex hull based method 50 , but they cannot be readily adapted to be used together with regularized SCA. These methods are promising though, and therefore require full attention in separate articles.</p><!><p>denote the kth data block with I rows representing subjects, objects, or experimental conditions measured on J k variables. One may notice that I does not have a subscript k, meaning that all K data blocks are to be analyzed jointly with respect to the same I subjects, objects, or experimental conditions. Each data block may have a different set of variables. Let</p><p>k k denote the concatenated data matrix, which is obtained by concatenating X k s with respect to rows (i.e.,</p><p>). Note that I may be much smaller than J k (i.e., high-dimensional data). Let</p><p>is the Lasso penalty 16 , and its corresponding tuning parameter is</p><p>is the Group Lasso penalty 17 , and its corresponding tuning parameter is λ G . Note that if λ = 0 L and λ = 0 G , Eq. 10 reduces to a least squares minimization problem. As a side note, before performing the regularized SCA analysis, all columns in X k may be mean-centered and scaled to norm one or to − J k 1/2 in order to give all blocks -even those that contain relatively few variables -equal weight; This procedure is referred to as data pre-processing. However, one may notice that in Eq. 10 the Group Lasso penalty is also weighted by J k . Thus, it is likely that, when data are scaled to − J k 1/2 , Eq. 10 would favor data blocks with fewer variables, because the Group Lasso penalty takes J k into account. In addition, because in this study we are interested in identifying the associations between (some) variables across data blocks, penalties are imposed on the component loading matrix 19,46 . T is assumed to be same for all K data blocks, and therefore it serves as a "bridge" linking all data blocks. Information shared among all data blocks or unique to some blocks, such as the loadings in</p><p>Let T denote the estimated component score matrix based on Eq. 10, and let Pk denote the estimated component loading matrix for the kth data block. Further, Let</p><p>denote the concatenated estimated component loading matrix, which is obtained by concatenating all P s k with respect to the columns (i.e.,</p><p>). The algorithm for estimating Eq. 10 requires an alternating procedure where T and PC are estimated iteratively. Given PC , T is obtained by computing = T VU T , where UΣV T is the SVD of P X</p><p>in Eq. 11 18 :</p><p>In Eq. 12, ⋅ S( ) denotes the soft-thresholding operator. The operator [x] + is defined as</p><p>, if > x 0, and</p><p>For details of the estimation procedure, see Algorithm 1 of Gu and Van Deun 18 . Information regarding the position of non-zero/zero loadings in P C may be known a priori. For example, Bolasso and stability selection procedures, which will be discussed shortly, can be used to identify the position of non-zero/zero loadings. Once the position of non-zero/zero loadings is identified, one uses regularized SCA with</p><p>to re-estimate the non-zero loadings in P C while keeping the zero loadings fixed throughout the estimation procedure. For details of the estimation procedure, see Algorithm 2 of Gu and Van Deun 18 .</p><p>Variable selection methods. The variable selection methods discussed in this article can be categorized into two groups. The first group, including CV with "one-standard-error" rule, rdCV, BIC criterion, and IS, aims at identifying the optimal λ L and λ G for Eq. 10. Once the optimal λ L and λ G are obtained, one re-estimates the model by using the optimal λ L and λ G . The second group, including the Bolasso with CV and stability selection,</p><!><p>The parent-child relationship data: Estimated component loading matrix generated by using regularized SCA with IS. Please be noted that the signs of components 1, 2, and 5 were manually changed from positive to negative. The signs of Component 3 were manually changed from negative to positive. Due to the invariance of signs of regularized SCA, changing signs do not influence the interpretation of loadings. Therefore, we changed the signs to make it easier for the reader to compare the www.nature.com/scientificreports www.nature.com/scientificreports/ aims at identifying the position of non-zero/zero loadings in P C through repeated sampling. Once the position of non-zero/zero loadings is identified, one re-estimates the non-zero loadings while keeping the zero loadings fixed at zero. In the remainder of this article, we assume that the number of components R is known. To identify R in practice, one may use the Variance Accounted For (VAF) method 9,10 and the PCA-GCA method 14 . Both methods are included in the R package "RegularizedSCA" 20 (for details on how to use the two methods, see Gu and Van Deun 18 ). We remind the reader that more research is needed for fully understanding how to identify R.</p><p>CV with "one-standard-error" rule. Given a set of λ L s (consisting of evenly spaced increasing values ranging from a value close to zero, say, 0.000001, to the smallest value making = P 0</p><!><p>), denoted by Λ L , and a set of λ G s (also consisting of evenly spaced increasing values ranging from a value close to zero to the smallest value making = P 0</p><p>), denoted by Λ G , the algorithm searches through a grid of λ L s and λ s s (i.e., the Cartesian product of Λ L and Λ G ). For each combination of λ L and λ G , denoted by λ λ ( , ) L G , the algorithm conducts K-fold CV. Take 10-fold CV for example, 10% of the data cells in X C are replaced with missing values, and afterwards, missing values in each column are replaced with the mean of that column. The algorithm then computes the mean squared prediction errors (MSPE) 51 for each λ λ ( , )</p><p>) is performed. Let X k q ( ) denote the data from the kth block for the qth fold. Let Pk q ( ) denote the estimated component loading matrix for the kth data block for the qth fold. Let T q ( ) denote the estimated component score matrix for the qth fold. Then MSPE is</p><p>). Let L G . Applying the "one-standard-error" rule 26 , the algorithm searches for the optimal pair, denoted by λ λ ( , )</p><p>As a side note, in the simulation, the algorithm searched the optimal pair whose MSPE was closest to (i.e., could be slightly larger or smaller than)</p><p>L G L G</p><p>. In the simulation, we used 5-fold CV.</p><!><p>The rdCV 27 , as its name would suggest, is an algorithm that performs double CV repeatedly. Double CV consists of two so-called "layers", and at each layer a CV is executed. Figure 11 presents a sketch of the algorithm. In the ρth repetition (ρ = … P 1, , repetition ), the concatenated dataset, X C , is randomly split into T segments with a (nearly) equal sample size; that is, each segment contains (roughly) the same number of subjects/objects/experimental conditions. The τth segment, denoted by are obtained from the traditional simultaneous component model without Lasso and Group Lasso penalties. We define the of the residual matrix given λ L and λ G , denoted by ∼ V , as</p><p>, where T and PC are obtained from Eq. 10. We define the degrees of freedom given λ L and λ G , denoted by λ λ df ( , ) L G , as the number of non-zero loadings in PC . Then the BIC criterion adjusted for regularized SCA, given λ L and λ G , based on Croux et al. is</p><p>and the BIC criterion adjusted for the regularized SCA method based on Guo et al. is</p><p>Notice that the BIC in Eq. 14 is exactly I times the BIC in Eq. 13. Thus, the two methods are in fact equivalent. Then, the optimal tuning parameter values, λ λ ( , ) ). Let # o denote the total number of zero loadings in PC . Then IS, according to Gajjar, Kulahci, and Palazoglu 28 and Trendafilov 29 , is Bolasso with CV. Bolasso, originally proposed by Bach 31 , has been extended to a hybrid procedure combining the original Bolasso with CV 32,33 for stably selecting variables in Lasso regression. Figure 12 presents the algorithm of the Bolasso with CV. In essence, the Bolasso is a bootstrapping procedure. For each bootstrap sample, regularized SCA with K-fold CV is executed, generating the optimal tuning parameters, λ λ ( , ) o . Let P repetition denote the total number of repetitions. Then in total P repetition P s C are generated. The algorithm then compares the P repetition P s C , checks which loadings have been estimated to be not zeros for P repetition times, and records the corresponding index set. As a result, an index set containing the position of non-zero loadings is obtained. Finally, PC and T are estimated given the index set. One may notice that because of the invariance of the regularized SCA solution under permutations of components 18 , the P s C must first be adjusted according to a reference matrix by using the Tucker congruence 42 (for details, see the R script provided in the supplementary material). As a side note, in the simulation, we used 5-fold CV and let = P 50 repetition .</p><p>Stability selection. Stability selection 25 was demonstrated for variable selection in regression analysis and graphical models based on the Lasso. To use this method for regularized SCA, we have made a few adjustments and present the algorithm in Fig. , the algorithm works as follows. First, 100 samples with ⌊ ⌋ I/2 subjects (i.e., rows) from X C are randomly drawn without replacement. For each sample created, regularized SCA with λ L s ( ) and λ = 0 G is applied. Therefore, the algorithm generates 100 P s C . Because of the invariance of regularized SCA solution under permutations of components, the P s C are adjusted according to a common reference matrix by using the Tucker congruence (for details, see the R script in the supplementary material). Then, the algorithm counts the number of times that the same loading is estimated to be a non-zero loading across the 100 P s C , which is then divided by 100, resulting in the selection probability for that loading (see Step 1(d) in Fig. 13). As a result, each component loading has a selection probability, which is then compared to a pre-defined selection probability threshold π thr , and the loadings for which the selection probabilities lower than π thr are constrained to be zero loadings. The error control theorem proposed by Meinshausen and Bühlmann 25 (Theorem 1, p. 7) adjusted for the regularized SCA model is  where EV denotes the expected number of falsely selected variables, Q denotes the expected non-zero loadings, and ∑ R J k k is the total number of loadings. We notice that, when Gu and Van Deun 19 applied stability selection in their study on regularized SCA, they failed to recognize the problem of Eq. 16: When used for regularized SCA, the lower bound for Q produced by Eq. 16 is not strict enough, making it difficult to tune Λ L . To explain, we use the first simulation study in the Results section as an example and consider the situation of , then ≥ Q 19, which is much smaller than 150 (i.e, the total number of non-zero loadings). Thus, using Eq. 16 to tune Λ L is likely to generate a component loading matrix that is too sparse. In this article, the algorithm tunes Λ L by using the number of expected non-zero component loadings Q, which is assumed known a priori (see Step 1(e) in Fig. 13). Thus, given λ L s ( ) , if the total number of loadings with selection probability not lower than π thr is equal to or larger than Q, then the algorithm ignores the remaining Lasso tuning parameter values () . The algorithm records the maximum selection probability across the s selection probabilities for each loading, ranks the loadings in descending order according to their associated maximum selection probabilities, and picks the loadings whose maximum probabilities belong to the first Q maximum probabilities (see steps 2, 3, and 4 in Fig. 13). Finally, the selected loadings are re-estimated, while the remaining loadings are fixed at zero. As a side note, in the simulation, we set π = . 0 6 thr . Also in the simulation, Q was known, which was the total number of non-zero loadings in P C true , but this is unrealistic in practice.</p>