Reproducible analysis of disease space via principal components using the novel R package syndRomics

Biomedical data are usually analyzed at the univariate level, focused on a single primary outcome measure to provide insight into systems biology, complex disease states, and precision medicine opportunities. More broadly, these complex biological and disease states can be detected as common factors emerging from the relationships among measured variables using multivariate approaches. ‘Syndromics’ refers to an analytical framework for measuring disease states using principal component analysis and related multivariate statistics as primary tools for extracting underlying disease patterns. A key part of the syndromic workflow is the interpretation, the visualization, and the study of robustness of the main components that characterize the disease space. We present a new software package, syndRomics, an open-source R package with utility for component visualization, interpretation, and stability for syndromic analysis. We document the implementation of syndRomics and illustrate the use of the package in case studies of neurological trauma data.


Introduction
The goal of the burgeoning field of precision medicine is to understand complex disease states and provide opportunities for deep patient phenotyping and highly targeted therapeutics. Precision medicine requires an understanding of multidimensional disease states. Yet, the analysis of biomedical data remains largely univariate, with response variables considered individually and reports involving several distinct analyses. This analytical approach limits our interpretation of the complexity of a disease by not considering the shared information across variables and potentially contributing to irreproducibility due to statistical limitations of multiple comparison testing. Understanding the full set of interrelated disease features through multivariate statistics is the goal of the growing domain of 'syndromics' (Ferguson et al., 2011). In particular, principal component analysis (PCA) and related multivariate statistics such as nonlinear PCA or factor analysis have been proposed as tools for extracting underlying factors or patterns (principal components [PCs]) reflecting disease states Haefeli et al., 2017a;Haefeli et al., 2017b;Kutcher et al., 2013;Nielson et al., 2014;Nielson et al., 2015;Panaretos et al., 2017;Rosenzweig et al., 2010;Rosenzweig et al., 2018;Rosenzweig et al., 2019;Zhang and Castelló , 2017). There are several other multivariate methods that could be used for multivariate pattern detection: other ordination and dimension reduction techniques, cluster analysis, discrimination analysis, or the plethora of more recent machine learning methods. The use of any of these methods has its advantages and pitfalls (Everitt and Hothorn, 2011). We focus on PCA as being one of the most widely used method for pattern detection. PCA is a multivariate statistical procedure that allows for the generation of new uncorrelated variables, called PCs, as a weighted combination of the original variables (Abdi and Williams, 2010;Hotelling, 1933;Jolliffe and Cadima, 2016). These components are ordered such that the first component explains the major source of variance in the data, the second component the second largest source of variance, etc. The extracted components reflect the interrelation between all the original variables or features, allowing for disease pattern detection, guiding in the interpretation of disease complex space and overcoming univariate analysis limitations.
Despite the extensive use of PCA in some subfields of biological research and the increasing use of PCA for disease pattern discovery, there is very limited information in the literature that can guide applied biomedical researchers about its implementation and interpretation. Here, we offer a practical guide to the application of PCA for the extraction of disease patterns that conform the disease space, with focus on reproducibility. By no means can we cover the extensive field of PCA in the present document. Rather, we aim to provide an introductory manual to extraction of reproducible disease patterns using multidimensional analytics, directed to biomedical researcher practitioners while pointing to additional relevant sources of information. We introduce a software package for the R programming language called syndRomics, implementing some of the tools described here. We will illustrate the analysis workflow and the use of the package in experimental data from case studies in neurotrauma.
The key steps in disease pattern detection by PCA are shown in Figure 1. The syndRomics package offers functionalities that aid in these steps, building on the extensive PCA framework developed by the R open-source community. The package implements a novel visualization tool, the syndromic plot first published by Ferguson et al., 2013, as well as functions to quickly generate two other publication-ready visualizations (a heatmap and a barmap). In addition, the package implements resampling strategies, providing data-driven approaches to analytical decision-making aimed to reduce researcher subjectivity and increase reproducibility. In particular, the package offers a function to extract metrics for component and variable significance by using nonparametric permutation methods (Landgrebe et al., 2002;Linting et al., 2011;Peres-Neto et al., 2003), to inform component selection and component interpretation. Finally, the package incorporates functions to study component stability toward understanding the generalizability and robustness of the analysis (Cattell and Baggaley, 1960;Cattell et al., 1969;Lorenzo-Seva and ten Berge, 2006).

Results
We will describe the general steps to use PCA for syndromic analysis and illustrate the use of the syndRomic package along the analytical steps with two case studies of neurotrauma data. Details of the usage and implementation of the package and functions are described in the Materials and methods section. The full code reproducing the analysis can be found in the supplementary material. Code boxes in the text provide snippets illustrating the main sections of the code. The first case study is used as a tutorial to illustrate the steps of analysis; the second case study is discussed at the end of the results section. For the first case, we used a publicly available preclinical dataset on the Open Data Commons for Spinal Cord Injury (odc-sci.org) (Callahan et al., 2017;Fouad et al., 2019). We selected a subset of the dataset with accession number ODC-SCI: 26  that has been previously used for deriving the so-called spinal cord injury (SCI) syndromics . The dataset contains 159 subjects (rats) that have been studied on different motor functional outcomes across time after cervical spinal cord injury. The subset chosen for the present analysis consists of 18 outcome variables measured at 6 weeks after injury. The included variables for this analysis are shown in Table 1. For additional details of these variables, see Ferguson et al., 2013. Step 1: Extracting PCA solution from the data There is extensive literature on performing PCA (Abdi and Williams, 2010;Jolliffe and Cadima, 2016;Zhang and Castelló , 2017). As a consideration, biomedical data aiming to capture the multivariate disease space usually contains variables of different types (i.e. categorical, continuous, etc.) and scales, known as 'mixed-type' data. Moreover, missing data is a common problem in biomedicine (Hollestein and Carpenter, 2017;Kaushal, 2014;Nielson et al., 2020) that needs to be solved Table 1. List of variables included in the first case study.

