# Large-scale correlation screening under dependence for brain functional connectivity network inference

Hanâ Lbath, Alexander Petersen, Sophie Achard

► To cite this version:

Hanâ Lbath, Alexander Petersen, Sophie Achard. Large-scale correlation screening under dependence for brain functional connectivity network inference. Statistics and Computing, 2024, 34 (2), pp.90. 10.1007/s11222-024-10411-x. hal-04884858

HAL Id: hal-04884858

https://hal.science/hal-04884858v1

Submitted on 13 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Large-Scale Correlation Screening under
Dependence for Brain Functional Connectivity
Network Inference

Hanâ Lbath¹*, Alexander Petersen ² and Sophie Achard ¹

¹*Univ. Grenoble Alpes, CNRS, Inria, Grenoble INP, LJK, 38000 Grenoble, France.

²Department of Statistics, Brigham Young University, Provo, Utah, 84602, USA.

Contributing authors: hana.lbath@inria.fr; petersen@stat.byu.edu; sophie.achard@univ-grenoble-alpes.fr;

Abstract

Data produced by resting-state functional Magnetic Resonance Imaging are widely used to infer brain functional connectivity networks. Such networks correlate neural signals to connect brain regions, which consist in groups of dependent voxels. Previous work has focused on aggregating data across voxels within predefined regions. However, the presence of within-region correlations has noticeable impacts on inter-regional correlation detection, and thus edge identification. To alleviate them, we propose to leverage techniques from the large-scale correlation screening literature, and derive simple and practical characterizations of the mean number of correlation discoveries that flexibly incorporate intra-regional dependence structures. A connectivity network inference framework is then presented. First, inter-regional correlation distributions are estimated. Then, correlation thresholds that can be tailored to one's application are constructed for each edge. Finally, the proposed framework is implemented on synthetic and real-world datasets. This novel approach for handling arbitrary intra-regional correlation is shown to limit false positives while improving true positive rates.

**Keywords:** brain functional connectivity, correlation screening, correlation threshold, network inference, rs-fMRI
# 1 Introduction

Large-scale network inference is a problem inherent to numerous fields, including gene regulatory networks, spatial data studies, and brain imaging. This work is motivated by an application to resting-state brain functional connectivity networks of single subjects. Such networks connect together correlated brain regions, which consist in groups of dependent voxels. These networks are key to providing insights into the diseased or injured brain (Achard et al, 2012; Richiardi et al, 2013; Malagurski et al, 2019). In this paper, the terms region and group will be used interchangeably, the former being associated with the motivating application, and the second with other data sources of similar structure to which the proposed methods also apply.

The goal of this work is to infer a binary network where nodes correspond to regions and edges are present only between nodes that are sufficiently highly correlated. The challenge is two-fold: not only does dependence between voxels within a region impact inter-regional correlation estimation, but it also affects inter-regional correlation threshold estimation, and thus edge detection. We propose a correlation screening approach, and tackle the problem of reliable large-scale correlation discovery between two groups of arbitrarily dependent variables.

In the context of brain functional connectivity, networks are often constructed from functional Magnetic Resonance Imaging (fMRI) data by spatially aggregating blood-oxygen-level-dependent (BOLD) time series within predefined brain regions, e.g., (Fallani et al, 2014). However, this may lead to overestimation of the inter-regional correlation, or inter-correlation for brevity (e.g., Halliwell (1962)), and hence incorrect edge detection. We propose a novel network inference framework that leverages, for each pair of regions, inter-correlation distributions instead of aggregation. To obtain the associated binary network we then present a thresholding step based on correlation screening. Existing approaches typically assume variables are independent within their region. Yet, as detailed in this work, any violation of this assumption markedly impacts
inter-correlation discovery, i.e., when the sample inter-correlation coefficient is greater
than a given threshold. As will be showcased later, high intra-regional correlation, or
intra-correlation for brevity, which corresponds to settings with homogeneous regions,
leads to lowered true positive rates (TPR). On the other hand, low intra-correlation, a
characteristic of inhomogeneous regions, leads to increased false positive rates (FPR).
In (Hero and Rajaratnam, 2011), a theoretical framework that accounts for arbitrary
dependence is presented. Their approach is nevertheless very difficult to implement in
practice and their empirical evaluation only covers the cases of independence or sparse
dependence. We hence introduce simple and practical expressions to characterize the
number of discoveries that flexibly incorporate dependence structures. These can then
be employed to find a correlation threshold per pair of regions that improves true
discovery rates under dependence, while limiting the number of false discoveries. The
main steps of the proposed pipeline are depicted in Figure 1 and are presented in
Sections 3, 4 and 5. We illustrate our work on synthetic data throughout this paper
and demonstrate the effectiveness of our framework on synthetic and real-world brain
rat imaging datasets in Section 6.

<figure><img src="image_2.png" /><figcaption>**Fig. 1:** Main steps of our proposed network inference pipeline. Each circle corresponds to a group of variables (represented by crosses). The sample inter-correlation estimation and thresholding steps are detailed for a pair of regions. Some edges were left out to improve readability.</figcaption></figure>
# 2 Preliminaries

In this section, we define the data model and its parameters that will be used in the rest of this paper.

## 2.1 Correlation coefficients

Let $a$ and $b$ be indices of two regions or groups consisting of $p_a$ and $p_b$ random variables, respectively. Denote by $\mathcal{R}_a$ the set of variables in region $a$ and $X_i^a$ the $i$th random variable in $\mathcal{R}_a$. Assume $n$ independent samples of $X_i^a$ are available and define the corresponding vector $\mathbf{X}_i^a = [X_{i,1}^a, \dots, X_{i,n}^a]^T$. $\mathcal{R}_b$, $X_j^b$ and $\mathbf{X}_j^b$ are similarly defined. As an illustration, in the context of brain functional connectivity, $X_i^a$ corresponds to voxel $i$ of brain region $a$, which is associated with an fMRI BOLD signal time series with $n$ time points. We define intra-correlation as the Pearson correlation between each pair of random variables within a given region. Inter-correlation is the Pearson correlation between pairs of random variables from two different regions.

Let $\rho_{ij}^{a,b}$ denote the true population inter-correlation coefficient between $X_i^a$ and $X_j^b$. For $a \neq b$, define the corresponding sample inter-correlation coefficient

$$R_{i,j}^{a,b} = \frac{\sum_{k=1}^{n} (X_{i,k}^{a} - \overline{X}_{i}^{a})(X_{j,k}^{b} - \overline{X}_{j}^{b})}{\sqrt{\sum_{k=1}^{n} (X_{i,k}^{a} - \overline{X}_{i}^{a})^{2} \sum_{k=1}^{n} (X_{j,k}^{b} - \overline{X}_{j}^{b})^{2}}}, \quad (1)$$

with $\overline{X}_i^a$, $\overline{X}_j^b$ the sample means. Denote the probability density, cumulative distribution, and quantile functions of $R_{i,j}^{a,b}$ by $f_{R_{i,j}^{a,b}}$, $F_{R_{i,j}^{a,b}}$, and $F_{R_{i,j}^{a,b}}^{-1}$, respectively. Population and sample intra-correlation coefficients and their distributions can analogously be defined by choosing $b = a$. Asymptotic closed-form expressions of the density of correlation can be obtained for Gaussian independently identically distributed (i.i.d.) variables $X_i^a$, $X_j^b$ (Muirhead, 2005). Note however that this work aims to tackle arbitrary dependence between variables, and in this context, to the best of our knowledge,
such explicit formulas have not been derived without defining a parametric dependence
structure.

