Let us consider two random vectors X = ( X 1 , … , X p ) ⊤ and Y = ( Y 1 , … , Y q ) ⊤ , where X contains the composition of p taxa and Y is a q -dimensional vector of non-compositional covariates. The nature of the composition makes X lie in a ( p − 1 ) -dimensional positive simplex. To address the compositionality, Aitchison and Bacon-Shone (1984) proposed applying the log-ratio transformation to compositional covariates resulting in Z / p = ( log ( X 1 / X p ) , … , log ( X p − 1 / X p ) ) , where X p is chosen as the reference component. CCA for compositional data can be formulated to find canonical coefficients a = ( a 1 , … , a q ) ⊤ and b − p = ( b 1 , … , b p − 1 ) ⊤ so that the correlation between a ⊤ Y and b − p ⊤ Z / p is maximized. Note that b − p ⊤ Z / p = ∑ j = 1 p − 1 b j ⁡ log ( X j / X p ) = ∑ j = 1 p b j ⁡ log ( X j ) with b p = − ∑ j = 1 p − 1 b j . To avoid the choice of a reference component, we can write the term b − p ⊤ Z / p in a symmetric form by noticing that b − p ⊤ Z / p = b ⊤ Z , where Z = ( log ( X 1 ) , … , log ( X p ) ) and b = ( b 1 , … , b p ) with b ∈ B p ≔ c = c 1 , … , c p ⊤ : ∑ j = 1 p c j = 0 . Therefore, the compositional CCA aims to find a ̃ and b ̃ such that a ̃ , b ̃ = argmax a ∈ R q , b ∈ B p Corr a ⊤ Y , b ⊤ Z = argmax a ∈ R q , b ∈ B p Cov a ⊤ Y , b ⊤ Z Var a ⊤ Y var b ⊤ Z . When the dimensions p and q are high (as compared to the sample size), regularization is required to encourage sparsity and to obtain a unique solution to the optimization problem. We let Σ Y Z = Cov ( Y , Z ) , Σ Y = Cov ( Y , Y ) and Σ Z = Cov ( Z , Z ) . Define the weighted l 1 norm based on a vector of non-negative weights w = ( w 1 , … , w p ) for a vector b as ‖ b ‖ 1 , w = ∑ j = 1 p w j | b j | . The compositional sCCA problem can then be formulated as max a ∈ R q , b ∈ R p a ⊤ Σ Y Z b s.t. a ⊤ Σ Y a ≤ 1 ,   ‖ a ‖ 1 ≤ C a ,   b ⊤ Σ Z b ≤ 1 ,   ‖ b ‖ 1 , w ≤ C b ,   b ∈ B p . (1)

Here C a , C b > 0 are some positive tuning parameters that control the global shrinkage level. The weight w j allows different penalization strengths according to the data or prior structure information. We will elaborate it in Section 2.3.

It has been shown that in high dimensions, treating the covariance matrix as diagonal can yield good results. Following the same strategy adopted by many of the existing high-dimensional CCA algorithms (e.g., Witten et al., 2009), we substitute in the identity matrix for Σ Y and Σ Z in the CCA formulation (Equation 1). Moreover, we write the (weighted) l 1 constraints on a and b in the Lagrangian form. Given a set of n samples { X i , Y i } i = 1 n , let Σ ̂ Y Z be the sample cross-covariance between Y and Z . We formulate the feasible CCA problem as min a ∈ R q , b ∈ R p − a ⊤ Σ ̂ Y Z b + λ a ‖ a ‖ 1 + λ b ‖ b ‖ 1 , w s.t. ‖ a ‖ 2 ≤ 1 ,   ‖ b ‖ 2 ≤ 1 ,   b ∈ B p , which can be solved by iteratively optimizing the objective function with respect to one parameter while fixing the other parameter. Specifically, we have the following two updating steps.

1. Update b : Fix a ( t ) and update b through

In this section, we discuss the updates in Equations 2, 3. Define the operator g h , λ , w = argmin b ∈ B p 1 2 ‖ h − b ‖ 2 2 + λ ‖ b ‖ 1 , w .

By exploring the Karush-Kuhn-Tuchker conditions, we obtain the following result.

Proposition 2.1. Set b ̆ = arg min b ∈ B p − ( a ( t ) ) ⊤ Σ ̂ Y Z b + λ b ‖ b ‖ 1 , w . The solution to (2) is given by b t = b ̆ , if   ‖ b ̆ ‖ 2 ≤ 1 , g Σ ̂ Y Z ⊤ a t , λ b , w ‖ g Σ ̂ Y Z ⊤ a t , λ b , w ‖ 2 , if   ‖ b ̆ ‖ 2 > 1 .

