Consider a collection of K data blocks/matrices, {Xk∈ℝn×pk:k=1,…,K}, where n is the number of subjects and p_k the number of features or variables in the kth dataset. Each data block can be written as X_k = J_k+A_k+E_k, where J_k represents the joint signal common to both data blocks, A_k represents the block-individual signal, which is unique to the k^th data block, and E_k represents full-rank isotropic noise. The JIVE model assumes that each rank-reduced signal matrix X_k−E_k [with rank r_k < min(n, p_k)] can be decomposed into a subspace of ℝⁿ that is common across X_k (the joint subspace) and a subspace that is unique to the kth dataset and orthogonal to the joint subspace (the individual subspaces) (Feng et al., 2018). In our presentation, we expand on one of three ways to represent the joint subspace, called the “common normalized score” representation in Feng et al. (2018). We emphasize this representation because it results in a correspondence between the joint components of each dataset, whereas the other representations are arguably less interpretable. The common basis, Z∈ℝn×rJ, is derived from joint analysis of all data blocks, and the other, Bk∈ℝn×rIk from the part that remains after joint analysis, where r_Ik = r_k−r_J. Let I_d be the d×d identity matrix and 0 a matrix of zeros. Furthermore, let the joint and individual signal matrices of the k^th data block take the form J_k = ZW_Jk and A_k = B_kW_Ik, respectively. Then the JIVE model corresponds to the matrix decomposition

We call Z joint subject scores and W_Jk joint variable loadings. Individual subject scores are given by B_k and individual variable loadings by W_Ik. Intuitively, this decomposition is similar to a singular value decomposition on each dataset but with part of the basis constrained to be equal in the two decompositions.

In this representation, we do not enforce orthogonality between B_k and Bk′. Later, we propose a permutation test for the joint rank, r_J, that determines when the correlation between signal is sufficiently large to be deemed joint, but allows insignificant correlation between individual subject scores. Our proposed approach will also result in an intuitive ordering of components by the strength of evidence that they are joint. Also note that in (1), the rows of the loadings matrices W_Jk are not orthogonal.

For the HCP network data that we examine in Section 4, we can translate each row of the score (Z, B_k) matrix into a low-dimensional vector summary of a participant's kth network data (e.g., FC). The joint scores Z summarize information that is common across modalities, while B_k comprise information unique to an individual modality. For instance, Section 4.3 shows that CJIVE joint scores are more strongly associated with a measure of fluid intelligence than individual scores. The lth row of the loading matrix W_Jk exhibits the magnitude with which network edges contribute to the lth column of the summary scores in Z. In Section 4.4, we examine variable loadings to develop insight into latent structures that are common within both modalities and those which are unique to each.

We conduct simulation studies to address the following gaps in the current understanding of the performance of R.JIVE and AJIVE: (1) accuracy when the joint signal strength is low vs. high; (2) rank selection when the number of joint components is >1; and (3) the impact of the initial signal rank selection on joint rank selection. We use a full factorial design with the following factors:

The joint rank was 3, individual ranks were 2, and n = 200 in all settings. The entries of the error matrices E₁ and E₂ were randomly drawn from a standard Gaussian distribution. The number of features in X₁ and the individual variation explained for both data blocks were held constant at p₁ = 200 and RI12=RI22=0.25, respectively.

Experimental factor 1 (i.e., p₂) allows us to assess the impact of p_k on the accuracy subspace estimation and r_J estimates. Factors 2 and 3 (i.e., RJ12 and RJ22) allow us to examine the impact of the joint signal's magnitude within each dataset.

For each simulation, the subject score matrix [Z, B₁, B₂] was drawn from a Bernoulli distribution, with probability 0.2 for Z and 0.4 for B_k. The use of two values is similar to the toy examples from Feng et al. (2018), which used ±1. Next, we defined loading matrices W_Jk and W_Ik with entries from independent, mean 0 multivariate Gaussian distributions with covariance matrices diag(9, 4, 1) and diag(4, 1), respectively. The values along the diagonals were chosen to ensure the strength of components within each joint/individual signal diminished from first to last. Note that this set-up results in approximately orthogonal A₁ and A₂. In R.JIVE, we use the option enforcing this orthogonality. This set-up may favor the rank-selection procedure in AJIVE since principal angles between A₁andA₂ are large and corresponding singular values are unlikely to exceed the Wedin and random bounds described in Section 2.1.2.

In order to achieve the desired values of RJk2 and RIk2, we rescale the joint and individual matrices such that X_k = d_kJ_K+c_kA_k+E_k for appropriate constants c_k and d_k, as described in Web Appendix B (Supplementary material).