In most of this paper, and for ease of calculation, we assume the joint distribution
of pairs of voxels *i*, *j* from a fixed region pair *a*, *b* are identically distributed. In such
cases, the sample inter-correlation coefficients $R_{i,j}^{a,b}$ are identically distributed, and the
*i*, *j* subscripts will be dropped, though we emphasize that independence within regions
is not assumed.

## 2.2 Synthetic data examples

We illustrate the different concepts introduced in this paper with data simulated as
follows. We consider two regions *a* and *b*, both containing *p* intra-correlated variables
following a multivariate normal distribution with a predefined Toeplitz covariance
structure. *n* independent samples of each of these *p* variables are generated. We hence
obtain data with a block diagonal covariance matrix of size 2*p* × 2*p*, where each block
corresponds to each region. The off-diagonal blocks correspond to the inter-correlation
coefficients, which are set to be constant across all pairs of voxels. The diagonal blocks
correspond to the intra-correlation coefficients, which follow a Toeplitz dependence
structure.

# 3 Inter-correlation estimation

## 3.1 Related work

Previous works on the estimation of inter-correlations have mostly focused on aggregating variables within predefined regions (Fallani et al, 2014; Dadi et al, 2019). In the context of brain functional connectivity network inference, some prefer techniques based on independent component analysis (ICA) (Calhoun et al, 2012), while most focus on summarizing all voxels within predefined brain regions by their average, e.g., (Achard et al, 2012, 2006; Di Martino et al, 2014; Malagurski et al, 2019). However,
such approaches suffer from loss of relevant information and can lead to statistical inconsistency and incorrect correlation estimation (Ostroff, 1993). In particular, the estimate of the average of weakly correlated time series, which corresponds to samples of a single variable in our data model, is poor (Wigley et al, 1984). Additionally, it has been observed on small samples that the correlation of averages is different than the average of correlations (Dunlap et al, 1983). This phenomenon can also be easily checked with arbitrary large samples. Furthermore, correlation of averages were empirically observed to overestimate the true correlation (Halliwell, 1962; Achard et al, 2011). Therefore, when correlating regional averages for binary network inference, one will tend to identify spurious edges. Some previous works attempted to improve false positive rate control utilizing multiple testing approaches (Drton and Perlman, 2007). However, in the context of arbitrary dependence structures, such methods cannot be straightforwardly applied. One alternative to aggregation is to measure the similarity, such as the Wasserstein distance or covariance (Petersen and Müller, 2019) between intra-correlation densities. However, this approach is not equivalent to that of the Pearson correlation. Indeed, while the Wasserstein distance may provide a first intuition about how regions are connected, it does not capture as much information about the relationship between the two regions as inter-correlations do.

## 3.2 On the impact of intra-correlation on inter-correlation estimation and detection

We first illustrate how intra-correlation affects the sample inter-correlation distribution in a simplified scenario, before considering a more general case. It has been known for some time in familial data studies that intra-correlations impact inter-correlation estimation (Rosner et al, 1977; Donner and Eliasziw, 1991). In the multivariate normal case, and under the assumption of within-group homoscedasticity, the asymptotic variance of the maximum-likelihood estimator of the inter-correlation, denoted as $R_{MLE}^{a,b}$,
was derived by Elston (1975). This estimator showcases similar behavior to the voxel-to-voxel sample inter-correlation coefficients $R_{i,j}^{a,b}$ and will help provide us with a first intuition about the impact of intra-correlation. We need to assume all variables $X_i^a$, $X_j^b$ have the same true inter-correlation $\rho^{a,b}$ and intra-correlation $\rho^{a,a}$ and $\rho^{b,b}$. This amounts to saying sample intra- and inter-correlation coefficients are identically distributed within their corresponding group, or pair of groups, respectively. Under these assumptions, and according to Elston (1975), the variance of the maximum-likelihood estimator is:

$$ \mathrm{Var}(R_{\mathrm{MLE}}^{a,b}) = \frac{1}{n} \left[ (\rho^{a,b})^2 - \frac{1}{p_a} [1 + (p_a - 1)\rho^{a,a}] \right] \times \left[ (\rho^{a,b})^2 - \frac{1}{p_b} [1 + (p_b - 1)\rho^{b,b}] \right] \\ + \frac{(\rho^{a,b})^2}{2n} \left[ \frac{p_a - 1}{p_a} (1 - \rho^{a,a})^2 + \frac{p_b - 1}{p_b} (1 - \rho^{b,b})^2 \right] \quad (2) $$

The expression in (2) shows that the variance of the sample inter-correlation coefficient explicitly depends on the true intra-correlation coefficients $\rho^{a,a}$ and $\rho^{b,b}$. When the number of samples $n$ is sufficiently large, the inter-correlation variance in the multivariate normal case hence increases when intra-correlation decreases. This observation implies that for a fixed threshold that does not depend on regional dependency structures, more false positive correlations are likely to be discovered.

This intuition is illustrated in the left hand-side of Figure 2 where the true inter-correlation is zero and no positive correlations are expected to be discovered. Conversely, for the same fixed threshold, when the true inter-correlation is positive (cf. right hand-side of Figure 2), increased intra-correlations, which lead to lower inter-correlation variance, may lead to decreased number of true positives. This phenomenon is observed regardless of the number of time points or variables (cf. supplementary materials).
<figure><img src="image_3.png" /><figcaption>**Fig. 2:** Effect of intra-correlation on sample inter-correlation distribution, for different population inter-correlation values. The correlation samples were computed between all pairs of variables from two groups. Each group contains $n = 150$ samples of $p = 500$ intra-correlated random variables following a multivariate normal distribution with Toeplitz intra-correlation (cf. Section 2.2). We can note the higher the intra-correlation, the lower the variance of the inter-correlation distribution.</figcaption></figure>

In fact, the impact that intra-correlation distributions have on the spatial average of sample inter-correlations can be quantified even without any distributional assumptions. The following result shows that, when the intra-correlation densities of two regions a and b are highly dissimilar, as quantified by a large Wasserstein distance, the average inter-correlation is upper-bounded. In particular, this phenomenon is by no means limited to the Gaussian case. We recall here the definition of the Wasserstein distance between two correlation densities (Petersen and Müller, 2019; Panaretos and Zemel, 2019): $d_W^2(f_{R^{a,a}}, f_{R^{b,b}}) = \int_0^1 [F_{R^{a,a}}^{-1}(c) - F_{R^{b,b}}^{-1}(c)]^2 dc$. The full proof is available in Appendix A.

**Proposition 1.** For any regions a, b with $p_a, p_b$ voxels, respectively, and $\mathcal{R}_a, \mathcal{R}_b$ the corresponding voxel sets, if there exists $A \in \mathbb{R}^+$ such that $d_W^2(f_{R^{a,a}}, f_{R^{b,b}}) \ge \min_{c \in [0,1]} (F_{R^{a,a}}^{-1}(c) - F_{R^{b,b}}^{-1}(c))^2 \ge A$, then, $\overline{R^{a,b}} = \frac{1}{p_a p_b} \sum_{i \in \mathcal{R}_a} \sum_{j \in \mathcal{R}_b} R_{i,j}^{a,b} \le 1 - \frac{\sqrt{A}}{2}$.

## 3.3 Proposed approach: inter-correlation distribution estimation