Remark 2.1. We employ the augmented Lagrangian method (ALM) to solve the optimization problem in g ( h , λ , w ) . Specifically, the ALM involves the following two steps of iterations b t + 1 ← argmin b 1 2 ‖ h − b ‖ 2 2 + λ b ‖ b ‖ 1 , w + μ 1 2 1 ⊤ b + d t 2 , d t + 1 ← d t + μ 2 1 ⊤ b t + 1 , where μ 1 , μ 2 > 0 are the step sizes in dual gradient ascent, which are set to be 1 in our numerical studies. The optimization problem in the first step can be solved using coordinate descent in an inner loop by iterating across the following p components b j t + 1 , r + 1 = 1 1 + μ 1 S h j − μ 1 ∑ i < j b i t + 1 , r + 1 + ∑ i > j b i t + 1 , r + d t , λ b w j , where S ( a , λ ) = sign ( a ) ( | a | − λ ) + denotes the soft-thresholding operator and h j is the j th component of h .

Using similar arguments, we can show that the solution to (3) is given by a t + 1 = a ̆ , if   ‖ a ̆ ‖ 2 ≤ 1 , S Σ ̂ Y Z b t , λ a S Σ ̂ Y Z b t , λ a 2 , if   ‖ a ̆ ‖ 2 > 1 , where a ̆ = a r g m i n a ∈ R q − a ⊤ Σ ̂ Y Z b ( t ) + λ a ‖ a ‖ 1 and S ( a , λ ) = ( S ( a 1 , λ ) , … , S ( a p , λ ) ) ⊤ . Set v = ( v 1 , … , v q ) ⊤ = Σ ̂ Y Z b ( t ) . The objective function in the definition of a ̆ becomes ∑ j = 1 q ( − a j v j + λ a | a j | ) . We see that a j ( t + 1 ) = 0 if λ a ≥ v j , and | a j ( t + 1 ) | = ∞ if λ a < v j . Therefore, we have a t + 1 = 0 , if   λ a ≥ ‖ v ‖ ∞ , S Σ ̂ Y Z b t , λ a S Σ ̂ Y Z b t , λ a 2 , if   λ a < ‖ v ‖ ∞ .

As our goal is to find two directions a and b to maximize Cov a ⊤ Y , b ⊤ Z , we only consider the updates when the l 2 constraints are binding, which leads to Algorithm 1 below.

Algorithm 1.

Compositional sCCA: compositional data versus non-compositional data.

1. Initialize a ( 0 ) as the first left singular vector with unit l 2 norm from the singular value decomposition of Σ ̂ Y Z .

To select the regularization parameters λ a and λ b , we consider a two-stage K -fold cross-validation (CV) method as motivated by Chen et al. (2013). We partition all the samples into M folds, and denote Y k = Y I k and Z k = Z I k , where I k are the indexes of samples in the k th fold for k = 1 , … , K . The K -fold cross-validation criterion is CV λ a , λ b = 1 K ∑ k = 1 K Corr a ̂ − k λ a , λ b ⊤ Y k , b ̂ − k λ a , λ b ⊤ Z k (6) where a ̂ − k ( λ a , λ b ) and b ̂ − k ( λ a , λ b ) are the solutions to the compositional sCCA problem based on the samples ( ∪ k = 1 K I k ) \ I k with the tuning parameters ( λ a , λ b ) . The parameter selection based on the CV criterion can be influenced by shrinkage problems arising from the sparsity penalty. To avoid this bias, we adopt the coefficients estimated from a two-stage approach when evaluating the CV criterion. In the first stage, we implement Algorithm 1 with the given tuning parameter pair ( λ a , λ b ) and exclude variables with zero coefficients. In the second stage, we recalculate the coefficients by applying Algorithm 1 with a tuning parameter pair (0,0). The optimal tuning parameter pair is chosen as the one that maximizes the CV value with these recalculated coefficients. Our approach accounts for the compositional structure in the second stage by restricting the coefficients of compositional data to B p . As will be shown below, the K -fold cross-validation performs reasonably well with K = 5 in our numerical studies.

