Rank-adaptive covariance testing with applications to
genomics and neuroimaging

In biomedical studies, comparing differences between covariance often offers scientific insights beyond what is inferred by comparing the mean differences. Especially, it can help determine if complex joint behavior differs between two groups of samples. In this paper, we present two specific applications of two-sample covariance testing within genomics and neuroimaging.

In genomics, gene expression networks, quantified by the covariance of expression levels of multiple genes, provides a better understanding of the genetic drivers of cellular behavior, particularly when the individual contribution from numerous genes is weak, but their joint co-expression levels are strong (Eisen \BOthers., \APACyear1998; Stuart \BOthers., \APACyear2003). Testing for differences in the gene expression networks between tissue types (e.g., a tumor tissue and normal tissue), molecular subtypes of a cancer (e.g., basal and HER2 subtypes in breast cancer), or cancer types (e.g., breast cancer versus ovarian cancer) can be used to identify important genomic biomarkers of complex diseases via these co-expression levels.

In the second motivating example, we consider the ‘batch effect’ problem in genomics, microbiome, and neuroimaging where non-biological variations are induced by collecting data from multiple lab/sites, platforms, preprocessing methods, or scanners. Within neuroimaging, a number of methods have been proposed to estimate and remove these effects (Hu \BOthers., \APACyear2023), and have the potential to increase reliability of scientific findings from the increased sample sizes and more diverse groups of subjects that come from combining datasets. Recently Chen \BOthers. (\APACyear2022), Zhang \BOthers. (\APACyear2023), and Zhang \BOthers. (\APACyear2024) have empirically observed batch effects reflected in the the heterogeneity of covariances of observations taken at different sites or using different scanners. Despite these observed differences, there has been limited work on testing whether these observed covariance heterogeneities are statistically significant, and hence whether preprocessing the data to mitigate batch effects in covariances is justified.

In both genomics and neuroimaging, techniques which leverage the low-rank structure of the often high dimensional data have been found to be useful in characterizing and understanding variations within the data. These low-rank structures are correspondingly reflected in the spiked structure of the singular values of the covariance matrix. Using tools from differential co-expression analysis (Chowdhury \BOthers., \APACyear2019), low-rank structures have been empirically observed in the differences between gene expression networks (Amar \BOthers., \APACyear2013; Meng \BOthers., \APACyear2018; Yu \BBA Luo, \APACyear2021). In addition, several model-based approaches in cancer genomics demonstrated the utility of leveraging low-rank structures in studying differences between tumor types (Park \BBA Lock, \APACyear2020; Lock \BOthers., \APACyear2022). As well, the batch effect induced covariance heterogeneities in neuroimaging appear to be driven by differences in low-rank structures (Veronese \BOthers., \APACyear2019; Carmon \BOthers., \APACyear2020; Chen \BOthers., \APACyear2022; Zhang \BOthers., \APACyear2023). When these low-rank structures differ between two groups, it is expected that a low-rank structure will explain most of the differences in covariances between two groups. Therefore, there is potential utility in constructing statistics which leverage the low-rank structure of the difference in covariance matrices between two groups in improving statistical power.