As previously discussed, aggregating variables within regions to estimate inter-correlation leads to loss of information and incorrect edge detection during the binary
network inference step. We have also brought to light the importance of taking intra-
correlation into account when manipulating inter-correlations. In addition, all the
previously cited approaches that aim to infer a binary network where nodes are groups
of variables only provide a single correlation threshold to be applied to all pairs of
regions. In this paper, we propose to derive a correlation threshold specific to each
pair of regions to better harness the particularities of the regional dependence struc-
tures. Instead of averaging variables within regions, we hence propose to estimate the
distribution of correlations measured between all pairs of variables from two differ-
ent regions. We then obtain an inter-correlation distribution per pair of region, which
then needs to be thresholded. To that end, we propose to leverage correlation screen-
ing. In that paradigm, an edge is said to be detected in the associated binary graph
if a sufficient number of sample inter-correlation coefficients of the corresponding pair
of regions are large enough. In the following sections, we derive simplified expres-
sions of the number of discoveries to propose a reliable method to threshold these
inter-correlation distributions.

4 Characterization of the number of discoveries
under dependence

Correlation screening (Hero and Rajaratnam, 2011), or independence screening (Fan and Lv, 2008), is often used in variable or feature selection problems. In such approaches, the goal is to discover sufficiently highly correlated variables. A practical method consists in defining a correlation threshold, above which correlation coefficients, and their associated variables, are said to be *detected* or *discovered*. Nonetheless, in high dimension, such approaches may suffer from a high number of false discoveries. In (Hero and Rajaratnam, 2011), the authors aim to mitigate this issue in the following way. They first propose the following maximum-based definition of the number of discoveries, pertaining to inter-correlation coefficients, with $\phi_{ij}^{ab}(\rho) = 1(|R_{i,j}^{a,b}| > \rho)$ for
all voxels i, j in regions a, b and correlation threshold $\rho \in [0, 1]$:

$$N^{ab}(\rho) = \sum_{i=1}^{p_a} \max_{j=1, \dots, p_b} \phi_{ij}^{ab}(\rho). \qquad (3)$$

The authors provide as well an approximation of the expected number of discoveries $E[N^{ab}]$ that depends on the number of variables $p$, the number of samples $n$ and a function of the joint distribution of a transformation of the variables. Then they employ the derived formula to compute critical threshold values based on a phase transition approach. Furthermore, the expected number of discoveries is used to control the number of false discoveries. It is then all the more essential to have an expression of the number of discoveries that is both interpretable and can easily be theoretically and empirically utilized. However, the expression for $E[N^{ab}]$ derived in (Hero and Rajaratnam, 2011), which depends on joint distributions, is difficult to compute, especially for single subject analysis where we have access to only a single sample of each signal. We hence provide simplified explicit expressions of the mean number of discoveries that still harness information contained in the intra-correlation distributions.

## 4.1 Maximum-based expression: $N^{ab}$

Empirically, intra-correlation has an impact on $N^{ab}$ and its average (cf. Figure 3, left). Indeed, for a given inter-correlation threshold, the smaller the intra-correlations, the larger the number of discoveries. This is in accordance with the observations from Figure 2. Additionally, we can remark that in this example the true inter-correlation is zero. Thus any discovery is a false positive. As the intra-correlation increases, a lower correlation threshold is then sufficient to maintain similar levels of false discoveries.

From (3), we can also conclude that for any $\rho \in [0, 1]$,

$$E[N^{ab}(\rho)] = \sum_{i=1}^{p_a} \left(1 - F_{|R_{i,1}^{a,b}|, \dots, |R_{i,p_b}^{a,b}|}(\rho, \dots, \rho)\right), \qquad (4)$$
<figure><img src="image_4.png" /><figcaption>Fig. 3: Normalized number of discoveries, $N^{ab}/p_a$ (Left) and $N_e^{ab}/p_a p_b$ (Right) as a function of the inter-correlation threshold for data simulated as described in Section 2.2, with $p_a = p_b = 500, n = 150$, true inter-correlation $\rho^{ab} = 0$ and Toeplitz intra-correlation with varying minimal intra-correlations. For each of the four intra-correlation values, 50 datasets were simulated and used to compute the number of discoveries (the colored curves), and their average (the black dotted curves). $N^{ab}/p_a$ decreases as the intra-correlation increases, while $N_e^{ab}/p_a p_b$ does not seem to be much impacted on average, even though its variability seems to increase with the intra-correlation value.</figcaption></figure>

with $F_{|R_{i,1}^a|, \dots, |R_{i,p_b}^a|}$ the joint distribution of the absolute values of the corresponding correlation coefficients, which will inherently take into account dependence structures between the inter-correlation coefficients. However, joint distributions are complicated to estimate and manipulate. Let $\tilde{\nu}^{ab} = E[N^{ab}]/p_a$. We then propose an approximate expression of $\tilde{\nu}^{ab}$, denoted by $\nu^{ab}$, that depends on the distribution of inter-correlations, and that is exact under some particular assumptions (cf. Appendix C):

$$ \nu^{ab}(\rho) = 1 - F_{|R^{a,b}|}(\rho)^{p_b}, \quad \rho \in [0, 1]. \qquad (5) $$

We can also derive the following inequality.

**Proposition 2.** Consider two regions a and b and a correlation threshold $\rho \in [0, 1]$. If all variables $X_i^a$ and $X_j^b$ in both regions follow a normal distribution and their sample inter-correlation coefficients are identically distributed, then for sufficiently large n,

$$ \nu^{ab}(\rho) \ge \tilde{\nu}^{ab}(\rho) \qquad (6) $$
Proof. As defined in (Lehmann, 1966), the random variables $T_1, T_2, \dots, T_p$ are Positively Quadrant Dependent (PQD) if for any positive number $t_1, t_2, \dots, t_p$,

$$P\left(\bigcap_{k=1}^{p} T_k \le t_k\right) \ge \prod_{k=1}^{p} P(T_k \le t_k). \qquad (7)$$

Under the assumption $X_i^a, X_j^b$ are normal for all $i, j$, the distribution of their sample correlation coefficients $R_{i,j}^{a,b}$ is asymptotically normal, e.g., (Ruben, 1966; Hotelling, 1953). Hence, according to Theorem 1 in (Šidák, 1967), when $n$ is large enough, $|R_{i,j}^{a,b}|$ are PQD. We can then note that equation (7) is equivalent to $E[\prod_{k=1}^{p} 1(|T_k| \le t_k)] \ge \prod_{k=1}^{p} E[1(|T_k| \le t_k)]$. Under the assumption the sample correlation coefficients are identically distributed, and setting $t_k = \rho$ and $T_k = R^{a,b}$ for all variables, we can thus write:

$$p_a \cdot \nu^{ab}(\rho) = \sum_{i=1}^{p_a} \left( 1 - \prod_{j=1}^{p_b} E[1(|R^{a,b}| \le \rho)] \right) \\ \ge \sum_{i=1}^{p_a} \left( 1 - E \left[ \prod_{j=1}^{p_b} 1(|R^{a,b}| \le \rho) \right] \right) = E[N^{ab}(\rho)] = p_a \cdot \tilde{\nu}^{ab}(\rho).$$

This concludes the proof. □

This result ensures $\nu^{ab}$ will provide thresholds that are at least as conservative as that of $\tilde{\nu}^{ab}$. Moreover, the assumptions needed in Proposition 2 are often reasonable in practice—and notably in functional brain connectivity applications where the signals associated with each voxel can appropriately be transformed (Whitcher et al, 2000). In addition, $\nu^{ab}$ can be estimated by $\hat{\nu}^{ab}(\rho) = 1 - \hat{F}_{|R^{a,b}|}(\rho)^{p_b}$, where $\hat{F}_{|R^{a,b}|}(\rho) = \frac{1}{p_b p_a} \sum_{i=1}^{p_a} \sum_{j=1}^{p_b} 1(|R_{i,j}^{a,b}| \le \rho)$ is the empirical cumulative distribution function (ecdf) of $|R^{a,b}|$. The result stated above can then be empirically observed in Figure 4. The curve of $\hat{\nu}^{ab}$ as a function of thresholds is also particularly close to that of the empirical values of $\tilde{\nu}^{ab}$ as long as both the inter-correlation and the
intra-correlation of region b are not too high. $\hat{\nu}^{ab}$ provides hence an approximation for the normalized expected number of discoveries that is easier to compute, while still accounting for the inter-correlation distribution. Moreover, in practice, the ecdf of the inter-correlation coefficients can be shown to depend on the intra-correlation structure (Azriel and Schwartzman, 2014) and hence allows us to account for it.