The chordal subspace norm is a distance metric for linear subspaces that has been generalized to matrices, say, true joint scores Zand estimated joint scoresZ^, of possibly different ranks (Ye and Lim, 2014) and can be calculated as

where q=mink[rank(Z),rank(Zk^)] and θ_m are the principal angles between the column space of Z and Z^. We use this metric in our simulation studies to describe the accuracy of JIVE estimates. Note when the column space of Z is contained in the column space of Z^, δ(Z,Z^)=0. Therefore, comparing results from different methods requires examination of rank estimates and subspace estimates.

We performed 100 simulations using the following methods: (1) AJIVE-r_k, where we used the true total number of components r_k (joint rank + individual rank) as input; (2) AJIVE-Over, where the total number of components was chosen to retain 95% of the variance; (3) R.JIVE-Oracle, which uses both the true r_k and r_J as input; and (4) R.JIVE-Free, with its permutation based algorithm for choosing ranks. We also defined (5) CJIVE-r_k and (6) CJIVE-Over using the same approach for total signal ranks and selecting the joint rank using our permutation test with n_perms = 500 and α = 0.05.

To investigate the prediction methods outlined in Section 2.2.3, the subjects for each simulation were randomly divided into training and test subjects, both with sample sizes n/2 = 100. AJIVE, CJIVE, and R.JIVE, all with true signal ranks used as inputs, were applied on the training datasets. Subject scores were predicted for new subjects, represented by the test datasets. We then assessed performance by calculating the Pearson correlation coefficient between predicted joint scores for the test datasets and true joint scores for the same datasets for each of the r_J joint score components.

Figures 2A,B show that CJIVE-r_k and AJIVE-r_k chose the correct joint rank in nearly 100% of simulations in all settings except for the low-signal lower-dimensional case. Further investigation indicated the joint rank selection in AJIVE tends to be driven by the random direction bound, rather than the Wedin bound (see Web Appendix A.1.2 in Supplementary material). AJIVE-Over and CJIVE-Over both routinely underestimated the number of joint components in all scenarios except lower dimensional high-signal case. When an estimate of r_k is very large, the correlation between permuted datasets can be very large, such that zero joint components are significant. The joint rank estimated in R.JIVE is equal to 2 in a majority of simulations when the joint signal in both datasets is relatively large (bottom-right panels in Figures 2A,B: RJ12=RJ22=0.5), while it is mostly 0 or 1 in the other scenarios.

CJIVE-r_k and AJIVE-r_k joint score subspace errors trended less than R.JIVE, CJIVE-Over, and AJIVE-Over in all settings, as shown in Figure 3. Although the chordal distances for joint loading subspaces from R.JIVE trended less than those from AJIVE-r_k, the lack of accurate joint rank estimates from R.JIVE may indicate that estimated subspaces partially lie within true subspaces.

To summarize, we find that CJIVE-r_k and AJIVE-r_k chose the joint rank correctly in most simulations. For both CJIVE-Over and AJIVE-Over, including too many initial signal components generally resulted in a noise-contaminated signal for each data matrix, which resulted in too few joint components or none at all. Moreover, CJIVE-r_k and AJIVE-r_k estimates of joint score and loading subspaces tended to be more accurate than both R.JIVE-Free and R.JIVE-Oracle. In simulations, AJIVE-r_k is equivalent to AJIVE-Scree plot because the signal rank is identified from the scree plots (Web Appendix Figure S2 in Supplementary material).

Out-of-sample subject score estimates were more accurate across joint components using CJIVE-prediction compared to G-inverse prediction (Figure 4). However, R.JIVE-prediction results were most accurate when the joint signal in at least one dataset was relatively strong, i.e., RJk2=0.5 for k = 1or2. The Pearson correlation coefficients tend to be close to 1, on average, for the first joint component of subject scores across all simulation settings. The third component was predicted poorly in all methods when the joint signal is relatively weak, i.e., RJ12=RJ22=0.05. Recall data were simulated so that the proportion of variance attributable to the jth joint component in X_k, k = 1, 2, j = 1, 2, 3 is given by RJk2=(3-(j-1)6). Therefore, components are ordered (from highest to lowest) by the proportion of joint variation that they contribute, which may contribute to the trend in poorer prediction as j increased.

Our data application uses measures of FC and SC from n = 998 study participants (532 females) in the young adult Human Connectome Project (HCP). Web Appendix Table S4 in Supplementary material provides demographics. We applied R.JIVE, AJIVE, CJIVE, and sCCA to examine multivariate relationships across brain networks as measured by Fisher z-transformed correlations from rs-fMRI (FC) and log-transformed streamline counts from dMRI (SC).