In both of the aforementioned application areas, a low-rank structure is likely to be useful in characterizing the differences in covariance matrices between populations, however to date few two-sample covariance testing methods, which test $H_{0}:\Sigma_{1}=\Sigma_{2}$, have been developed which explicitly are able to leverage this type of low-rank structure. Schott (\APACyear2007) and Li \BBA Chen (\APACyear2012) consider test statistics based on the squared Frobenius norm of $\Sigma_{1}-\Sigma_{2}$. Srivastava \BBA Yanagihara (\APACyear2010) develops an estimator based on the trace of $\Sigma_{1}$ and $\Sigma_{2}$ as well as $\Sigma_{1}^{2}$ and $\Sigma_{2}^{2}$. These tests which involve the Frobenius norm or trace can be seen as a function of all of the singular values of $\Sigma_{1}$ and $\Sigma_{2}$ or $\Sigma_{1}-\Sigma_{2}$, and hence do not exploit the data’s potential low-rank structure. Tests which are more powerful in specific low-rank settings include: Cai \BOthers. (\APACyear2013) a test powerful against sparse alternatives based on the maximum of standardized elementwise differences of sample covariance matrices between two groups, Danaher \BOthers. (\APACyear2015) a biological pathway inspired test which uses the leading eigenvalues and trace of the sample covariance matrices, Zhu \BOthers. (\APACyear2017) which uses the sparse leading eigenvalue of $\Sigma_{1}-\Sigma_{2}$, He \BBA Chen (\APACyear2018) whose test statistic is based on differences in superdiagonals between $\Sigma_{1}$ and $\Sigma_{2}$, which is particularly powerful when $\Sigma_{1}$ and $\Sigma_{2}$ have a bandable structure, and Ding \BOthers. (\APACyear2024) which compares a subset of eigenvalues between the two sample covariance matrices.

A selection of the two-sample methods listed in Section 1.3 leverage some form of low-rank structure in either the covariance matrices, or difference in covariance matrices, between groups; however, they individually all make strong assumptions of what the form of the low-rank structure that exists is. This can be problematic in the two application areas this paper examines where it is unclear in advance what low-rank structures exist, let alone which low-rank structures to use in order to maximize power for two-sample covariance matrix testing. In this paper we bridge this gap and propose a two-sample covariance testing method which is able to adapt to the form of the low-rank structures which exist within the data, without placing stringent prior assumptions on their form. Specifically, we construct a test statistic that is adaptive to the form of the low-rank structure found in the difference of two sample covariance matrices, and utilize a permutation scheme to ensure strict Type I error control in the finite sample setting.

The remainder of the paper is organized as follows. In Section 2 the adaptive test statistic underlying RACT is presented, and its permutation procedure introduced. Section 3 investigates the asymptotic behavior of the individual test statistics included in RACT’s adaptive test statistic, and discusses the signal-to-noise trade-off of these individual test statistics. A simulation study is presented in Section 4, demonstrating RACT’s adaptivity to various forms of covariance differences. We then apply RACT to the two application areas of interest in Section 5, and compare its performance to a selection of competing two-sample covariance testing methods. Finally in Section 6 we discuss potential extensions, and limitations, of RACT.

Denote $X_{1}^{(1)},\dots,X_{n_{1}}^{(1)}$, $X_{1}^{(2)},\dots,X_{n_{2}}^{(2)}\in\mathbb{R}^{p}$ as observed features from two groups of sample sizes $n_{1}$ and $n_{2}$, respectively, with equal population means (for practical purposes as a preprocessing step the data can be centered using the sample mean for each group), and let $(\Sigma_{1},\Sigma_{2}),(\hat{\Sigma}_{1},\hat{\Sigma}_{2})$ be the population and sample covariances of these groups. In addition, we define $n=n_{1}+n_{2}$ to be the total number of samples and $\hat{\Sigma}=((n_{1}-1)\hat{\Sigma}_{1}+(n_{2}-1)\hat{\Sigma}_{2})/(n-1)$ to be the pooled covariance. In this section, we test for differences in the covariance matrix, however if one preferred to test for differences in the correlation matrix this could be accomplished by redefining $(\Sigma_{1},\Sigma_{2}),(\hat{\Sigma}_{1},\hat{\Sigma}_{2})$ as population and sample correlation matrices. We define our null and alternative hypotheses as $H_{0}:\Sigma_{1}=\Sigma_{2}$ versus $H_{1}:\Sigma_{1}\neq\Sigma_{2}$.