<figure><img src="image_5.png" /><figcaption>Fig. 4: Normalized number of discoveries $\hat{\nu}^{ab}$, $\hat{\nu}_e^{ab}$, $N^{ab}/p_a$ and $N_e^{ab}/p_a p_b$ as a function of the correlation threshold for different minimal intra-correlation values $\rho_{min}^{aa}$, $\rho_{min}^{bb}$ and a zero inter-correlation for data simulated as described in Section 2.2, with $p_a = p_b = 500$, $n = 150$. We can note $\hat{\nu}_e^{ab} \le \hat{\nu}^{ab}$. We can also observe the asymmetric behavior of $N^{ab}/p_a$ and how it is close to $\hat{\nu}^{ab}$ when $\rho_{min}^{aa} \gg \rho_{min}^{bb}$, unlike when $\rho_{min}^{bb} \gg \rho_{min}^{aa}$.</figcaption></figure>

## 4.2 Sum-based expression: $N_e^{ab}$

We now present another, and more intuitive, way to characterize the number of discoveries $N_e^{ab}$ (Hero and Rajaratnam, 2011). It represents the total number of discoveries and will enable us to propose less conservative thresholds:

$$N_e^{ab}(\rho) = \sum_{i=1}^{p_a} \sum_{j=1}^{p_b} \phi_{ij}^{ab}(\rho), \quad \rho \in [0, 1]. \qquad (8)$$
However, it is not straightforward to derive a critical correlation threshold from this expression. We propose the simplified expression below:

$$ \hat{\nu}_e^{ab}(\rho) = 1 - \hat{F}_{|R^{a,b}|}(\rho), \quad \rho \in [0, 1]. \qquad (9) $$

We can also remark in Figure 4 that $\hat{\nu}_e^{ab}$ and $N_e^{ab}/p_a p_b$ look indistinguishable.

## 4.3 Link between $N_e^{ab}$ and $N^{ab}$

We have presented so far two ways to characterize the number of discoveries. We will now discuss how they relate to one another. We can remark the following inequality.

**Proposition 3.** For all $\rho \in [0, 1]$,

$$ \hat{\nu}_e^{ab}(\rho) \le \hat{\nu}^{ab}(\rho). \qquad (10) $$

The proof is straightforward and can be found in Appendix B. This result can notably be observed in Figure 4. Thus $\hat{\nu}^{ab}$ is more conservative than $\hat{\nu}_e^{ab}$, which in some circumstances may be desirable. Nonetheless, when the inter-correlation is zero we expect no discoveries. In this case, a critical correlation threshold can hence be defined as the minimum correlation such that the number of discoveries is zero. In such cases, using $\hat{\nu}_e^{ab}$ then seems to be preferable, since it provides a lower correlation threshold for a similar number of false discoveries.

# 5 Correlation threshold definition

Now we have better characterized the number of discoveries, we can use it to construct correlation thresholds tailored to one's data and that ensure, to a certain extent, a restricted number of false discoveries and improved number of true discoveries. We
present in this section two possible correlation threshold definition approaches. The
idea behind correlation threshold definition is to ensure that it is very unlikely for any
discovery to correspond to a correlation value that could have happened at random.
This amounts to a setting where the true inter-correlation is zero, which may not be
true in practice. Surrogate data defined such that the population inter-correlation is
zero can hence be utilized to estimate the correlation thresholds. We denote $\hat{F}_{0,|R^{a,b}|}^{-1}$
the corresponding quantile function.

**5.1 FWER-based threshold**

Correlation thresholds with family-wise error rate (FWER) theoretical control can be
derived for specific dependence structures. These approaches control the probability
of making at least one false discovery. Analogously to Proposition 2 in (Hero and
Rajaratnam, 2011), it can be shown, under a weak dependence condition, that $N_e^{ab}$
converges to a Poisson random variable when $p_a, p_b \to \infty$ and $P(N_e^{ab} > 0) \to 1 -$
$\exp(-E[N_e^{ab}])$. Our proposed expression $\hat{\nu}_e^{ab}$ can then be used to compute correlation
thresholds $\rho_\alpha^{ab}$ that guarantee a FWER at level $\alpha$. Nevertheless, the weak dependence
assumption upon which this approach hinges is often not reasonable in practice.

**5.2 Quantile-based threshold**

The correlation threshold can also be defined such that the False Positive Rate (FPR)
is guaranteed to be less than a given level $\alpha$. The FPR is the ratio between the number
of false positives (FP) and the total number of ground truth negatives, that is $p_a \cdot p_b$
in the $\rho^{ab} = 0$ case. Controlling the FPR at level $\alpha$ is thus equivalent to ensuring the
number of discoveries is less than $FP = \alpha \cdot p_a \cdot p_b$. Since in our setting ($\rho^{ab} = 0$) any
discovery is a false positive, we can set $\hat{\nu}_e^{ab} \cdot p_a \cdot p_b = \alpha \cdot p_a \cdot p_b$, which leads to the
threshold $\rho_{q,\alpha}^{ab} = \hat{F}_{0,|R^{a,b}|}^{-1}(1-\alpha)$. We can remark that when $\alpha = 0$, the chosen threshold
is larger than any of the observed absolute correlations, ensuring there will be no
discoveries. Additionally, this threshold will depend on the intra-correlation, as does
the ecdf (Azriel and Schwartzman, 2014). We can also remark this threshold guarantees a FWER at level $\alpha = 0$ under the previous weak dependence assumption. $\hat{\nu}^{ab}$ can also be used to similarly derive a critical threshold, although stringent conditions on the region sizes would then need to be verified when $\alpha \neq 0$.

<figure><img src="image_6.png" /><figcaption>Fig. 5: Comparison of different critical thresholds for 50 replicates of data simulated as described in Section 2.2 with $p = 150$, $n = 100$, $\rho^{a,b} = 0$ and varying constant intra-correlation values. The FWER- and quantile-based thresholds $\rho_{\alpha}^{ab}$ were computed for $\alpha = 0$.</figcaption></figure>

## 5.3 Numerical results

We compare in Figure 5 the two correlation thresholds defined above with two other approaches: the critical thresholds $\rho_{\text{hero}}$ obtained in (Hero and Rajaratnam, 2011), and a simple method where the threshold is set to $\rho_{\text{poli}}^{ab} = \hat{\mu}_{0}^{ab} + \hat{\sigma}_{0}^{ab}$ where $\hat{\mu}_{0}^{ab}$ and $\hat{\sigma}_{0}^{ab}$ are the sample mean and standard deviation of the sample inter-correlation of the surrogate data (Poli et al, 2015). We observe that, when intra-correlation is lower than 0.5, both our proposed approaches provide similar thresholds to Hero and Rajaratnam (2011). Nonetheless, our FWER- and quantile-based methods provide less conservative thresholds when the intra-correlation is high. The lower thresholds imply an increase in the true positive rate. We can note as well that $\rho_{\text{poli}}^{ab}$ is lower than all three other
thresholds for all intra-correlation values, which could lead to a large number of false discoveries, as will be shown in the next section. It is also decreasing when intra-correlation increases in accordance with the observations about the effect of intra-correlation on the distribution of sample inter-correlation in Section 3.