HCP rs-fMRI data comprise two left-right phase encoded and two right-left phase encoded 15-min eyes-open rs-fMRI runs (Glasser et al., 2013). Each run used 2-mm isotropic voxels with 0.72 s repetition time. For each run, we calculated the average time series for each of the 68 cortical regions of interest (ROIs) from Desikan et al. (2006) plus the 19 subcortical gray-matter ROIs from Glasser et al. (2013). For each participant and pair of ROIs, the Pearson correlation was calculated, Fisher z-transformed, and then averaged across the four runs. The lower diagonal of each subject's connectivity matrix was vectorized, resulting in p₁ = 3, 741.

For each HCP participant, three left-right and three right-left phase-encoded runs of dMRI from three shells of b = 1, 000, 2, 000, and 3, 000 s/mm² with 90 directions and 6 b₀ acquisitions interspersed throughout were acquired (Glasser et al., 2013). Whole-brain tractography for each participant was conducted using probabilistic tractography as detailed in Zhang et al. (2018). On average, around 10⁵ voxels occurring along the white matter/gray matter interface were identified as seeding regions for each participant. Sixteen streamlines were initiated for each seeding voxel, resulting in ~10⁶ streamlines for each participant. Nodes of the SC networks were defined from the same ROIs as the rs-fMRI. Edges were represented by the number of viable streamlines between ROIs, with viability determined by three procedures: (1) each gray matter ROI is dilated to include a small portion of white matter region; (2) streamlines connecting multiple ROIs were cut into pieces such that no streamlines pass through ROIs; and (3) apparent outliers were removed. Finally, edges where at least 99% of subjects had zero streamlines were removed, and the remaining streamline counts were log transformed. There were p₂ = 3, 330 edges in the resultant SC data matrix. Plots of the mean FC and SC appear in the Web Appendix Figure S5 (Supplementary material).

Both AJIVE and CJIVE with the scree-plot method for choosing total ranks estimated two joint components (Table 2), which implies that results from these methods are equivalent. Similarly, both AJIVE and CJIVE estimated 0 joint components when the total ranks were chosen to result in retaining 95% of the variation. R.JIVE with its permutation tests estimated 1 joint component.

The canonical correlations were ρ₁ = 0.31 and ρ₂ = 0.21 using 1,000 permutations in CJIVE-Scree plot. n CJIVE-Scree plot, we also examined the breakdown of the joint variances by component: the proportion of variation attributable to joint component 1 was 0.094 in FC and 0.017 in SC (Table 2). For component 2, the values were 0.018 and 0.015, respectively.

In addition to the previous analysis, we performed an irregular grid search to examine the impact of the signal rank selection on the estimation of the joint rank when using CJIVE and AJIVE. Specifically, we examined {2, 5, 7, 10, 15, 20, 25, 30, 39, 40, 50, 75, 100, 200, 225, 500, 990} for FC, where 39 and 225 capture 80 and 95% of the variance, respectively, and {2, 5, 7, 10, 15, 20, 25, 30, 40, 50, 75, 100, 200, 330, 500, 683, 990} for SC, where 330 and 683 capture 80 and 95% of the variance. The proportion of variation explained by the corresponding number of PCs are given in Web Appendix Table S3 (Supplementary material). The joint rank estimates for each pair of signal ranks are displayed in Figure 5. The joint rank in CJIVE and AJIVE tended to increase initially. When 80% of the variance was retained in FC and SC (r_FC = 39 and r_SC = 330), CJIVE and AJIVE estimated r_J= 3 and 4, respectively. The joint ranks were maximized at r_FC = 225 and r_SC = 75, with CJIVE selecting r_J = 11 and AJIVE r_J = 12 (Figure 5). The joint rank estimated by AJIVE and CJIVE depends on the choice of total signal rank, but hereafter we focus on the more parsimonious representation from Figure 5, which is easier to interpret.

Joint subject scores from CJIVE-Scree plot, R.JIVE (R.JIVE using permutation tests for both joint and individual signal ranks), and sCCA, and individual scores from CJIVE-Scree plot and R.JIVE were examined for associations with fluid intelligence (gF). In the HCP, gF was measured as the number of correct responses to the Penn Progressive Matrices Test. We selected this variable as it has previously been examined in Finn et al. (2015), Smith et al. (2015), and our prior study Risk and Gaynanova (2021), and no other behavioral variables were examined. Here, AJIVE-Scree plot results are equivalent to CJIVE-Scree plot, since both methods chose two joint components. We used the R-package rsq to calculate the adjusted partial R-squared from the multiple regression predicting fluid intelligence from the joint and individual scores (Zhang, 2022), and then take the square root to obtain the partial correlation coefficients. We also estimated two pairs of canonical variables with sCCA. In order to compare results from sCCA to CJIVE, we averaged canonical variables across datasets to obtain a single subject score vector for each joint component. In sCCA, permutations tests resulted in sparsity parameters equal to 0.1 using the PMA R package (Witten et al., 2009).