Testing $H_{0}:\Sigma_{1}=\Sigma_{2}$ requires quantifying the difference between $\hat{\Sigma}_{1}$ and $\hat{\Sigma}_{2}$, for example by using a probabilistic model. As noted in Section 1.3, several existing two-sample covariance testing methods are based on test statistics which utilize the Frobenius norm or leading singular values for detecting differences in covariance structures. However, using either of these test statistics individually would be limited for certain covariance differences and lead to underpowered results. For example, a test statistic based on only a single singular value would be underpowered when several singular values drive the difference in covariance. On the other hand, the Frobenius norm, constructed by $||\hat{\Sigma}_{1}-\hat{\Sigma}_{2}||_{F}=\sqrt{\sum_{r=1}^{p}\sum_{s=1}^{p}(\hat{\Sigma}_{1}[r,s]-\hat{\Sigma}_{2}[r,s])^{2}}$, is well-suited when each entry of $\hat{\Sigma}_{1}-\hat{\Sigma}_{2}$ has non-zero expectation, but it could be underpowered when the entries with non-zero expected values are sparse. Also, even when the signals are dense, considering that the Frobenius norm of a matrix is equivalent to the square root of the sum of squares of all singular values of the matrix, it may be underpowered when $\hat{\Sigma}_{1}-\hat{\Sigma}_{2}$ is low-rank.

For a fixed $k\in\mathbb{N}$, one way to characterize the difference between sample covariance matrices is through the use of a Ky-Fan($k$) norm defined by

where $\|A\|_{(k)}=\sum_{l=1}^{k}\sigma_{l}(A)$ and $\sigma_{l}(A)$ is the $l$th largest singular value of a matrix $A$. If $\Sigma_{1}-\Sigma_{2}$ is low-rank, there likely will exist a $k<p$ such that this Ky-Fan($k$) norm captures most of the variation of the signal, and by excluding the bottom $p-k$ singular values, ignores noise introduced through finite sample variability.

Although $T_{k}$ is a well-motivated test statistic with a prespecified $k$, it is unclear in advance what value of $k$ will maximize $T_{k}$’s power. If $k$ is chosen too small, then $T_{k}$ may fail to capture the signal which exists outside of the top-$k$ singular values of $\hat{\Sigma}_{1}-\hat{\Sigma}_{2}$. Alternatively, if $k$ is chosen too large then some of the singular values included in $T_{k}$ will be noise, decreasing the signal to noise ratio of $T_{k}$. Section 3 further discusses this signal-to-noise trade-off in the asymptotic setting.

The problem of selecting the ‘optimal’ $k$ motivates the proposed method, rank-adaptive covariance testing (RACT), which considers the adaptive test statistic:

where $\mathcal{K}=\{1,\dots,K\}$ represents the collection of Ky-Fan($k$) norms from 1 to $K$, and $\text{E}_{H_{0}}[T_{k}]$, $\text{Var}_{H_{0}}[T_{k}]$ are the expected value and variance of $T_{k}$ under $H_{0}$.

As discussed in Section 3, for normal data in the asymptotic setting, we show that the signal-to-noise ratio of $T_{k}$ is a function of $k$, $\Sigma_{1}$, and $\Sigma_{2}$. Therefore, when $\Sigma_{1}$ and $\Sigma_{2}$ are unknown it is difficult to ascertain in advance which $k$ will attain the maximum signal-to-noise ratio for $k\in\mathbb{N}$. By taking a maximum across different values of $k$ for an appropriately normalized $T_{k}$, under the alternative the behavior of $T_{\text{RACT}}$ will resemble the behavior of the signal-to-noise maximizing $T_{k}$, and similarly the power of $T_{\text{RACT}}$ will be close to the power of this $T_{k}$ (our simulations in Section 4 reflect this adaptivity). By including a diverse set of norms in $T_{\text{RACT}}$, an investigator is freed from having to make a potentially consequential decision as to what single matrix norm should be used for testing $H_{0}$, while at the same time potentially benefiting from the inclusion of norms which are sensitive to certain types of structures in $\Sigma_{1}-\Sigma_{2}$.