# 6 Network inference results

In this section we provide an illustration of our network inference approach on synthetic and real-world data and compare it to several methods.

## 6.1 Comparison to other methods

As mentioned in Section 3.1, most single-subject fMRI studies that use anatomical parcellations estimate the inter-correlation by computing the correlation coefficient between spatial regional averages of the signals. We will refer to the correlation of averages approach by CA and our proposed correlation screening method by CS. Various methods can then be employed to define the correlation thresholds. They usually belong to either of these two paradigms: (i) *relative thresholding*, that is, estimation of the binary network by extracting a fixed proportion of edges, e.g., (van den Heuvel et al, 2017), and (ii) *absolute thresholding* where one chooses, more or less arbitrarily, a fixed threshold that will be applied to all edges—as opposed to our proposed edge-specific thresholds. These two approaches are implemented in several popular packages, such as the Brain Connectivity toolbox (Rubinov and Sporns, 2010) and CONN (Whitfield-Gabrieli and Nieto-Castanon, 2012), with little guidance on the choice of threshold. Relative thresholds will always detect the same proportion of edges regardless of the true inter-correlation value, and as such can be disregarded in this work. Both the thresholds proposed by Hero and Rajaratnam (2011) and Poli et al (2015) are absolute thresholds. The latter was used in (Boschi et al, 2021) to threshold CA-based functional connectivity. In (Becq et al, 2020b), the authors apply
to all edges of a CA-based functional connectivity network a fixed threshold $\rho_{becq}$ determined according to a multiple testing approach (see Appendix G of their work for more details). In practice, the values of $\rho_{becq}$ are close to that of $\rho_{hero}$. The thresholds $\rho_{poli}^{ab}$, $\rho_{\alpha}^{ab}$ and $\rho_{q,\alpha}^{ab}$ are estimated using surrogate data where the true inter-correlation is zero and the intra-correlation is constant and equal to the average sample intra-correlation.

## 6.2 Synthetic data results

We generated synthetic datasets with ten inter-connected regions. For each dataset, 10 regions are simultaneously simulated, each region containing $p = 150$ intra-correlated variables following a multivariate normal distribution. 100 independent samples of each of these variables are obtained. For each region, a Toeplitz intra-correlation is used with the same minimal intra-correlation value across all ten regions. There are 41 true positive (constant true inter-correlation $\rho^{ab} = 0.2$) and 4 true negative edges (constant inter-correlation $\rho^{ab} = 0$) in the ground-truth network. The different simulation parameters are chosen to ensure the population covariance matrix of the ten regions is positive semidefinite. To identify the edges of the binary network, the pairwise thresholds are applied to the corresponding distributions of the absolute value of the sample inter-correlation. In particular, pairs of regions where the inter-correlation is larger than the threshold with a probability at most 0.05 are not identified as edges.

Table 1 displays the false positive and true positive rates (FPR and TPR, respectively) of the different methods, for varying minimal intra-correlation values. The FPR and TPR are defined as follows: FPR = FP/(FP+TN) and TPR = TP/(TP+FN). FN stands for false negatives (i.e., an edge is undetected when it actually exists). FPR is expected to be close to 0 and TPR to 1. Results in Table 1 showcase that, as expected from the previous section, using $\rho_{poli}^{ab}$ leads to high FPRs, while $\rho_{hero}$ and $\rho_{becq}$ lead to decreasing TPRs as intra-correlation increases. Additionally, while the FPR is slightly
<table><caption>**Table 1:** Comparison of the mean (standard deviation) TPR and FPR of several network inference methods on synthetic data across 100 repetitions. The ground-truth networks consist of 10 regions, with 4 real negative edges (true inter-correlation $\rho^{ab} = 0$), 41 real positive edges (true inter-correlation $\rho^{ab} = 0.2$), five minimal intra-correlation values $\rho^{aa}_{min}$, $n = 100$, and $p = 150$.</caption><thead><tr><th rowspan="2">Method</th><th colspan="2">ρ<sup>aa</sup><sub>min</sub> = 0.5</th><th colspan="2">ρ<sup>aa</sup><sub>min</sub> = 0.6</th><th colspan="2">ρ<sup>aa</sup><sub>min</sub> = 0.7</th><th colspan="2">ρ<sup>aa</sup><sub>min</sub> = 0.8</th><th colspan="2">ρ<sup>aa</sup><sub>min</sub> = 0.9</th></tr><tr><th>FPR</th><th>TPR</th><th>FPR</th><th>TPR</th><th>FPR</th><th>TPR</th><th>FPR</th><th>TPR</th><th>FPR</th><th>TPR</th></tr></thead><tbody><tr><td>CA + ρ<sup>ab</sup><sub>poli</sub></td><td>0.51 (0.23)</td><td>0.97 (0.02)</td><td>0.46 (0.22)</td><td>0.90 (0.07)</td><td>0.45 (0.25)</td><td>0.69 (0.12)</td><td>0.43 (0.26)</td><td>0.33 (0.14)</td><td>0.30 (0.24)</td><td>0.04 (0.03)</td></tr><tr><td>CA + ρ<sup>becq</sup></td><td>0.13 (0.13)</td><td>0.47 (0.20)</td><td>0.04 (0.09)</td><td>0.22 (0.13)</td><td>0.03 (0.08)</td><td>0.11 (0.09)</td><td>0.01 (0.04)</td><td>0.05 (0.06)</td><td>0.01 (0.04)</td><td>0.04 (0.04)</td></tr><tr><td>CS + ρ<sup>ab</sup><sub>poli</sub></td><td>1.00 (0.00)</td><td>1.00 (0.00)</td><td>1.00 (0.00)</td><td>1.00 (0.00)</td><td>1.00 (0.00)</td><td>1.00 (0.00)</td><td>1.00 (0.00)</td><td>1.00 (0.00)</td><td>1.00 (0.00)</td><td>1.00 (0.00)</td></tr><tr><td>CS + ρ<sup>hero</sup></td><td>0.03 (0.08)</td><td>0.13 (0.12)</td><td>0.03 (0.08)</td><td>0.13 (0.11)</td><td>0.04 (0.09)</td><td>0.12 (0.09)</td><td>0.02 (0.07)</td><td>0.11 (0.09)</td><td>0.01 (0.06)</td><td>0.09 (0.07)</td></tr><tr><td>CS + ρ<sup>ab</sup><sub>poli</sub></td><td>0.06 (0.10)</td><td>0.23 (0.17)</td><td>0.08 (0.12)</td><td>0.32 (0.17)</td><td>0.13 (0.13)</td><td>0.42 (0.17)</td><td>0.15 (0.16)</td><td>0.52 (0.16)</td><td>0.23 (0.17)</td><td>0.68 (0.13)</td></tr><tr><td>CS + ρ<sup>aa</sup><sub>poli</sub></td><td>0.06 (0.11)</td><td>0.23 (0.18)</td><td>0.07 (0.11)</td><td>0.32 (0.17)</td><td>0.14 (0.13)</td><td>0.42 (0.16)</td><td>0.15 (0.16)</td><td>0.53 (0.16)</td><td>0.25 (0.18)</td><td>0.67 (0.13)</td></tr></tbody></table>
increased, correlation screening methods with FWER- or quantile-based thresholds markedly improve the TPR when the intra-correlation is high, and should be preferred in that case. Indeed, when intra-correlation is 0.9 all other methods (except CS + $ρ_{poli}^{ab}$) have a TPR close to zero, while the FWER- and quantile-based thresholds have a TPR close to 0.7. The CS+$ρ_{poli}^{ab}$ method displays a FPR of 1 for all intra-correlation values and should thus to be avoided. Since the FWER- and quantile-based thresholds are empirically equivalent, from now on we will be using $ρ_{q,α=0}^{ab}$, which, unlike $ρ_{α=0}^{ab}$, has a theoretical control over false positives that is valid for any dependence structure.