Among the joint scores, CJIVE-Scree plot resulted in the highest partial correlation coefficient (r = 0.251). Partial correlation coefficients for individual scores (r = 0.248) and the overall correlation of total scores (joint + individual, r = 0.363) were highest in R.JIVE (Table 3). R.JIVE contained a total of 151 components while CJIVE-Scree plot included 15 components.

In all three methods, only the first joint component and no individual components were significantly associated with fluid intelligence after correction for multiple comparisons (CJIVE-Scree plot: first joint component p = 10⁻¹², Bonferroni corrected for 15 comparisons; R.JIVE: p = 10⁻¹¹ for the joint component, corrected for 153 comparisons; sparse-CCA: p = 10⁻³ for the first joint component, corrected for 2 comparisons).

Since edges from FC and SC networks comprise the features in our input data blocks, loadings are imposed onto symmetric matrices. The sign indeterminacy of the joint loadings for each component was chosen to result in positive skewness. In Figure 6A, we see that there were strong positive loadings throughout the FC. Overall, there was no clear spatial correspondence between FC and SC, and the correlation between loadings was −0.04. Instead, overall higher FC was associated with higher SC in many regions, particularly frontal-frontal and frontal-subcortical, with SC loadings in the opposite direction in certain connections between occipital, parietal, temporal, and subcortical. Plots of the joint loadings for the second component and individual loadings appear in Web Appendix Figures S3–S5 (Supplementary material).

Taking the L1 norm of each row within each loading matrix reduces the number of features to the number of nodes, which provides a more detailed examination of the patterns. In this analysis, we are particularly interested in L1 norms that are large in both the left and right hemispheres, which suggests the loadings are capturing meaningful biological structure. In the FC loadings, Figure 6B shows that the most prominent cortical regions in the first joint component correspond to ROIs from the frontal, occipital, and temporal lobes, with extensive left-right hemispheric correspondence. In the SC loadings, we again see left-right hemispheric correspondence, this time in the parietal and temporal lobes, as well as regions that did not exhibit hemispheric correspondence. L1 norms of subcortical regions (not shown) were large in the left and right accumbens, left caudate, and left putamen in both modalities. Additionally, the right putamen and right caudate were prominent in FC, while both left and right hippocampus were prominent in SC. FC and SC loadings for the individual components are depicted in Web Appendix Figures S4, S5 (Supplementary material). FC individual component 2 has large loadings on cortical to subcortical edges, and component 5 has large subcortical to subcortical loadings. SC individual component 3 has prominent loadings in both subcortical-subcortical and cortical-subcortical edges.

Subjects from the HCP data were split into two sets of equal sample size (n = 499) to examine the reproducibility of our results. We will refer to the first sub-sample as “sample A” and the second as “sample B.” CJIVE-Scree plot found r_J = 1 for both samples, while AJIVE-Scree plot found r_J = 2 for sample A and r_J = 1 for sample B. The correlations between the joint loadings from sample A and B were equal to 0.61 for FC and 0.65 for SC (CJIVE-Scree plot and AJIVE-Scree plot are equivalent). When a second joint component was estimated, the correlation of the FC loadings was 0.29 and the SC loadings was 0.38.

We evaluated the three prediction methods from Section 2.2.3. We compared the predicted joint scores (using the out-of-sample loadings) to those from the scores extracted from a separate analysis of sample B (Figure 7). Pearson correlations between the G-inverse predicted subject scores and CJIVE-Scree plot subject scores were 0.52 and 0.15 for components 1 and 2, respectively. Pearson correlations between subject scores estimated on sample B and those predicted for sample B using R.JIVE-predict were −0.02 and −0.03 for components 1 and 2, respectively. Using CJIVE-prediction, Pearson correlations were 0.67 and 0.22 for components 1 and 2, respectively. Similar results were achieved when CJIVE loadings from sample B data were used to predict subject scores for sample A. Recall that in simulations R.JIVE-prediction tended to outperform other methods when the joint signal was relatively large in at least one data set (Figure 4). Future research should explore the conditions that may favor R.JIVE-prediction vs. CJIVE-prediction.

The computation time of CJIVE-Scree plot including the rank permutation test was 99 s. AJIVE-Scree plot with its joint rank selection took 157 s. Using pre-specified scree plot ranks and joint rank = 2, the run time for R.JIVE was over 4 h. These computation times mirrored those in our simulation study, where, on average, CJIVE was twice as fast as AJIVE and ranged from 2 to 50 times faster than R.JIVE (Web Appendix Table S1 in Supplementary material).