The maximum $K$ can be chosen either as $\min(n,p)$, or a smaller value that reflects prior knowledge of the data or computational considerations (since truncated SVD for the top $K$ singular values is an $O(K\times\min(n,p)^{2})$ operation). In this paper, we choose $K$ to be the smallest $K\leq\min(n,p)$ such that the variation of $\hat{\Sigma}$ explained by its top $K$ singular values exceeds 80%. This is done as the bottom 20% of the singular values are unlikely to provide much of the signal for our test statistic in the high-dimensional setting, and by including observations from both groups we ensure $K$ is the same for all permutations of the data.

We use permutation to both estimate $\text{E}_{H_{0}}[T_{k}]$ and $\text{Var}_{H_{0}}[T_{k}]$, as well as to set a critical value for conducting hypothesis testing based on $T_{\text{RACT}}$. The use of permutation to set a critical value is attractive due to the dependencies inherent in the $T_{k}$ statistics. Even if the asymptotic distribution $T_{\text{RACT}}$ was known, the critical values from this distribution could be inaccurate for Type I error control for a finite $n$. We create $B$ permuted datasets (permuting subjects between the two groups) and take the empirical means and standard deviations of $T_{k}$ for all $k\in\mathcal{K}$ to estimate $\text{E}_{H_{0}}[T_{k}]$ and $\text{Var}_{H_{0}}[T_{k}]$. Then, using these same $B$ permuted datasets, we set the critical value at level $\alpha$ of $T_{\text{RACT}}$ for testing $H_{0}$, as the $(1-\alpha)$th quantile of $\{T_{\text{RACT}}^{(1)},\dots,T_{\text{RACT}}^{(B)}\}$ where $T_{\text{RACT}}^{(b)}$ is calculated in the same way as in (2) except each $T_{k}$ is calculated using the permuted dataset. We then reject $H_{0}$ if $T_{\text{RACT}}$ is larger than this critical value.

We conducted extensive simulation studies to better understand the performance of RACT. All simulations, unless otherwise noted, were run in the setting of $n_{1}=n_{2}=25,\;p=250$. For simulations related to the null hypothesis, the covariance matrix was $\Sigma_{1}$ for all observations; otherwise $\Sigma_{1}$ and $\Sigma_{2}$ represent the covariance matrices for group 1 and group 2. We simulated data from covariance matrices with changing low-rank, and low-rank block structures to reflect the block structures commonly found in biomedical data. In the data simulation settings S1-S3 below, we let $\Sigma_{1}=I+UU^{\top}$ and $\Sigma_{2}=I+VV^{\top}$ where $U$ and $V$ are low-rank matrices, so that $\Sigma_{1}-\Sigma_{2}$ is low-rank. We define $w_{U},w_{V}$ as the ranks of $U$ and $V$, and in our simulations either $w_{U}=w_{V}=2$ or $w_{U}=w_{V}=5$. To randomly generate appropriately sized blocks of the low-rank components of the covariance matrix $U_{1}U_{1}^{\top},U_{2}U_{2}^{\top},V_{1}V_{1}^{\top}$, we generated the columns of $U_{1},U_{2},V_{1}$ independently using the first $w$ singular vectors from randomly generated matrices $A_{1}A_{1}^{\top},A_{2}A_{2}^{\top},A_{3}A_{3}^{\top}$, where $A_{1},A_{2},A_{3}$ are appropriately sized matrices with i.i.d. standard normal entries. Figure 1 is provided for graphical illustrations of our simulation settings.

1. S1.

   (LowRank): We set $UU^{\top}=\tau^{2}U_{1}U_{1}^{\top},VV^{\top}=\tau^{2}V_{1}V_{1}^{\top}$ with $U_{1}U_{1}^{\top},V_{1}V_{1}^{\top}\in\mathbb{R}^{p\times p}$. Here, all elements of the covariance matrix experience a change with high probability.

2. S2.

   (LowRankBlockLarge): We set

   $\displaystyle UU^{\top}=\begin{bmatrix}\tau^{2}U_{1}U_{1}^{\top}&0\\ 0&U_{2}U_{2}^{\top}\end{bmatrix}\;\;\;VV^{\top}=\begin{bmatrix}\tau^{2}V_{1}V_{1}^{\top}&0\\ 0&U_{2}U_{2}^{\top}\end{bmatrix}\;\;\;U_{1}U_{1}^{T},V_{1}V_{1}^{T},U_{2}U_{2}^{\top}\in\mathbb{R}^{p/2\times p/2}.$