Variable Definition
wtChng Change of animal weight (grams) from day of Injury to 6 weeks post-injury to be able to apply most standard PCA algorithms. Therefore, some pre-processing transformations are usually applied before performing PCA. For example, linear PCA is sensitive to the scale of variables, thus when applying a linear PCA to continuous variables of different units or scales, a common practice is to scale the data to unit variance first (i.e. equivalent to performing the PCA on the correlation matrix). The use of the package to conduct syndromics analysis from linear PCA is illustrated on the first case study. In cases of datasets with mixed data types and/or non-linear relationships between variables, nonlinear PCA with optimal scaling transformation (Linting et al., 2007a;Mair and Leeuw, 2019) has been previously used for disease pattern analysis (Rosenzweig et al., 2018;Rosenzweig et al., 2019). We used the syndRomics package to analyze patterns from a nonlinear PCA in the second case study. In cases with missing data, strategies such as data imputation or the use of PCA algorithms allowing missing values might be needed (Dray and Josse, 2015).
While missing values analysis and dealing with missingness is an extensive topic that is not covered in detail here (Rubin, 1976;Buuren, 2018), the chosen case studies do contain missing values and illustrate how the package can help to determine the stability of the PCs when imputing missing values (see component stability section). Another consideration is selecting which variables to include in the analysis. For PCA of experimental data where there are stratifying factors (e.g. control vs. treatment), it is important to leave out variables that directly capture the variance of these factors, which would bias PCA results toward separating the experimental groups. This bias is problematic since in syndromic analysis, the goal is to find the relationship between variables describing different diseases states in an unsupervised (i.e. not guided by our design) manner. For instance, if treatment indicators are included and the variance between treatment groups is high, the PCA solution would directly capture the experimental design and confound the multivariate patterns.
The disease components can be used in subsequent analysis as multivariate outcomes or predictor indicators (Haefeli et al., 2017a;Nielson et al., 2015;Rosenzweig et al., 2018;Rosenzweig et al., 2019). PCA is used to extract the correlation structure between variables, generating new independent variables as linear combinations. Beyond the use of PCs as proxies for disease patterns, the PCs can help mitigate issues that might appear when analyzing several variables such as multicollinearity, overfitting, and multiple testing (Altman and Krzywinski, 2018;Johnson et al., 1973;Lever et al., 2017).
The reader is referred to some materials of interest on considerations and limitations when conducting PCA and related methods for biomedical research (Jiang and Eskridge, 2000;Konishi, 2015;Nguyen and Holmes, 2019;Zhang and Castelló , 2017).
Case study: In the first case study, the goal is to run a linear PCA to study the motor function components 6 weeks after cervical spinal cord injury. This will summarize all motor function variables as a small set of independent components explaining different aspects of the motor behavior after an SCI. The data contains missing values (Figure 4-figure supplement 1), and therefore we performed missing values analysis before continuing with the workflow. Typically, the first step in missing values analysis is to determine patterns of missingness and classify missing values as missing completely at random (MCAR), missing at random (MAR) or missing not at random (MNAR) (Rubin, 1976). The type of missingness will guide the decision on which is an acceptable procedure to deal with missing values. For instance, deleting all subjects that contain at least one missing observation is common practice (aka listwise deletion or complete-case analysis), but it is only acceptable if missing values are MCAR. Otherwise, the robustness and proper estimation of the missing values can not be guaranteed (Schafer and Graham, 2002;Buuren, 2018). In the example data, subjects have been pooled together from different experiments. We know that the observed pattern of missingness (Figure 4-figure supplement 1) is due to a set of animals where some of the outcome measures were not studied, suggesting that missing values are MNAR. We confirmed that missing values are not MCAR using a previously described test of MCAR (Jamshidian and Jalal, 2010) implemented in the MissMech package in R (Jamshidian et al., 2014), which rejected the hypothesis of MCAR missingness in our data. Thus, excluding subjects from the analysis is not justified. Instead, we have used multiple imputation through the mice R package (Buuren and Groothuis-Oudshoorn, 2011) to generate 50 imputed datasets and pooled them using the mean of each observation. We will illustrate on the component stability section how the syndRomics package can be used to determine the robustness of multiple imputation for disease pattern analysis. We extracted the PCA solution of the pooled imputed data using the prcomp() function in R after centering and scaling the data to unit variance (R code box 1). Other similar functions in R or other software can be used.
Step 2: Component selection: how many components to keep?
The first question, after running PCA for extracting the disease components is usually to determine how many PCs are relevant. As a general consideration, the PCs with lower eigenvalues (i.e. explain less variance) have a higher chance of representing noise in the data (Jolliffe and Cadima, 2016), questioning their generality and value. The goal is to determine the minimal set of components that can be used to describe the disease space. Importantly, there is not a single, specific rule for this determination. A common method in PCA and related methods is the Scree test by Cattell, 1966, where all PCs are ordered in descending rank by their eigenvalues, and PCs above the 'elbow' are retained. Another criterion is the eigenvalue greater than one rule which is applied to standardized PCAs (from the correlation matrix) with the criteria of only keeping PCs with an eigenvalue (i.e. the variance of a component) above 1 (Guttman, 1954;Kaiser, 1960). A more thorough description of these and others methods can be found elsewhere (Glorfeld, 1995;Horn, 1965;Vitale et al., 2017;Zwick and Velicer, 1986). Simulations have shown these methods (specially the eigenvalue greater than one rule) to be less robust than a re-sampling approach for selecting the number of relevant components (Zwick and Velicer, 1986). The syndRomics package incorporates a nonparametric permutation test approximated through Monte Carlo re-sampling of the total 'variance accounted for' (VAF) of each PC to aid in the selection of relevant PCs (Buja and Eyuboglu, 1992;Glorfeld, 1995;Horn, 1965;Landgrebe et al., 2002). The permutation test can also assist in component interpretation by studying the contribution of each variable to the PCA solution (Buja and Eyuboglu, 1992;Linting et al., 2011) as we will see in the next section.
The goal of the permutation test is to determine whether the extracted PCs can be considered to be generated not-at-random. This method has been shown to outperform parametric tests for PCA in situations similar to biomedical data where sample sizes are relatively small and the data rarely comply with the assumptions of the models (Buja and Eyuboglu, 1992;Horn, 1965;Zwick and Velicer, 1986). In that regard, a hypothesis test is defined as: H (null) :PC VAF is indistinguishable from a random generation H (alternative) :PC VAF is different from random The p values are calculated by: where q is the number of times the chosen metric is higher in the permuted distribution than in the original PCA solution and P is the number of permutations (Buja and Eyuboglu, 1992). Rejecting the null hypothesis is interpreted as evidence of the tested PC being generated from true signal and not by random noise. This sets a lower bound for which PCs to consider 'important' above noise, but does not indicate the magnitude of the 'importance', which is represented by VAF. Importantly, for datasets with several directions of variance and high signal-to-noise ratio, PCs with low VAF can still be statistically significant. The value of interpreting such PCs must be judged by the researcher in the context analysis in question. It is also important to consider how big P needs to be when performing re-sampling, such as with the permutation test incorporated in the package. The reader should note that the lowest p value that can be calculated is dependent on P. For example, if P is set to a value of 10 (a relatively low value), the smallest p value that can be detected is 0.09, which occurs when q ¼ 0. Accordingly, P should be set high enough to reach the desired minimum p value. Moreover, simulation studies have shown that P under 99 have low power and a minimum of 499 permutations is recommended (Buja and Eyuboglu, 1992;Abdi and Williams, 2010;Linting, 2007). By default, we have set the number of permutations to 1000 (smallest p value approximately equal to 0.001) as this has been shown to produce good results (Landgrebe et al., 2002;Linting et al., 2011). Users of the package should keep in mind that higher numbers of permutations will increase computation time with potentially only a small gain on the approximation. Our simulations indicate that between 500 and 1000 permutations provide a good compromise between computing time and precision in estimating confidence intervals, depending on the data volume ( Figure 4-figure supplement 1). The package implements a single permutation strategy for VAF, the so-called permD (permutation of the entire data set) (Buja and Eyuboglu, 1992;Linting et al., 2011) where variables are permuted independently and concomitantly ( Figure 2A) opposed to permV (permutation of a single variable) (Linting et al., 2011) where variables are permuted one at the time ( Figure 2B). These methods are further discussed on the component interpretation section.
R Code Box 2.
Case study: After performing a PCA, we first determined the number of components that can be regarded as informative. Several criteria can be used as mentioned earlier. Here, we opted for the permutation test of VAF, computed using the permut_pc_test() function (R Code Box 2). We have applied this test to the data using 10,000 permutations. The results show that the three first PCs (PC1, PC2, and PC3) are significantly different from random at an alpha of 0.05 adjusting the p value ( Figure 3), and therefore we will keep these three PCs for subsequent analysis. PC1 accounts for 32.9% of the variance, PC2 18.3% and PC3 9.8%.
Step 3: Component interpretation: what do these components mean?
A key part of the analytical workflow is the interpretation of the main components, where the most relevant PCs can be used to represent the correlation between the original variables as a proxy for multivariate disease patterns. Each component is composed of a weighted combination of all the variables. Some components might be explained by only a few variables with high importance, whereas others might have several variables with important contributions to them. There are a few metrics that can be used for interpreting the relation between the original variables and the PCs (Abdi and Williams, 2010). In the syndRomics package, we use the standardized loadings or correlation vector coefficients (Jackson and Hearne, 1973), and the communalities, which are the sum of squared loadings for each variable across selected PCs representing how much of the variance of each variable can be explained by the total number of kept components. Loadings can be interpreted as the Pearson's r correlation coefficient between a PC and a variable, and it is used to assess the contribution of individual variables on each PC and the direction on which the variable moves along the PC (i.e. opposite or same direction as in the interpretation of a correlation). Communalities can be interpreted as the global impact of a variable in the chosen PCA solution.
In general, the strategy consists of determining a threshold for the absolute value of loadings or the communalities above which variables are considered to have important contribution in the definition of a component or across the chosen PCs. For example, if a threshold of |loading| > 0.2 is chosen, all variables for a given PC with a loading > 0.2 or a loading < À0.2 will be considered to contribute on the PC (aka salient variable). The matter then turns to determining an appropriate threshold. Some somewhat arbitrary rules of thumb for the loadings have been established. However, those have a strong determination in psychological studies and whether they are appropriate in biomedical research has yet to be verified. An alternative 'quasi-inferential' method is to use permutation test as discussed above for PC VAF but testing for metrics of variable contribution such as independently and concomitantly for each variable. For each P sample, a PCA is run and either the loadings, communalities or VAF are calculated. All P PCA solutions form the null distribution for non-parametric hypothesis testing of loadings or VAF. (D) The permutation test algorithm for loadings under permV is performed with and extra step of Procrustes rotation between each of the P samples to the parent component loadings. The P rotated loadings will then form the null distribution for each variable.
loadings (Buja and Eyuboglu, 1992;Peres-Neto et al., 2003) or communalities (Linting et al., 2011). Using resampling strategies, these permutation methods offer data-driven determination of variable importance and contribution, which might reduce subjective biases. Thus, rejecting the null hypothesis for a given metric, variable and PC, suggest that such variable has a contribution onto the construction of the component that is above what is expected by random noise. As in the case of VAF, this establishes a lower bound for |loadings| or communalities below which they should be The graph shows the original VAF for the first five PCs and the average and 95% confidence interval VAF of the permuted PCA distribution (p=10000) using the permD method. * Statistical difference for the non-parametric test at alpha = 0.05 and adjusted p value by BH. The three first PCs were selected for the subsequent analysis. (B) Barmap of the original communalities (bars) and the permuted distribution (permV, p=3000) for each variable calculated over the first three PCs. (C) Barmap of the original loadings (bars) and the permuted distribution (permV, p=3000) for each variable and each of the first three PCs. Solid dotes represent the mean of the permuted distribution and error bars represent the 95% CI. The online version of this article includes the following source data and figure supplement(s) for figure 3: Source data 1. csv file containing the source data for panel A in Figure 3. Source data 2. csv file containing the source data for panel B in Figure 3. Source data 3. csv file containing the source data for panel C in Figure 3.  considered noise. In situations with stable solutions and high signal-to-noise ratio, low |loadings| or communalities might still be statistically significant, but the contribution of the variable should be gauged respect to other variables. In the package, we have incorporated permutation test of the loadings as in Buja and Eyuboglu, 1992;Peres-Neto et al., 2003 that can serve to determine the loading threshold, where the variables are permuted independently and concomitantly ( Figure 2A and C). Linting et al., designed and tested an strategy for the communalities where only one variable is permuted at the time, showing great results in determining the contribution of variables using communalities (Linting et al., 2011). This method has resulted in better determination of the significant contribution of variables on the PCA solution with higher statistical power and proper type I error, and therefore has been incorporated in the package as the default method for both the communalities and the loadings ( Figure 2B and D). Following Linting et al., terminology, users can specify the permutation strategy for the loadings as one variable at the time (permV, as in [Linting et al., 2011]) or as all the variable together (permD, as in [Buja and Eyuboglu, 1992;Linting et al., 2011;Peres-Neto et al., 2003]). See Materials and methods for details on the permutation algorithms. In addition to permutation strategies, the package implements bootstrapping methods for constructing confidence intervals of component loadings and communalities that can also facilitate PCs interpretation (see component stability).
The selection of number of permutations in this case follows similar rationale as described above for the VAF. It is important to note that the minimal number of permutation needed to have enough statistical power and precision will depend on the size of the dataset, both on the number of variables and samples (Buja and Eyuboglu, 1992; Figure 4-figure supplement 1). There is also the Table 2. List of variables included in the second case study.  (Benjamini and Hochberg, 1995) method. As a rule of thumb, these researchers advised to only use multiple testing correction (for FDR) when the data contains at least 20 variables and 100 observations or subjects, and to use the uncorrected p-values otherwise (Linting et al., 2011). p-Value adjustment has been incorporated in the permutation function on the package, with controlling for FDR by BH as default. The reader should be cautioned against overinterpreting or misinterpreting the meaning of a PC. The interpretation can be subjective, and unconscious biases can be reflected on the interpretation of PCs. The tools offered by the package help mitigate potential subjective biases, although data biases will affect the results. Another consideration is that it is possible that some of these metrics seem to 'contradict' each other. For example, there is the possibility that a component has an important contribution to the variance of the data (high VAF) and yet all the loadings be small. Contrary, a component with a small set of high loadings could be considered to be insignificant by permuting its VAF (Buja and Eyuboglu, 1992). As in any analytical approach, domain knowledge is critical for the interpretation of disease components.
Case study: After deciding to keep three components, we studied the communalities and loadings to determine their identity. Here, we applied the permut_pc_test() function (R Code Box 3) setting the argument statistic = 'commun' or 's.loadings' and the perm.method = 'permV' and using the BH method for controlling for FDR. The results of the permutation test on the communalities can be seen in Figure 3B and in Table 3. We can appreciate that all variables are significantly represented by the three chosen PCs, although there are five variables with communality less than 0.5, indicating that the retained PCs only explain 50% of the variance on these variables. In PCA, communalities can suggest which variables do or do not contribute to the extracted components altogether. Considering the loadings, the results for PC1, PC2, and PC3 are shown in Figure 3C and in Tables 4, 5 and 6, respectively. One can appreciate that the cutoff loading for significance at alpha 0.05 using the adjusted p value is approximately |0.21| for PC1, |0.25| for PC2 and |0.4| for PC3. This behavior of different thresholds for significance has been previously described (Buja and Eyuboglu, 1992) and reflects the fact that PCs accounting for less variance might contain more random noise, thus needing a higher loading for a variable to be considered as an important contributor. Loadings are indicative of both strength of association between a variable and a PC and the direction in which they interact. For reading on the interpretation of the loadings and components in this case study, see Ferguson et al., 2013. R Code Box 3 permut_pc_test (pca, pca_data, p=1000, ndim = 3, statistic = 's.loadings', perm. method = 'permV'). permut_pc_test (pca, pca_data, p=1000, ndim = 3, statistic = 'communa', perm. method = 'permV').
Step 4: Component stability: how robust are the components?
The presence of a syndrome or disease pattern, represented by a component, should hold true regardless of variations in experiments or metrics meant to measure that same pattern. For example, two experiments with different subjects but the same collected variables should result in inferentially equivalent components if they are true features of the disease and not experimental artifacts. The sensitivity of PCs to experimental, metric, or other forms of variation is termed 'component stability'.
Components from different PCAs (from different experiments as an example) that are extremely similar are considered to be a stable, and characterizing component stability is important to determine the robustness of the initial PCA (Guadagnoli and Velicer, 1988;Linting, 2007). A robust PC would be largely unaffected by data variations (i.e. low sensitivity). The goal of the stability analysis is to determine such sensitivity. Given that performing multiple replication experiments in biomedical research is not always possible, component stability can be approximated by resampling techniques such as bootstrapping (Babamoradi et al., 2013;Linting et al., 2007a;Timmerman et al., 2007;Zientek and Thompson, 2007). Bootstrap methods for component stability have been extensively studied, but users should be aware of the limitations and advantages of these methods and their performance for component stability depending on the use case (Babamoradi et al., 2013;Guadagnoli and Velicer, 1988;Linting et al., 2007b;Timmerman et al., 2007;Zientek and Thompson, 2007). The package implements functionalities to help study the component stability affected by data selection variability by implementing bootstrapping methods (Babamoradi et al., 2013;Linting et al., 2007b;Timmerman et al., 2007;Zientek and Thompson, 2007) and stability metrics. The default method used in the package is the simple or ordinary bootstrap consisting of generating a new sample that has the same size (i.e. same number of subjects or observations) and same variables as the original data, but where the subjects have been randomly selected from the data with replacement ( Figure 4A). This process is repeated several times (here referred as b times) to generate a sample of bootstrapped data. In the first case example, each of the b bootstrapped samples contain 159 subjects and 18 variables, but one subject might appear more than once and another subject might not show up in a specific sample.
Component stability can be studied at the whole component level or at the level of the individual variables through the loadings and communalities. The package implements component similarity indexes (aka factor matching indexes) (Cattell and Baggaley, 1960;Cattell et al., 1969;Guadagnoli and Velicer, 1991) as metrics to study the stability of PCs. These metrics can be used to determine the similarity between the different bootstrapped samples, to test the validity of the extracted component under two or more experimental conditions, to assess the multidimensional equivalence of two or more replication experiments, or to determine the impact of imputing missing values.
Case study: To understand the sensitivity of variables and components to experimental variations, we used the pc_stability() function with b = 1000 bootstrapped samples (R Code Box 4), setting the sim argument to 'balanced' to perform balanced bootstrapping. The function will return the average of the loadings and the specified similarity metric across all the b samples as well as the specified confidence interval. For this example, the 95% CI (accelerated and bias-corrected, see Materials and methods) and the bootstrapped average can be seen in Figure 5. In general, the original loadings are close to the bootstrapped average which indicates that the results are unbiased. Moreover, the confidence regions for most higher value loadings are reasonably small, suggesting that these loadings are stable to experimental variation. In addition, the similarity metrics for the three PCs suggest component stability, meaning that the composition of the components is also stable.
The accepted values for these metrics indicating stability might vary by field and the metric of interest. Some indicative values are mentioned in the respective method section, but the user should be aware of subjective biases when determining a threshold for considering stability (Lorenzo-Seva and ten Berge, 2006). Finally, the function will return the average and percentile CI of the communalities, which can be used to assess which variables are more stable in the selected PCA solution.
R Code Box 4.
Assessing the impact of imputation methods on introducing noise when dealing with missing data should be considered. As described earlier, we used multiple imputation to generate the dataset for analysis. Multiple imputation generates m complete datasets where the imputed values might vary, but the observed values are the same (in the case study m = 50). We used the stability analysis described above to determine the sensitivity of the PCA solution to variations introduced by imputing missing values. We calculated the similarity metrics between all the 50 imputed datasets for the first three PCs, as well as the loadings. We observed high similarities between the PCs obtained from the imputed datasets (Table 7) and the loadings showed narrow CIs (Figure 3-figure supplement 1). We concluded that multiple imputation has produced stable solutions with acceptable impact on both the component and variables. A future version of the package might include more robust methods for pooling and testing multiple imputation in PCA context (van Ginkel and  Step

5: Component visualization
Communicating the analysis is a necessary part of the workflow. Although we have included this at the end of the use case, visualization can be also used for aiding in component selection, component interpretation and component stability analysis. There are several ways a PCA solution can be visualized. Here, we describe the plots implemented in the syndRomics package. We have coded three types of plots (syndromics plot, heatmap, and barmap) using the grammar of the graphics framework (Wilkinson, 2005) implemented in R by the ggplot2 package. This allow users to customize the plots using the rich landscape of the ggplot2 universe. The syndromic plot was first published by Ferguson et al., 2013 and represents PCs as the center of a Venn diagram ( Figure 1A), consisting of (1) a middle convex triangle displaying the 'variance accounted for' (VAF) for a given PC and (2) radial arrows pointing to the center of the triangle for each variable with a standardized loading above a certain threshold ( Figure 6A-C). The width of each arrow and the color saturation are proportional to the magnitude of the standardized loading they represent. The color of each arrow additionally differentiates between positive or negative loadings (e.g. blue represents a loading of +1, red represents a loading of À1, and white represents a loading of 0).  Syndromic plots are especially useful for conveying PC identity in an easy to understand, concise way for publication. Heatmap and barmap plots are alternative visualizations of the loadings beyond the syndromic plot. The major difference between these two plots and the syndromic plot is that both the barmap ( Figure 5A) and heatmap ( Figure 6D) plots display all variables (or a manually selected subset) instead of only the ones with loadings above a given threshold. The absolute loadings that exceed a cutoff threshold can be noted (e.g. by a star *). Moreover, in the case of barmap plots, the cutoff is represented in the graph by vertical lines. This is particularly useful when there are too many above-threshold variables, which would crowd the syndromic plot visualization, or when comparing loadings between PCs more easily. In addition, barmaps are useful for documenting the results of the resampling procedures since error bars can be used to represent the variation of the metrics over the resampling. The permut_pc_test() and pc_stability() functions return such plots.
Case study (R Code Box 5): We can visualize the three selected PCs using the plotting functions in syndRomics (Figure 6). In this case, we chose to represent PC1, PC2 and PC3 using the syndromics plots ( Figure 6A, B and C, respectably) using a cutoff threshold of |0.45|. Notice that this is higher than the threshold for significance found by the permutation analysis given the high number of variables. The full loading pattern of the three first PCs can be visualized by a heatmap (Figure 6D), where we have chosen a different cutoff for each PC (0.21, 0.25, and 0.4 for PC1, PC2, and PC3 respectively), or a Barmap ( Figures 3B and 5A). Barmaps can be obtained for the loadings (barmap_loadings()) or for the communalities (barmap_commun()).

Case study 2
In the second case study, we used selected variables from the Transforming Research And Clinical Knowledge in Traumatic Brain Injury (TRACK-TBI) pilot study (Yue et al., 2013) that were analyzed previously and made publicly available . The released dataset version contains 586 de-identified human subjects who were enrolled in the TRACK-TBI pilot study and the 26 selected variables previously analyzed . These variables are a subset of brain imaging results, outcome metrics and genetic polymorphism ( Table 2). The goal is to describe patterns of association between these three categories of variables. A noticeable difference between this dataset and the one used in the first case study is that here we are dealing with a mixed type dataset, where some variables are continuous, some nominal and some ordinal. Therefore, we performed a version of nonlinear PCA that allows for the extraction of patterns in these kinds of data. The syndRomics package has been programmed to work with the results of the princals() function from the Gifi R package. The code for this analysis is found in the supplementary material.
Missing data analysis showed an overall 21.2% missingness distributed between the outcomes and genetic polymorphism variables (Figure 7-figure supplement 1). With the exception of 'MRI results' that has high missingness (61.7% of the observations), all imaging variables are complete. 'MRI results' variable was excluded from the analysis. The subsequent test for MCAR suggest that there are 17 different patterns of missingness and that the hypothesis of MCAR can be rejected overall (p-value<0.001). Thus, excluding subjects from the analysis is not justified (Schafer and Graham, 2002;Buuren, 2018). We instead performed 50 multiple imputations using the mice R package as in the first case study. The 50 imputed datasets where then aggregated to perform nonlinear PCA using princals() (see Materials and methods for details).
Permutation test of PC VAF suggests that the first 6 PCs contain information that can be regarded as significant above random chance. Although a deep analysis of these six PCs might be of interest, the first three PCs explain the major variance (25.6%, 10.6%, and 9.8%, respectively). Therefore, we focused on interpreting these for illustration purposes (Figure 7 and Tables 8-10). The first PC significantly loaded highly on two genetic variants in opposite directions (SNP_DRD2 loading = À0.677, SNP_ANKK1_Gly318AR loading = 0.661) as well as outcomes of neuropsychological function at 6 months after TBI (CVLT_long loading = À0.614, CVLT_short loading = À0.542) ( Figure 7A-C). All other variables also significantly loaded on to PC1, but with |loadings|~0.3 (Tables 8-10) suggesting that their contribution in PC1 identity is less important. Lower values in CVLT (California Verbal Learning Test) suggest learning and memory impairments, which are well known after TBI. Given the negative loadings for the included CVLT measures (short and long recall), negative values in PC1 might reflect better CVLT outcomes at 6 months after TBI. The stability of the PC1 pattern to multiple imputation is relatively low, with higher loadings showing high variation (Figure 7-figure supplement 1, Table 11), emphasizing the importance of studying stability of components to multiple imputation. Nonetheless, the bootstrapped loadings were stable ( Figure 7B), and decay in CVLT performance after TBI has been previously associated to polymorphisms in DRD2 and ANKK1 genes (Failla et al., 2015;McAllister et al., 2008;Nielson et al., 2017;Yue et al., 2017), providing literature validation of PC1. The variables with higher positive loadings in PC2 were related to the imaging findings and negative loadings with global function outcomes at 3 and 6 months after TBI (GOSE score) (Figure 7. A, D). Lower scores in GOSE are indicative of lower global function and positive values in imaging findings are suggestive of a bigger or more noticeable brain damage. PC2 presented the higher stability to both resampling and to multiple imputation (Figure 7-figure supplement 1, Table 11). Altogether, PC2 might be interpreted as a surrogate for 'TBI severity', where higher positive values would indicate higher brain damage with less function at 3 and 6 months after injury, a signature described in the previous analysis of this data . Finally, given the instability of PC3, with most loadings being considered non-significant by the permutation test and the high variance to multiple imputation (Figure 7-figure supplement 1, Table 11), PC3 can not be interpreted with certainty, and we should not attempt its explanation.

Discussion
Biomedical research needs more multivariate analytics to help realize the potential of precision medicine. While multiple variables are collected in typical preclinical experiments and clinical trials, univariate statistics continue to be the major analytical and decision-making approaches across the different biomedical fields, narrowing our understanding of the complexity of any disease. With the advent of 'omics', analytical approaches for high-dimensional data have started to become more prevalent for the analysis of biological data. Yet, outside the realm of medical bioinformatics, biomedical research continues to be, for most part univariate. The lack of multivariate approaches in analyzing biomedical data can cause biases and constraints to the interpretation of the results and contribute to the lack of reproducibility and bench-to-bedside translation (Ferguson et al., 2011;Huie et al., 2018).
The extraction of disease space, through the use of multivariate methods, can increase our understanding of complex relationships commonly present in biomedical data while preventing some of the issues associated with an excessive use of univariate analytics such as multiple comparison testing and associated false discoveries by chance (Benjamini and Hochberg, 1995;Ferguson et al., 2013;Krzywinski and Altman, 2014). For example, it is common in biomedicine to measure several behavioral and histopathological outcomes that are analyzed independently at the univariate level. This approach increases the chance of false-positive results due to the accumulation of type I testing errors (Benjamini and Hochberg, 1995;Dunn, 1961;Krzywinski and Altman, 2014). Although there are methods to correct for errors when running numerous tests such as multiple-testing correction, their use in biomedicine outside of bioinformatic analysis is scarce. Even when correcting for multiple testing, performing several univariate analyses limits our understanding since univariate analysis does not allow us to study and infer the relationship between measures that might capture different aspects of the matter of study. In our first example case study, several functional tests can be used to study the recovery of forelimb motor function after cervical spinal cord injury in animal models. Each test further contains multiple measures about particular aspects of recovery. Knowing the relationship between these measures through multivariate approaches can increase our understanding of the matter of study while reducing the burden of multiple testing . Importantly, it is also possible that a single univariate test that does not produce significant results misses true biological effects, while a multivariate analysis including the same variables can find patterns and relationship between variables that are significant. Syndromic analysis is, therefore, a framework that uses multivariate analysis of biomedical data in a holistic way, aiming to reveal interactions within complex (patho)-physiological niches, that would be otherwise challenging to discern. Applying syndromic analysis to biomedical data will help uncover the complex relationships of variables and features that constitute different disease and biological states and ultimately accelerate research toward precision medicine.
The syndRomics package implements several functionalities for the visualization, the interpretation, and the analysis of the stability of principal components to facilitate the extraction and analysis of disease patterns. We have demonstrated its usage, showing the potential of the package to support PCA-based analysis in understanding disease complexity. Although the core functionalities of the package are included, future versions might also implement outputs from other PCA functions   as inputs, such as those from the PCA functions in the FactoMineR package (Lê et al., 2008) or the psych package (Revelle, 2017), allowing for better integration to the PCA landscape in R. In addition, other algorithms of bootstrapping and permutation methods for PCA solutions could be incorporated to increase the options and better adapt to the specifics needs of the user (Hong et al., 2006;Linting et al., 2011;Vitale et al., 2017;Zientek and Thompson, 2007).
Here, we emphasize guidance and tools for robust determination of PCA-based disease patterns. We have incorporated resampling methods aiming to reduce subjective biases and to study the stability and generality of the analysis. Although we have shown the use of these functions in different contexts along the process, much more work can be done to extend syndRomics. For example, we demonstrated the stability of our analysis under multiple imputation, and future research could investigate number of multiple imputations or missing conditions necessary for stable disease pattern detection. In addition, visualization features of syndRomics may be extended to help interpret disease patterns resolved by other multivariate or machine learning tools involving structure coefficients or feature impact scores. The syndRomics resampling methods could also be used to estimate the sample size required for stable PCs in the context of syndromic analysis, allowing for sample planning. The implementations in the package are thus positioned to empower both biological and statistical research toward understanding complex biology and diseases.

Availability and requirements
The code to reproduce this analysis can be found in the supplementary material. The data for the first use case comes from the ODC-SCI (Open Data Commons for Spinal Cord Injury, RRID:SCR_ 016673, http://odc-sci.org), ODC-SCI:26 dataset (https://scicrunch.org/odc-sci/about/odc-sci_26).

Package implementation
The syndRomics package offers two major functionalities for the purpose of aiding in the process of syndromics analysis: (1) visualization functions and (2) functions incorporating resampling methods to determine stability and inference of PCs.

Visualization functions
The visualization functions are: syndromic_plot(), heatmap_loadings(), barmap_loadings(), barmap_commun() and VAF_plot(). For the visualization functions, the user can pass an R data.frame object with the standardized loadings (or other metrics) obtained by running PCA and related multivariate methods in their preferred software. We opted for this approach to avoid requiring specific implementations of PCA. Loadings obtained from any PCA solution can be easily formatted for usage with the syndRomics visualization functions. All functions in the package that takes a data frame as argument use the same format (Table 12): variables are organized as rows, and the first column is called 'Variables' and contains the names of the respective variables. The other columns contain the PC loadings and are named 'PC1', 'PC2', etc. Alternatively, the visualizations can also be obtained from the output of the prcomp() function in the stats package in R (linear PCA) or from the output of the princals() function in the Gifi package in R (non-linear PCA by categorical PCA). Finally, the results from pc_stability() and permut_pc_test() can be passed to the plot() generic function in R as the package incorporate the corresponding S3 method for 'syndromics' class object. The list of arguments for the syndromic_plot() function are presented in the package manual. The syndromic_plot() function will internally call extract_syndromic_plot() function (see utility functions) and return a list of ggplot2 objects containing the syndromic plot for the first ndim PCs. For example, if ndim = 5, a syndromic plot for PCs 1 to 5 will be generated. Another important argument is the cut_off, which determines the threshold of absolute standardized loadings to consider for plotting. This argument is chosen by the user and is required (with no default). Another required argument is VAF in case the syndromic_plot() function is called using a data.frame input. If the output of the prcomp() or princals() functions is used, the syndromic_plot() function extracts VAF internally and the user-defined VAF will be ignored. When required, VAF is a character vector of the form 'XX %","XX%", etc., where XX is the VAF for each PC to plot, starting with the first PC, followed by the second, etc. (e.g. c('60.1%","25.3%") for PC1 and PC2, respectively). An issue we found during the implementation is that the arrow visualization does not display correctly in the R graphical device on Windows machines. Rendering the plot into *.pdf format, for instance using the ggsave() function from the ggplot2 package, solves the problem.

heatmap_loadings(), barmap_loadings() and barmap_commun()
Most of the functionalities described for the syndromic_plot() function also apply for the heatma-p_loading(), the barmap_loading(), and the barmap_commun() functions. A noticeable difference in barmap_loading() is that the function will plot the PCs specified in ndim instead of the first ndim components. For example, if ndim = c (3,4,5), components 3, 4, and 5 will be plotted. This allows for more flexibility on which components to plot, such as isolating a single component (e.g. ndim = 3 will only plot component 3).

VAF_plot()
This function can be used to plot a VAF plot from a prcomp() or princals() output. There are two style options, 'line' or 'reduced'.

Resampling functions
There are two major functions using resampling methods, the permut_pc_test() function that implements nonparametric permutation test for either PC VAF for aiding in component selection or PC loadings and communalities for aiding in component interpretation, and the pc_stability() function that implements bootstrapping of PC loadings for stability analysis. These functions take as input the output of the prcomp() or the princals() functions in R as well as the original dataset used on these functions as inputs. The specific call of prcomp() or the princals() used to obtain the original PCA solution is passed down to the resampling functions in the syndRomics package, ensuring that the same arguments are used for resampling (with the exception of the data argument on the original prcomp() or the princals() call, that will be internally changed for each resampling iteration).

permut_pc_test()
In the syndRomics package, the null distribution for the permutation test is generated by permuting the values of each variable independently and concomitantly several times (permD) or permuting one variable at the time (permV) and re-running the PCA on each permuted sample (Figure 2; Buja and Eyuboglu, 1992;Glorfeld, 1995;Linting et al., 2011). When permV method is selected to measure the impact of permuting on loadings, a step of Procrustes rotation of each loading matrix toward the original loading matrix is added to resolve sign reflection, rotation indeterminacy and component translocation ( Figure 2D, see pc_stability for detailed explanation). This step is not performed when the analysis is performed on the communalities since are invariant to such PCA resampling issues (Linting et al., 2011). Confidence intervals of the permuted distribution (null distribution) are calculated using the (1-a)x100% (percentile) of the distribution (Buja and Eyuboglu, 1992). The function calls the permut_pca_D() or permut_pca_V() utility generic function internally to generate the permuted distribution of the selected metric (either "VAF","s.loadings" or "comuna") using either the prcomp() function for linear PCA or the princals() function for nonlinear PCA implemented as S3 R method for the class "prcomp" or "princals". If "VAF" is specified the permD permutation will be used, ignoring the input of the user on the perm.method argument, returning a matrix containing the VAF for the original PCs, as well as the average and the CI of the permuted VAF distribution. In case "s. loadings" or "communa" are specified, the specified permutation method will be considered (i.e. permD or permV) and the function will return a list of matrices, one for each selected PC, with the original loadings, and the average and CI of the permuted loadings distribution. In both cases, p values are calculated as described in the main text and returned. Adjusted p values using the specified method in the adjust.method argument are also returned.

pc_stability()
Component stability can be studied at the whole component level, known as factor invariance, or at the level of the individual loadings. We have implemented both options in the package. By default, the pc_stability() function returns the average and the accelerated and bias-corrected 95% confidence intervals (CI) of the loadings of the bootstrap distribution (Efron, 1987). Depending on the sample size and the number of chosen resamples, the bias-corrected CI will fail and the percentile (1-a)x100% CI will be returned (with corresponding notification). In addition, component similarity or factor matching metrics can be computed by setting the test_similarity=TRUE, which will call the component_similarity() function. For each of the specified similarity metrics, this function returns the average of the metric and its confidence interval (95% CI by default) by the percentiles of the bootstrap distribution. The confidence level and the CI method for the loadings can be changed by changing the conf and ci_type arguments. The function uses the boot() function for generating the bootstrapped samples and the boot.ci() function for extracting the confidence intervals of the loadings. Both boot() and boot.ci() are from the boot package in R. This allows the use of different bootstrapping strategies such as simple or ordinary bootstrapping (by default) or balanced bootstrapping. The reader is referred to the boot package documentation for more details on the different sim methods. A major problem of performing resampling procedures in PCA is what is known as indeterminacies that can invalidate comparing between bootstrapped samples (Babamoradi et al., 2013;Chan et al., 1999;Linting, 2007;Timmerman et al., 2007;Zabala and Pascual, 2016). Sign reflection refers to the change of sign on the component loadings in a PC given slight variation of the data. In addition, slight data variation can also cause component/factor translocation, the change in the position of a component in the PCA solution (e.g. PC1 shifts to the position of PC2), especially when two components have similar VAF. Another problem on performing PCAs with variations in the data is the possibility of rotation indeterminacy when the PCA solution of a resampled data presents with a different rotation of the original PCA solution. These issues generate artificially biased bootstrapped distributions, potentially invalidating the procedure (Timmerman et al., 2007;Zientek and Thompson, 2007). We have implemented a step of procrustes rotation between the original loadings (target) and the bootstrapped sample, as has been previously demonstrated to be a reasonable method to deal with such issues (Timmerman et al., 2007;Zientek and Thompson, 2007). The Procrustes rotation is obtained by the procrustes() function from the pracma package. The algorithm for bootstrapping the PCA solutions is represented in Figure 4 and implemented in the utility function boot_pca_sample(). The number of bootstrap samples is set to 1000 by default. The user must be careful on setting the number too low, reducing the performance of the approximation (Efron, 1987). However, setting the number of bootstrap samples too high might unnecessarily increase computing time with little gain (Figure 4-figure supplement 1).

Indexes of component similarity
We have included several component similarity indexes for determining component/factor invariance in syndRomics. The function component_similiarity () returns the specified similarity metrics as well as their summary statistics (average and standard deviation, if applicable) from a list of loading matrices (load.list). The argument s_cut_off is used in the calculation of the Cattell's s index (see below) and ndim is used to limit the number of components from which to compute the indexes from. Each index has been programmed in a separate utility function for convenience. Although they are not meant to be manually called, users can call them to calculate any of these metrics for a given set of two component loadings. The similarity_metric argument takes a single character or a vector of characters to specify which metrics to compute. These can be: 'cc_index', 'r_correlation', 'rmse' and/or 's_index'. The user can also specify 'all' to get all metrics. Their definitions are documented below.

Congruence coefficient (CC, 'cc_index')
First suggested by Burt, 1948, it was popularized by Tucker, 1951 and therefore is also known as Tucker's congruence coefficient. It is calculated as (2) where x i and y i are the loadings of the variable i on the component or factor x and y respectively.
x; y ð Þ is equivalent to the cosine of the angle between two vectors and is also referred to as the cosine similarity metric. CC is a measure of proportional similarity between two components, and technically the index has a range from -1 (perfect negative congruence) to 1 (perfect positive congruence). In practice, because the all the loadings of a PC can be multiplied by -1 without changing the meaning of the PC, the absolute value of CC is considered, which correspondingly ranges from 0 to 1. The closer to 1, the more similar the two components are. Chan et al. discussed the 0.9 rule of thumb as an indicator of good matching between PCs (Chan et al., 1999). The application of CC as a similarity metric for factor invariance has been extensively studied (Chan et al., 1999;Lorenzo-Seva and ten Berge, 2006).

Pearson's correlation coefficient (r, 'r_correlation')
The calculation of r between two vectors of component loadings has also been used as a pattern matching metric (Guadagnoli and Velicer, 1991). It is computed as (3): In the syndRomics package, the Pearson's correlation coefficient is calculated using the cor() function of the stats package.
Root mean square error (RMSE, 'rmse') RMSE has also been used as a metric for factor matching (Guadagnoli and Velicer, 1991). It is calculated as the square root of the average squared difference of the loadings of the variables as (4): where n is the number of variables in both components x and y. A RMSE of 0 determines a perfect matching, and therefore the smaller the RMSE is, the more equivalent the two components x and y are.

Cattell's s index ('s_index')
The s index was first suggested by Cattell and Baggaley, 1960;Cattell et al., 1969. It is based on the factor mandate matrix (Cattell and Baggaley, 1960) where loadings are either one if a component is considered to act on a variable, called a salient variable, or 0 if not (forming the hyperplane space). Cattell's suggested an arbitrary ± 0.1 cutoff where variables with loadings outside the cutoff range are removed from the hyperplane and considered to be salient variables. In practice, one might want to alter the threshold depending on the experimental conditions. Any loading inside the cutoff range is then interpreted as having been produced by chance. The s index is calculated from the cross-classification of the common variables of two components/factors: where PS = positive salient variable; H = hyperplane variable; NS = negative salient variable; f ij is the joint frequency. Positive and negative salient variables are variables outside the cutoff range with positive or negative loadings respectively. Pattern matching is determined by comparing the cell frequencies in the cross-classification table. Here we implement the simplified form of calculating s (5): The reader is referred to Cattell and Baggaley, 1960;Cattell et al., 1969;Guadagnoli and Velicer, 1991 for details on reasoning and calculations. s ranges from 1 (perfect similarity) to À1 permut_pca_V.prcomp(), permut_pca_V.princals(). It returns a list of the results of permuting the data, conducting a PCA and extracting either the VAF or the standardized loadings for each P as in Figure 2.

Plot.syndromics()
This function implement the S3 method for plotting 'syndromics' class objects generated by pc_stability() and permut_pc_test() functions using the R generic plot(). It returns specific plots calling the visualization functions implemented in the package.

Nonlinear PCA
Nonlinear PCA by optimal scaling and alternating least square was obtained using the princals() function from the {Gifi} package in R. We specified to analyze all variables with nominal restriction scaling, allowing for non-monotonic transformations, and set a restriction of 3 degrees in polynomial transformations for nonlinearity. The corresponding instruction was: princals(nlpca_data, ndim = ncol (nlpca_data), ordinal = FALSE, degrees = 3, knots = knotsGifi(nlpca_data, type='E')), where nlpca_data is the imputed dataset for case study 2 (see supplementary code for more details).

Missing data analysis
Details on the code are available as supplementary material. Data wrangling for the two case studies was performed using R packages included in the Tidyverse package. Missing pattern visualization were obtained using the naniar (Tierney et al., 2020) R packages. Test for MCAR was performed using the TestMCARNormality() function from the MissMech package (Jamshidian et al., 2014). Multiple imputation was performed using predicting mean matching method available in the mice (Buuren and Groothuis-Oudshoorn, 2011) R package, setting the number of imputations to m = 50. A list of 50 complete datasets were then obtained and processed by PCA as specified in the main text. For each m dataset, the loadings where extracted and rotated using Procrustes rotation (pracma package) toward the average of the imputed datasets. The distributions of loadings and component similarities for the first three PCs where calculated using the syndRomics package as described above.