We briefly describe an extension to the case, where both X and Y are compositional. For example, we want to associate the composition of the bacterial taxa with that of the fungi taxa. Let Y = ( Y 1 , … , Y q ) ⊤ be the relative abundances of another set of compositional features. In this case, we let U = ( log ( Y 1 ) , … , log ( Y q ) ) ⊤ and define Σ ̂ U Z as the sample covariance between U and Z . Following the derivation in Section 2.1, we formulate the two-sided compositional sCCA problem as min a ∈ R q , b ∈ R p − a ⊤ Σ ̂ U Z b + λ a ‖ a ‖ 1 , w 1 + λ b ‖ b ‖ 1 , w 2 s.t. ‖ a ‖ 2 ≤ 1 ,   ‖ b ‖ 2 ≤ 1 ,   a ∈ B q , b ∈ B p , where w j = ( w 1 j , … , w p j ) for j = 1,2 are non-negative weights. This problem can be solved by Algorithm 2 below. We use the two-stage CV criterion to select the tuning parameters in the same way as described in Section 2.3.

Algorithm 2.

Compositional sCCA: compositional data versus compositional data.

1. Initialize a ( 0 ) as the first left singular vector with unit l 2 norm from the singular value decomposition of Σ ̂ U Z .

Based on the setups described in Section 2.1, our procedure is to translate the auxiliary information into different penalization strengths through the weights w . Our framework is general enough to incorporate different types of external structures. For instance, co-expressed genes can be classified into the same group to reflect their biological relationships. In our study, we focus on leveraging the taxonomic grouping structure among taxa. Taxonomically related taxa, such as multiple species within the same genus, tend to share biological traits. Consequently, we anticipate that these taxa will exhibit similar relationships with omics features, leading to the expectation of comparable CCA coefficients. We translate the grouping structure information into different restrictions on the weights. Specifically, we divide the taxa into different groups according to their taxonomy such as phylum, family, and genus. We next consider the set of weights: M Group = { w ∈ 0 , C U p : w i = w j   if   i , j ∈ S d   for   i , j ∈ 1,2 , … , p   and   d ∈ 1,2 , … , D } , where D represents the number of groups and C U denotes an upper bound on the weights.

Following Pramanik and Zhang (2020), we impose a penalty term on the weights and propose an algorithm to jointly estimate weights and parameters. Specifically, we define h w j ; γ = exp w j 1 − 1 γ 1 − 1 γ , if  0 < γ < 1 , w j , if  γ = 1 .

We estimate ( a , b ) and w jointly by solving the following problem min a ∈ R q , b ∈ R p , w ∈ M − a ⊤ Σ ̂ Y Z b + λ a ‖ a ‖ 1 + λ b ∑ j = 1 p w j | b j | − log ⁡ h w j , γ s.t. ‖ a ‖ 2 ≤ 1 ,   ‖ b ‖ 2 ≤ 1 ,   b ∈ B p .

The design of the function h is to reduce our method to the classic (iterative) adaptive Lasso when there is no external information. The readers are referred to Pramanik and Zhang (2020) for more discussions on the motivation.

Next we introduce the algorithm to solve the above problem. We focus on the update for w as the updates for a and b remain the same as in Section 2.2. In particular, we update w through w t + 1 ← argmin w ∈ M ∑ j = 1 p w j | b j t | − log ⁡ h w j , γ .

When M = M Group , it is straightforward to verify that w j t + 1 = C U  if  b j t = 0  for all  j ∈ S d , 1 | S d | − 1 ∑ j ∈ S d | b j t | γ  otherwise .

If we do not have any prior structural information on b that we can take advantage of, we take M to be [ 0 , C U ] p . In this case, we have w j t + 1 = C U  if  b j t = 0 , | b j t | − γ  otherwise .

Algorithm 3 summarizes the implementation details of the structure adaptive Compositional sCCA. The selection of tuning parameter pair ( λ a , λ b , γ ) follows a similar approach as described in Section 2.3, with the difference of using ( λ a , λ b , γ ) in the first stage and ( 0,0 , γ ) in the second stage.

Algorithm 3.

Structure-Adaptive Compositional sCCA: compositional data versus non-compositional data.

1. Initialize a ( 0 ) as the first left singular vector with unit l 2 norm from the singular value decomposition of Σ ̂ Y Z .

We first consider the CCA problem between compositional data (i.e., microbiome data) and non-compositional data (e.g., metabolomics data) following a similar setting considered in Chen et al. (2013). To capture the dependence between the two sets of high-dimensional data, we use a latent variable model to generate the compositional variables { X i } (log scale) and non-compositional variables { Y i } (original scale), where the dependence between these two sets of variables is governed by a latent variable ν . Specifically, we assume that log X i = ν i ω X + ε X , i , Y i = ν i ω Y + ε Y , i , where ν i ∼ N ( 0 , σ ν 2 ) and ε X , i , ε Y , i follow N ( 0 p , σ ε 2 I p × p ) and N ( 0 q , σ ε 2 I q × q ) , respectively. The coefficients ω X ∈ R p and ω Y ∈ R q control the relative contributions of individual variables to the overall association. The ratio σ ν / σ ε determines the overall association strength between log ( X ) and Y , with a larger value indicating stronger association. For the dimensions, we set ( p , q ) = ( 100,100 ) , ( 200,200 ) . We consider two setups for ω X .