## 6.3 Real-world data results

We applied our framework on functional Magnetic Resonance Imaging (fMRI) data acquired on both dead and live rats, anesthetized using Isoflurane (Becq et al, 2020a,b). The datasets are freely available at https://dx.doi.org/10.5281/zenodo.7254133. The scanning duration was 30 minutes with a time repetition of 0.5 second so that 3600 time points were acquired. After preprocessing as explained in (Becq et al, 2020b), based on an anatomical atlas, 51 groups of time series, corresponding to the rat brain regions, were extracted for each rat. Due to insufficient signal, four regions were excluded. Each time series captures the functioning of a given voxel. The dead rats provide experimental data where the ground-truth network is empty. Indeed, no legitimate functional activity should be detected, whereas for the live rat under anesthetic, we expect non-empty graphs as brain activity keeps on during anesthesia. We can note that no ground-truth is available for the live rat networks. As expected, the networks of the dead rats estimated using our proposed correlation screening method are empty, i.e. it does not detect any false positive edges, with the exception of one edge in one rat. However, the CA + $ρ_{poli}^{ab}$ approach (Boschi et al, 2021; Poli et al, 2015) detects over 300 false positive edges, and Becq et al (2020b) (later denoted B2020) detect between one and four false positive edges (cf. Table 2). While our approach is more
conservative than the other two, important edges are still detected in the live rats, mainly in motor regions (M1 and M2) and somatosensory regions (S1 and S2), as shown for instance in Figure 6.

<table><caption>**Table 2:** Comparison of the number of edges in the networks obtained via our proposed network inference approach and two methods from the literature (B2020 and CA + $\rho_{poli}^{ab}$) for dead (Top) and live (Bottom) rat brain fMRI data. In the dead rat brain networks, any detected edge is a false positive.</caption><thead><tr><th rowspan="2">DEAD RATS ID</th><th colspan="3">NUMBER OF EDGES</th></tr><tr><th>CS + ρ<sup>ab</sup><sub>q,α=0</sub></th><th>B2020</th><th>CA + ρ<sup>ab</sup><sub>poli</sub></th></tr></thead><tbody><tr><td>20160524_153000</td><td>1</td><td>4</td><td>316</td></tr><tr><td>20160609_161917</td><td>0</td><td>4</td><td>317</td></tr><tr><td>20160610_121044</td><td>0</td><td>1</td><td>325</td></tr></tbody></table>

<table><thead><tr><th rowspan="2">LIVE RATS ID</th><th colspan="3">NUMBER OF EDGES</th></tr><tr><th>CS + ρ<sup>ab</sup><sub>q,α=0</sub></th><th>B2020</th><th>CA + ρ<sup>ab</sup><sub>poli</sub></th></tr></thead><tbody><tr><td>20160615_103000</td><td>25</td><td>647</td><td>820</td></tr><tr><td>20160614_095825</td><td>411</td><td>847</td><td>910</td></tr><tr><td>20160615_121820</td><td>116</td><td>477</td><td>692</td></tr><tr><td>20160421_133725</td><td>83</td><td>591</td><td>910</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Fig. 6:** Brain functional connectivity network of a dead and two live rats (anesthetized with Isofluorane) inferred using our proposed correlation screening framework with the quantile-based threshold (CS + $\rho_{q,\alpha=0}^{a,b}$).</figcaption></figure>
# 7 Discussion

We have presented a novel approach to infer connectivity networks when nodes represent groups of correlated variables. We have formally established the importance of leveraging dependence structures to reliably discover inter-correlations. Our method consists in estimating, for each pair of groups, an inter-correlation distribution before deriving a tailored threshold based on a correlation screening approach. In particular, we proposed simplified expressions for the mean number of discoveries that allow for easier theoretical and empirical manipulation, and flexibly take into account dependence within groups. Motivated by a real-world application, we have demonstrated the feasibility of our approach on a real dataset of rat brain images.

This work has several possible theoretical extensions. First, while we provide a method with theoretical FPR control for any setting, we provide theoretical FWER control guarantees only under a weak dependence assumption. Nevertheless, this assumption is often unrealistic and relaxing it is difficult and would be an interesting direction to explore. Additionally, FWER approaches may sometimes be too conservative. On the other hand, the false discovery rate (FDR) enables to control the average number of FPs, which is often sufficient. A procedure to define correlation thresholds was proposed in (Cai and Liu, 2016) that leverages a quantity linked to the sum-based mean number of discoveries $E[N_e^{ab}]$. While they provide FDR control, it is only valid under some particular dependence conditions. It would nonetheless be interesting to extend their work to arbitrary dependence.

In this paper, the aim was to reliably detect one edge at a time. It would then be interesting to build upon the proposed edge-centric correlation thresholds to develop a multiple testing framework so as to provide theoretical control over the estimation of the connectivity of all pairs of region, perhaps by leveraging existing bootstrapping techniques (Cai and Liu, 2016).
Finally, it would be valuable to provide practitioners with a way to quantify edge detection uncertainty. For instance, confidence intervals for each edge of the whole inferred network could be defined. Some work has been done to determine confidence intervals for correlation coefficients in the bivariate case, both for underlying normality, e.g., (Ruben, 1966; Muirhead, 2005), and unknown distributions (Hu et al, 2020). It would be worth exploring how these methods could build upon our approach to extend them to a more general case in order to account for dependence.

**Supplementary information.** Proofs of the propositions are available in the appendix. Additional discussions and details can be found in the supplementary materials. Source code, including a notebook detailing how to reproduce the figures of this paper, is also available at: https://gitlab.inria.fr/q-func/csinference.

**Acknowledgments.** This work was supported by the project Q-FunC from Agence Nationale de la Recherche under grant number ANR-20-NEUC-0003-02 and the National Science Foundation under grant IIS-2135859.

## References

Achard S, Salvador R, Whitcher B, et al (2006) A resilient, low-frequency, small-world human brain functional network with highly connected association cortical hubs. The Journal of Neuroscience 26(1):63–72. 00831

Achard S, Coeurjolly J, Marcillaud R, et al (2011) fMRI functional connectivity estimators robust to region size bias. In: Proc. IEEE Workshop on Statistical Signal Processing (SSP), Nice, France, pp 813–816

Achard S, Delon-Martin C, Vértes P, et al (2012) Hubs of brain functional networks are radically reorganized in comatose patients. Proc Natl Acad Sci 109(50):20608–20613
Azriel D, Schwartzman A (2014) The empirical distribution of a large number of correlated normal variables. Journal of the American Statistical Association 110:00-00

Becq GGJPC, Barbier E, Achard S (2020a) Brain networks of rats under anesthesia using resting-state fmri: comparison with dead rats, random noise and generative models of networks. Journal of Neural Engineering

Becq GJP, Habet T, Collomb N, et al (2020b) Functional connectivity is preserved but reorganized across several anesthetic regimes. NeuroImage 219:116945

Boschi A, Brofiga M, Massobrio P (2021) Thresholding functional connectivity matrices to recover the topological properties of large-scale neuronal networks. Frontiers in Neuroscience 15

Cai TT, Liu W (2016) Large-scale multiple testing of correlations. Journal of the American Statistical Association 111(513):229-240. PMID: 27284211

Calhoun VD, Sui J, Kiehl K, et al (2012) Exploring the psychosis functional connectome: aberrant intrinsic networks in schizophrenia and bipolar disorder. Frontiers in psychiatry 2:75