3. S3.

   (LowRankBlockSmall): A similar setup to S2 except $U_{1}U_{1}^{\top},V_{1}V_{1}^{\top}\in\mathbb{R}^{10\times 10}$ and $U_{2}U_{2}^{\top}\in\mathbb{R}^{(p-10)\times(p-10)}$.

4. S4.

   (OffDiagonal): Here $\Sigma_{1}=\begin{bmatrix}A_{1}&0\\ 0&I\end{bmatrix},\Sigma_{2}=\begin{bmatrix}A_{2}&0\\ 0&I\end{bmatrix}$ where $A_{1},A_{2},I\in\mathbb{R}^{p/2\times p/2}$. $A_{1}$ has an equicorrelation structure in that $a_{r,r}=1$ for $r\in 1,\dots,p/2$ and $a_{r,s}=\tau^{2}$ for $r\neq s$. $A_{2}$ is equal to $A_{1}$ except the covariances between dimensions $(1,\dots,\lceil p/4\rceil-1)$ and $(\lceil p/4\rceil,\dots,p/2)$ are set to $-\tau^{2}$.

To evaluate Type I error, we generated 10000 datasets and used 2000 permutations for each dataset. For power analysis we used 1000 datasets and 1000 permutations.

The Cancer Genome Atlas dataset (TCGA) contains 11,000 samples from 33 tumor types for a range of data modalities and has enabled a better understanding of the molecular changes which cancer causes (Hutter \BBA Zenklusen, \APACyear2018). Of the data modalities available from TCGA, we focus on gene expression, which serves as a vital link in understanding the mechanism through which DNA mutations can lead to cancer. The study of gene expression, recently aided by the development of new sequencing technologies such as RNA-Seq (Wang \BOthers., \APACyear2009), has led to a better understanding of the differences between cancer and normal cells (and therefore serves as a potential biomarker), and has helped in the identification of new disease subtypes (Gerstung \BOthers., \APACyear2015). Specifically, we analyze gene expression networks in two types of lung cancers, lung squamous cell carcinomas (LUSC) which is a common type of lung cancer (The Cancer Genome Atlas Research Network, \APACyear2012), and lung adenocarcinomas (LUAD) which is a leading cause of cancer death (The Cancer Genome Atlas Research Network, \APACyear2014). Gene expressions for protein coding genes for these tumors was downloaded using the BiocManager R package (Morgan \BBA Ramos, \APACyear2024), for which $n_{1}=553$ LUSC tumor samples and $n_{2}=600$ LUAD tumor samples were available from 19962 genes. We transform the data by taking $\log_{2}(1+\text{count})$ and keep the $p=500$ genes which exhibit the most variability across both tumor types. Finally within each tumor type, we regressed out age and sex from the data, and then used the residuals for our evaluations.

Using all samples we see in the difference of covariance matrices between tumor types some evidence of a low-rank structure. Using all samples, the first singular value and the first 93 singular values, represent 16% and 80% respectively of the total sum of all singular values. This suggests that when testing the equality of these covariance matrices much of the signal can be found in a limited number of low-rank structures. This is reflected in Figure 4 where we examine the empirical power of various Ky-Fan($k$) norms for a fixed subsample size of $n_{1}=n_{2}=15$. A rapid increase in power when $\alpha=0.05$ is seen as $k$ increases, before peaking at $k=13$ and decreasing modestly thereafter. The adaptivity of RACT is demonstrated in that its empirical power, which takes into account multiple $k$, settles near the maximum power across all $k$.

For comparing RACT to other methods, we first note that for moderate subsample sizes of between 15 and 50 we see LC and HC exhibit slightly inflated Type I error of approximately 10%, while CLX and DHW experience more severe Type I error inflation. Hence we use a permutation-based version of LC, HC, and CLX, and exclude DHW from our power comparison. RACT’s Type I error is controlled close to $\alpha=0.05$, while SY is slightly conservative. Figure 4 shows RACT exhibits increased power compared to SY, LC, CLX, and is comparable to that of HC.