S1 ω X = 0.85 10 × ( 1,1,1,1,1,1,1,1,1 − 9 , 0 p − 10 ) ⊤ ;

The first two groups contain both zero and nonzero entries reflecting the fact that the external structure information is imperfect and noisy. We set ω Y = 0.85 × ( 0.08 , 0.084 , 0.089 , … , 0.12 , 0 q − 10 ⊤ ) ⊤ . Next, we fix σ ε = 1 and vary σ ν within { 1,2 , … , 8 } to control the strength of the canonical correlation. We report the true positive rate (TPR), false positive rate (FPR), Matthew’s correlation coefficient (MCC), and Precision to measure the performance of different methods. Here, TPR = TP TP + FN , FPR = FP FP + TN , MCC = TP × TN − FP × FN TP + FP TP + FN TN + FP TN + FN , Precision = TP FP + TP , where TP, FP, TN, and FN represent the true positives, false positives, true negatives, and false negatives, respectively. The TPR, FPR, MCC, and Precision are computed by averaging over 100 simulation replicates. Denote the estimated canonical coefficients by a ̂ and b ̂ . Their estimation targets are ω X / ‖ ω X ‖ 2 and ω Y / ‖ ω Y ‖ 2 , respectively, where this normalization is to ensure comparability. The estimation accuracy is evaluated using the root mean square error (RMSE).

We compare the performance of the following four methods.

1. sCCA: sCCA without considering the compositional effect;

For AC-sCCA and SAC-sCCA, we also apply adaptive weights on a with M = [ 0 , C U ] p in implementation. The value of C U is set to 1 0 5 .

Table 1 summarizes the results for the above four methods when fixing σ ν = 4 . For both Setups S1 and S2, C-sCCA outperforms sCCA in terms of all four measures, especially in reducing the false positive rates and increasing the precision in identifying relevant compositional components, demonstrating the advantage of taking into account the compositional constraint. Compared to the first two methods, AC-sCCA and SAC-sCCA further reduce the FPR and thus lead to higher MCC in estimating a and b . For identifying b in Setup S1, SAC-sCCA outperforms the other three methods by exhibiting higher TPR, nearly zero FPR, and thus higher MCC because of incorporating grouping information. Figure 1 is in general consistent with these findings. As association strength increases, the TPR, MCC, and Precision of C-sCCA, AC-sCCA, and SAC-sCCA increase, whereas their FPRs show a declining trend. When σ ν ≥ 3 , the precision and FPR in estimating b of sCCA becomes worse as association strength increases, which means sCCA identifies more true variables at the cost of including more false variables. Figure 2 presents the RMSE in estimating the canonical coefficients, which decreases as the association strength σ ν increases. By accounting for the compositional effect, the C-sCCA, AC-sCCA, and SAC-sCCA outperform sCCA, with SAC-sCCA providing the most accurate estimation. This demonstrates that methods accounting for the compositional nature yield more accurate estimations than those that do not, and considering structural information can further enhance performance.

Finally, we examine the scenario where p = 100 and q = 200 to assess the performance of our method on unbalanced datasets. As presented in Table 2, the four methods exhibit similar performance to that observed in the case where p = q . The results indicate that our method successfully handles unbalanced dimensions.

In this section, we study the performance of the proposed compositional sCCA for the association between two compositional datasets, for example, bacterial taxa abundance vs fungi taxa abundance. We modify the setting in Section 4.1 by considering the following models log X i = ν i ω X + ε X , i , log Y i = ν i ω Y + ε Y , i , where ν i ∼ N ( 0 , σ ν 2 ) and ε X , i , ε Y , i follow N ( 0 p , σ ε 2 I p × p ) and N ( 0 q , σ ε 2 I q × q ) , respectively. We again consider two setups for ω X and ω Y .

S3 ω X = 0.85 10 × ( 1 9 ⊤ , − 9 , 0 p − 10 ) ⊤  and  ω Y = 0.85 × ( 0.08 , 0.085 , 0.09 , … , 0.12 , − 0.9 , 0 q − 10 ⊤ ) ⊤ ;