Dadi K, Rahim M, Abraham A, et al (2019) Benchmarking functional connectome-based predictive models for resting-state fmri. NeuroImage 192:115-134

Di Martino A, Yan CG, Li Q, et al (2014) The autism brain imaging data exchange: towards a large-scale evaluation of the intrinsic brain architecture in autism. Molecular psychiatry 19(6):659

Donner A, Eliasziw M (1991) Methodology for inferences concerning familial correlations: a review. J Clin Epidemiol 44(4/5):449-455
Drton M, Perlman MD (2007) Multiple Testing and Error Control in Gaussian
Graphical Model Selection. Statistical Science 22(3):430 – 449

Dunlap W, Jones M, Bittner A (1983) Average correlations vs. correlated averages.
Bulletin of the Psychonomic Society 21:213–216

Elston RC (1975) On the correlation between correlations. Biometrika 62(1):133-140

Fallani FDV, Richiardi J, Chavez M, et al (2014) Graph analysis of functional brain
networks: practical issues in translational neuroscience. Philosophical Transactions
of the Royal Society B: Biological Sciences 369(1653):20130521

Fan J, Lv J (2008) Sure independence screening for ultrahigh dimensional feature
space. Journal of the Royal Statistical Society: Series B (Statistical Methodology)
70(5):849–911

Halliwell JW (1962) Dangers inherent in correlating averages. The Journal of
Educational Research 55(7):327-329

Hero A, Rajaratnam B (2011) Large scale correlation screening. Journal of the
American Statistical Association 106:1540–1552

Hotelling H (1953) New light on the correlation coefficient and its transforms. Journal
of the Royal Statistical Society Series B (Methodological) 15(2):193–232

Hu X, Jung A, Qin G (2020) Interval estimation for the correlation coefficient. The
American Statistician 74(1):29–36

Lehmann EL (1966) Some Concepts of Dependence. The Annals of Mathematical
Statistics 37(5):1137 – 1153
Malagurski B, Péran P, Sarton B, et al (2019) Topological disintegration of resting state functional connectomes in coma. *NeuroImage* 195:354–361

Muirhead RJ (2005) Aspects of Multivariate Statistical Theory. Wiley-Interscience

Ostroff C (1993) Comparing correlations based on individual-level and aggregated data. *Journal of Applied Psychology* 78:569–582

Panaretos VM, Zemel Y (2019) Statistical aspects of wasserstein distances. *Annual Review of Statistics and Its Application* 6(1):405–431

Petersen A, Müller HG (2019) Wasserstein covariance for multiple random densities. *Biometrika* 106(2):339–351

Poli D, Pastore VP, Massobrio P (2015) Functional connectivity in in vitro neuronal assemblies. *Frontiers in Neural Circuits* 9

Richiardi J, Achard S, Bunke H, et al (2013) Machine learning with brain graphs: Predictive modeling approaches for functional imaging in systems neuroscience. *IEEE Signal Process Mag* 30(3):58–70

Rosner B, Donner A, Hennekens C (1977) Estimation of interclass correlation from familial data. *Applied Statistics* 26:179–187

Ruben H (1966) Some new results on the distribution of the sample correlation coefficient. *Journal of the Royal Statistical Society Series B (Methodological)* 28(3):513–525

Rubinov M, Sporns O (2010) Complex network measures of brain connectivity: Uses and interpretations. *NeuroImage* 52(3):1059–1069. *Computational Models of the Brain*
van den Heuvel MP, de Lange SC, Zalesky A, et al (2017) Proportional thresholding
in resting-state fmri functional connectivity networks and consequences for patient-
control connectome studies: Issues and recommendations. NeuroImage 152:437-449

Whitcher B, Guttorp P, Percival DB (2000) Wavelet analysis of covariance with application to atmospheric time series. Journal of Geophysical Research 105(D11)(14):941-962

Whitfield-Gabrieli S, Nieto-Castanon A (2012) Conn: A functional connectivity toolbox for correlated and anticorrelated brain networks. Brain Connectivity 2(3):125-141. PMID: 22642651

Wigley TML, Briffa KR, Jones PD (1984) On the average value of correlated time series, with applications in dendroclimatology and hydrometeorology. Journal of Applied Meteorology and Climatology 23(2):201-213

Šidák Z (1967) Rectangular confidence regions for the means of multivariate normal distributions. Journal of the American Statistical Association 62(318):626-633

# Appendix A Proof of Proposition 1

## A.1 U-scores

Before proving Proposition 1, we need to introduce U-scores. *U-scores* are an orthogonal projection of the *Z-scores* of random variables. They are confined to an (*n* − 2)-sphere centered around 0 and with radius 1, denoted *S*<sub>*n*−2</sub>, with *n* the number of samples. We refer to (Hero and Rajaratnam, 2011) for a full definition. U-scores namely provide a practical expression of the correlation coefficient as an inner product of U-scores: *R*<sub>*i*,j*</sub><sup>ab</sup> = *U*<sub>*i*</sub><sup>*a*</sup><sup>*T*</sup>*U*<sub>*j*</sub><sup>*b*</sup> = 1 − ||*U*<sub>*i*</sub><sup>*a*</sup> − *U*<sub>*j*</sub><sup>*b*</sup>||<sup>2</sup>/2, where *U*<sub>*i*</sub><sup>*a*</sup>, *U*<sub>*j*</sub><sup>*b*</sup> are the random variables of the U-scores of voxels *i* and *j* in regions *a* and *b*, respectively, and ||*·*||<sup>2</sup> is
the squared Euclidean distance. Consequently, when U-scores are close to one another on $S_{n-2}$, they are associated with a high correlation.

## A.2 Intuition behind Proposition 1.

Roughly speaking, Proposition 1 means that if the Wasserstein distance between the densities of sample intra-correlation coefficients is sufficiently large, which means that the two distributions are highly different, then the average Euclidean distance between U-scores from the two regions is large too. Hence the average of inter-correlations is quite low. This phenomenon is illustrated in Figure A1 where two regions with different intra-correlation densities are depicted when $n = 3$.

<figure><img src="image_8.png" /><figcaption>Fig. A1: $S_{n-2}$ with $n = 3$ and U-scores from two regions (red diamonds and orange discs) that have a high intra-correlation density Wasserstein distance. Recalling that a high Euclidean distance between U-scores implies a low correlation, we can intuitively observe that the average inter-correlation is upper-bounded.</figcaption></figure>

## A.3 Proof of Proposition 1

Let us first remark:

$$ \inf_{c \in [0,1]} (F_{R^{a,a}}^{-1}(c) - F_{R^{b,b}}^{-1}(c))^2 \le d_W^2(f_{R^{a,a}}, f_{R^{b,b}}) \le \sup_{c \in [0,1]} (F_{R^{a,a}}^{-1}(c) - F_{R^{b,b}}^{-1}(c))^2 \quad (A1) $$

and as $(F_{R^{a,a}}^{-1} - F_{R^{b,b}}^{-1})^2$ is continuous on $[0,1]$, it attains its supremum and infimum.
Additionally, we can notice that for each $c \in [0, 1]$, there exist two points $U, V \in S_{n-2}$ such that $F_{R^{a,a}}^{-1}(c) = 1 - \|U - V\|^2/2$, and similarly for region b. Moreover, for all $x, y \in \mathcal{R}_a$, there exists $c \in [0, 1]$ such that, with their corresponding U-scores denoted $U_x, U_y$ (which are in $S_{n-2}$), $R_{x,y}^{a,a} = 1 - \|U_x - U_y\|^2/2 = F_{R^{a,a}}^{-1}(c)$, and analogously for region b.