The second dataset we analyze is from the Social Processes Initiative in Neurobiology of the Schizophrenia(s) (SPINS) group (Hawco \BOthers., \APACyear2021). This dataset consists of diffusion tensor imaging (DTI) measurements of fractional anisotropy (FA) and mean diffusivity (MD). DTI is a non-invasive technique which measures the diffusion of water molecules within neural tissue. MD relates to the amount of diffusion which is occurring, while FA is a summary statistic related to the architecture of the tissue being examined which is influenced by how uniform diffusion is for different directions (O’Donnell \BBA Westin, \APACyear2011; Van Hecke \BOthers., \APACyear2016). DTI has been found to be useful in a clinical setting in providing biomarkers for amyotrophic lateral sclerosis, and as a diagnostic tool for Parkinson’s disease (Tae \BOthers., \APACyear2018). Zhang \BOthers. (\APACyear2023) investigated the biases that are introduced for DTI imaging in the SPINS study when measurements are taken at different sites or when using different scanners, specifically inter-scanner covariance heterogeneity. In the SPINS study most sites began with General Electric 3T (GE) scanners, and all ended with Siemens Prisma 3T (SP) scanners. In the below analysis we provide further evidence of inter-scanner covariance heterogeneity in the SPINS study. Along the lines of Zhang \BOthers. (\APACyear2023) we use linear regression to remove the effect of age, age², gender, diagnosis, age $\times$ gender, and age $\times$ diagnosis, and then use these residuals to test for differences in the covariance matrix of observations for each scanner type. For both FA and MD each observation has $p=73$, and in total we have $n_{1}=130$ observations from GE scanners and $n_{2}=195$ observations from SP scanners.

An examination of the difference in covariance matrices using all samples reveals a low-rank structure for both FA and MD measurements. For FA the first singular value and the sum of the first 19 singular values represent 30% and 80% of the total sum of all singular values. The low-rank structure is more pronounced for MD where the first singular value represents 59% of the total sum, and the sum of the top 6 singular values represents 81% of the total sum. Using a fixed subsample size of 30 and 15 for FA and MD respectively we see in Figure 4 for individual Ky-Fan($k$) norms power is maximized at $k=12$ for FA and $k=14$ for MD. Reflecting the lower-rank structure of the MD data, we see a less steep increase in power as $k$ increases relative to FA. Similar to the TCGA data we see across both FA and MD datasets inflated Type I errors for LC, CLX, HC, and DHW; therefore we use permutation-based versions of LC, CLX, and HC and exclude DHW from the power analysis. RACT’s Type I error control is accurate, and SY is very conservative (Type I error of 0%). Relative to other methods we see similar results to what was seen in the TCGA dataset, with HC once again performing most comparably to RACT.

Across TCGA, SPINS FA, and SPINS MD we note that RACT’s best relative performance appears for the SPINS MD dataset, where its empirical power is strictly higher than other methods for subsample sizes up to 35. This potentially reflects the relatively lower-rank structure of the SPINS MD data, and correspondingly RACT’s ability to leverage this structure better than other methods. Also it is notable that RACT’s performance is consistently strong across all datasets whereas the performance of SY, CLX, and LC is more variable; this would be expected given RACT’s adaptive test statistic as compared to the test statistics of other methods which will be more sensitive to the specific form the difference in covariance takes. The method performing most similarly to RACT is HC. However, HC’s superdiagonal-based test statistic leads its performance to vary based on the ordering of the dimensions, and for an arbitrary reordering of the dimensions for these datasets will produce different power results. This could particularly affect HC’s results on SPINS FA, which we observe in Figure 3 to have several blocks along the diagonal of the difference in covariance matrices. This contrasts with RACT’s test statistic which is invariant to permutation of the dimensions.