Therefore, under the assumption $\min_{c \in [0,1]} (F_{R^{a,a}}^{-1}(c) - F_{R^{b,b}}^{-1}(c))^2 \ge A$, there exist $U_{x_a}, U_{y_a}, U_{x_b}, U_{y_b} \in S_{n-2}$ such that

$$ \min_{c \in [0,1]} (F_{R^{a,a}}^{-1}(c) - F_{R^{b,b}}^{-1}(c))^2 = \frac{1}{4} (\|U_{x_b} - U_{y_b}\|^2 - \|U_{x_a} - U_{y_a}\|^2)^2, $$

and it follows for all $v_a, w_a \in \mathcal{R}_a, v_b, w_b \in \mathcal{R}_b$,

$$ A \le \frac{1}{4} (\|U_{x_b} - U_{y_b}\|^2 - \|U_{x_a} - U_{y_a}\|^2)^2 \le \frac{1}{4} (\|U_{v_b} - U_{w_b}\|^2 - \|U_{v_a} - U_{w_a}\|^2)^2. $$

Thus, expanding the term on the right and applying the triangle inequality, followed by the reverse triangle inequality,

$$
\begin{aligned}
2\sqrt{A} &\le (\|U_{v_b} - U_{w_b}\| + \|U_{v_a} - U_{w_a}\|) \cdot \|U_{v_b} - U_{w_b}\| - \|U_{v_a} - U_{w_a}\| \\
&\le (\|U_{v_b} - U_{v_a}\| + \|U_{v_a} - U_{w_b}\| + \|U_{v_a} - U_{w_b}\| + \|U_{w_b} - U_{w_a}\|) \cdot \|U_{v_b} - U_{w_b} - (U_{v_a} - U_{w_a})\| \\
&\le (\|U_{v_b} - U_{v_a}\| + \|U_{v_a} - U_{w_b}\| + \|U_{v_a} - U_{w_b}\| + \|U_{w_b} - U_{w_a}\|) \cdot (\|U_{v_b} - U_{v_a}\| + \|U_{w_b} - U_{w_a}\|) \\
&\le (\|U_{v_b} - U_{v_a}\|^2 + \|U_{v_b} - U_{v_a}\| \cdot \|U_{w_b} - U_{v_a}\| + \|U_{w_b} - U_{v_a}\| \cdot \|U_{w_b} - U_{w_a}\| \\
&\quad + \|U_{v_b} - U_{v_a}\| \cdot \|U_{w_b} - U_{w_a}\|) \\
&\quad + (\|U_{w_b} - U_{w_a}\|^2 + \|U_{w_b} - U_{v_a}\| \cdot \|U_{v_b} - U_{v_a}\| + \|U_{w_b} - U_{w_a}\| \cdot \|U_{w_b} - U_{v_a}\| \\
&\quad + \|U_{v_b} - U_{v_a}\| \cdot \|U_{w_b} - U_{w_a}\|).
\end{aligned}
$$
We can then notice

$$
\begin{align*}
\overline{\|U^b - U^a\|^2} &= \left( \frac{1}{p_a p_b} \sum_{v_a \in \mathcal{R}_a} \sum_{v_b \in \mathcal{R}_b} \|U_{v_b} - U_{v_a}\| \right)^2 \\
&= \frac{1}{(p_a p_b)^2} \sum_{h_a \in \mathcal{R}_a} \sum_{h_b \in \mathcal{R}_b} \|U_{h_b} - U_{h_a}\|^2 + \\
&\quad \frac{1}{(p_a p_b)^2} \sum_{h_a \in \mathcal{R}_a} \sum_{h_b \in \mathcal{R}_b} \sum_{k_b \in \mathcal{R}_b - \{h_b\}} \|U_{h_b} - U_{h_a}\| \cdot \|U_{k_b} - U_{h_a}\| + \\
&\quad \frac{1}{(p_a p_b)^2} \sum_{h_a \in \mathcal{R}_a} \sum_{h_b \in \mathcal{R}_b} \sum_{k_a \in \mathcal{R}_a - \{h_a\}} \|U_{h_b} - U_{h_a}\| \cdot \|U_{h_b} - U_{k_a}\| + \\
&\quad \frac{1}{(p_a p_b)^2} \sum_{h_a \in \mathcal{R}_a} \sum_{h_b \in \mathcal{R}_b} \sum_{k_a \in \mathcal{R}_a - \{h_a\}} \sum_{k_b \in \mathcal{R}_b - \{h_b\}} \|U_{h_b} - U_{h_a}\| \cdot \|U_{k_b} - U_{k_a}\|.
\end{align*}
$$

Thus $\overline{\|U^b - U^a\|^2} \ge \frac{1}{(p_a p_b)^2} \cdot \frac{(p_a p_b)^2}{2} \cdot 2\sqrt{A} = \sqrt{A}$. From the Cauchy-Schwarz inequality,

$$
\overline{R^{a,b}} \le 1 - \frac{\overline{\|U^b - U^a\|^2}}{2},
$$

which completes the proof.

**Appendix B Proof of Proposition 3**

Let us recall Proposition 3.

**Proposition 3.** For all $\rho \in [0, 1]$,

$$
\hat{\nu}_e^{ab}(\rho) \le \hat{\nu}^{ab}(\rho). \qquad (B2)
$$

Proof. Since for all $\rho \in [0, 1]$, $0 \le \hat{F}_{|R^{a,b}|}(\rho) \le 1$, then $\hat{F}_{|R^{a,b}|}(\rho)^{p_b} \le \hat{F}_{|R^{a,b}|}(\rho)$. Thus, $\hat{\nu}^{ab}(\rho) = 1 - \hat{F}_{|R^{a,b}|}(\rho)^{p_b} \ge 1 - \hat{F}_{|R^{a,b}|}(\rho) = \hat{\nu}_e^{ab}$ $\square$
# Appendix C Additional insights on $\nu^{ab}$

We can derive the following proposition.

**Proposition 4.** If, for a fixed $i = 1, \dots, p_a$, all sample inter-correlation coefficients $R_{i,j}^{a,b}$ are i.i.d., $\nu^{ab} = \tilde{\nu}^{ab}$.

*Proof.* We can first remark $\max_{j \in \mathbb{R}_b} \phi_{ij}^{ab} = 1 - \prod_{j=1}^{p_b} (1 - \phi_{ij}^{ab})$. Thus,

$$ 
\begin{aligned}
E[N^{ab}] &= \sum_{i=1}^{p_a} E[1 - \prod_{j=1}^{p_b} (1 - \phi_{ij}^{ab})] \\
&= \sum_{i=1}^{p_a} (1 - \prod_{j=1}^{p_b} E[(1 - \phi_{ij}^{ab})]) \text{ under the assumption of independence.}
\end{aligned}
$$

For a fixed $i = 1, \dots, p_a$, under the assumption all $|R_{i,j}^{a,b}|$ are identically distributed, then, for all $j, l = 1, \dots, p_b$, $F_{|R_{i,j}^{a,b}|} = F_{|R_{i,l}^{a,b}|}$. Denote, $F_{|R^{a,b}|}$ the distribution function such that $F_{|R_{i,j}^{a,b}|} = F_{|R^{a,b}|}$ for all $i = 1, \dots, p_a$, $j = 1, \dots, p_b$. Thus, $p_a \cdot \tilde{\nu}^{ab} = E[N^{ab}] = p_a \cdot [1 - F_{|R^{a,b}|}^{p_b}] = p_a \cdot \nu^{ab}$. $\square$

# Appendix D Further information about our implementation

Our implementation is based on R 4.2.3. All experiments were performed on a laptop running on Ubuntu 18.04 with eight 1.8GHz 64-bits Intel Core i7-10610U CPUs, 32 GB of memory and a 1 TB hard drive.