# Inter-regional correlation estimators for functional magnetic resonance imaging

Sophie Achard, Jean-François Coeurjolly, Pierre Lafaye de Micheaux, Hanâ

Lbath, Jonas Richiardi

► To cite this version:

Sophie Achard, Jean-François Coeurjolly, Pierre Lafaye de Micheaux, Hanâ Lbath, Jonas Richiardi. Inter-regional correlation estimators for functional magnetic resonance imaging. NeuroImage, 2023, 282, pp.120388. 10.1016/j.neuroimage.2023.120388. hal-04242995

HAL Id: hal-04242995

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

Submitted on 15 Oct 2023

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 CC BY 4.0 - Attribution - International License
Inter-regional correlation estimators for functional
magnetic resonance imaging

2

3

4

Sophie Achard<sup>a</sup>, Jean-François Coeurjolly<sup>a</sup>, Pierre Lafaye de Micheaux<sup>b,c,d</sup>,
Hanâ Lbath<sup>a</sup>, Jonas Richiardi<sup>e</sup>

5

6

7

8

9

10

11

12

<sup>a</sup> Univ. Grenoble Alpes, CNRS, Inria, Grenoble-INP, LJK, 38000 Grenoble, France

<sup>b</sup>AMIS, Université Paul Valéry Montpellier 3, France

<sup>c</sup>PreMeDICaL - Precision Medicine by Data Integration and Causal Learning, Inria
Sophia Antipolis, France

<sup>d</sup>Desbrest Institute of Epidemiology and Public Health, Univ Montpellier, INSERM,
Montpellier, France

<sup>e</sup>Department of Radiology, Lausanne University Hospital and University of Lausanne,
Switzerland

13

Abstract

14

15

16

17

18

19

Functional magnetic resonance imaging (fMRI) functional connectivity be-
tween brain regions is often computed using parcellations defined by func-
tional or structural atlases. Typically, some kind of voxel averaging is per-
formed to obtain a single temporal correlation estimate per region pair. How-
ever, several estimators can be defined for this task, with various assumptions
and degrees of robustness to local noise, global noise, and region size.

20

21

22

23

24

25

26

In this paper, we systematically present and study the properties of 9
different functional connectivity estimators taking into account the spatial
structure of fMRI data, based on a simple fMRI data spatial model. These
include 3 existing estimators and 6 novel estimators. We demonstrate the em-
pirical properties of the estimators using synthetic, animal, and human data,
in terms of graph structure, repeatability and reproducibility, discriminabil-
ity, dependence on region size, as well as local and global noise robustness.
We prove analytically the link between regional intra-correlation and inter-region correlation, and show that the choice of estimator has a strong influence on inter-correlation values. Some estimators, including the commonly used *correlation of averages* (CA), are positively biased, and have more dependence to region size and intra-correlation than robust alternatives, resulting in spatially-dependent bias. We define the new *local correlation of averages* estimator with better theoretical guarantees, lower bias, significantly lower dependence on region size (Spearman correlation 0.40 vs 0.55, paired t-test $T=27.2, p = 1.1e^{-47}$), at negligible cost to discriminative power, compared to the CA estimator.

The difference in connectivity pattern between the estimators is not distributed uniformly throughout the brain, but rather shows a clear ventral-dorsal gradient, suggesting that region size and intra-correlation plays an important role in shaping functional networks defined using the CA estimator, and leading to non-trivial differences in their connectivity structure. We provide an open source R package and equivalent Python implementation to facilitate the use of the new estimators, together with preprocessed rat time-series.

*Keywords:* functional connectivity, correlation, aggregated data, familial correlations, serial correlations

# 1. Introduction

Functional connectivity of the brain is estimated from observations using non invasive techniques such as electroencephalography (EEG), magnetoencephalography (MEG) or functional Magnetic Resonance Imaging (fMRI).
Each recording provides time series associated to spatial locations within regions of the brain. Functional connectomes, that is, graphs representing the estimated connectivity, are then constructed by computing dependence measures between the time series. These connectomes are used in fundamental neuroscience, for example to study development (Fan et al., 2021; Tooley et al., 2021), and in clinical neuroscience, to characterize psychiatric (Fornito et al., 2017) or neurological (Fornito et al., 2015) disorders. They form a compact yet expressive representation of brain activity that can be used for downstream analysis tasks such as diagnosis (Castellanos et al., 2013) or more generally machine learning approaches (Richiardi et al., 2013; Dadi et al., 2019).

For graph construction, typically, each region of the brain, defined by a structural or functional parcellation, is associated to a given set of voxels amongst the thousands for which a signal is recorded. The idea is then to extract a representative of the set of voxels to attach one time series to each region. When structural atlases are used, the most common approach is to take the average of the voxel time series at each time point. Indeed, almost 70% of papers on PubMed that used the Human Connectome Project dataset to conduct functional-connectivity-related studies in the last five years (e.g., Ogawa (2021); Figueroa-Jimenez et al. (2021); Bolt et al. (2017); Zhang et al. (2016)) use this method (cf. Appendix Appendix A). While there are numerous other approaches to connectivity estimation (including the related techniques for estimating parcellation from connectivity, see, e.g., Eickhoff et al. (2015), or using regional medians (Braun et al., 2012) or eigenvectors (Büchel and Friston, 1997; Braun et al., 2012) instead of means), we
focus on correlations between averaged regional time-courses, and argue that
this technique introduces bias in the estimation of the functional connec-
tomes.

The choice of acquisition sequence and hardware, physiological noise (Caballero-Gaudes and Reynolds, 2017), preprocessing (Braun et al., 2012), and subject motion all impact correlation estimators. Acquisition effects are site-dependent, causing heterogeneity problems in multi-site studies, although various harmonization techniques have been proposed as mitigationistrillon et al. (2014); ?. In addition, it has been shown that computation of connectomes is affected by three main parameters: the length of the acquisition (Whitlow et al., 2011; Van Dijk et al., 2010), the number of regions (Stanley et al., 2013; Cao et al., 2019) and the chosen frequency band (Cordes et al., 2001; Salvador et al., 2005; Braun et al., 2012; Chen and Glover, 2015). Finally, sample size will also play a role in terms of group comparisons (Termenon et al., 2016).

Aggregation across voxels is often used because one wishes to increase
the signal to noise ratio. This approach is also common in other areas of
statistical analysis, for instance because the data are collected in different
groups, organizations, or regions. However, difficult challenges arise due to
the presence of correlations within the collected datasets.

*Measurement errors* can impact inter-region correlations. They have been
well studied in fMRI, and are related to both the hardware and the sub-
ject (Greve et al., 2013). They are known to impart correlation structure
to the data that is not linked to neural activity (Jo et al., 2010; Murphy
et al., 2013), at various spatial scales. This problem is also common in other areas of statistical analysis. For example, Ostroff (1993) studied correlations between the score variables job satisfaction and technology (i.e., perception of the amount of standardization of tasks performed) both at the individual and organizational levels (individuals grouped into organizations such as companies). When organization-level estimates of correlation (i.e., correlations based on aggregated data) are obtained by averaging individual-level estimates of correlations, they showed that the ratio of individual to organizational correlations varied widely depending on measurement errors as well as other factors. This is an instance of the fallacy of the wrong level, when “correlations at a more macro level are used to make inference about individuals, or vice versa” Ostroff (1993).

Region size can also influence inter-region correlations. In fMRI, the dependence of inter-region correlations on brain region size has been noted (Achard et al., 2011), showing a positive relationship between region size and correlation values to the rest of the brain. This is particularly problematic because many atlases have some dependency between region size and spatial location (e.g. some deep gray matter structure may be parcellated into smaller regions or subregions than cortical structures). It has also been shown that temporal autocorrelation increases with region size, both for volume-based and surface-based parcellation (Afyouni et al., 2019), and that at a small scale regional homogeneity (Zang et al., 2004) – also called local connectivity and measured by the Kendall correlation coefficient of small neighbourhood time series – is larger for 9-voxel than for 27-voxels neighbourhoods (Jiang and Zuo, 2016).
In studies of familial data (Rosner et al., 1977), specific characteristics are obtained for different families with different sizes: *correlation between the children and parents* and the *average of correlations between all children and parents* are not equal in the majority of cases, due to differing number of children per family.

Finally, the *spatial aspect* of the data also complicates correlation estimation, in particular due to spatial autocorrelation between voxels (more simply named spatial correlation in the rest of the paper). Spatial correlation means that observations in neighbouring voxels are not independent. However, independence is an assumption many statistical estimators rely on to simplify hypothesis testing by enabling the application of the central limit theorem, leading to false positives and artificially low *p*-values. **Spatial correlation has been identified to be present in fMRI data previously, in particular in activation studies. Even if the methods to take into account the spatial correlation are different from the functional connectivity, it is worth detailing these specific approaches here. The common point between our problem in this paper and the activation studies is that spatial correlation has an impact on the design of the methodological approaches. Indeed, two classical approaches in fMRI activation studies are to either assume voxel independence, or conversely to smooth the data (Hartvig and Jensen, 2000). Smoothing itself can lead to location and amplitude mis-estimation (Descombes et al., 1998); alternatively, estimating smoothness from data allows adjusting effective number of degrees of freedom and reducing false positives, but can itself lead to variability in *p*-values from hypothesis tests (Poline et al., 1995). Spatial correlation has also been shown to deflate *p*-values in structural imaging (Burt
et al., 2020), and has long been recognized as an issue in functional connectivity, for instance with early voxelwise (PET) connectivity approaches removing correlation “between neighbouring voxels which can be attributed to spatial correlation” (Cao and Worsley, 1999). More recently, methods from spatial statistics have been applied for clustering fMRI data (Ye et al., 2009, 2011), and spatial correlation-preserving null models have been proposed to compare functional network maps (Alexander-Bloch et al., 2018; Markello and Misic, 2021) and thus avoid deflated p-values due to spatial correlation. Computing correlations is also common when geostatistics is applied to ecology, geography, climate studies, and more. The data collected in these fields are attached to a spatial position and usually with spatial correlation. The problem was first reported by Student (1914), and studied in (Clifford et al., 1989) for two spatial processes. Applications of these methods can be found for example in the study of meteorological data (Gunst, 1995).

In all these fields of application, the main difficulty is to take into account
spatial correlation when the goal is to construct estimators of correlation and
to build testing procedures when the averaged variables are not independent.

In light of the above, preferable inter-correlation estimators should exhibit
at least four properties i) face validity, ii) high repeatability, iii) preservation
of the differences between individuals (discriminative power), and iv) inde-
pendence from region size. In this paper, we question the default choice of
using correlations of averages of voxel timecourses, and examine in details
the assumptions of various methods and their robustness to various types of
noise. We propose first a simple definition of a spatial model of fMRI with
intra-correlations within brain regions. Then, computations of correlations are described and we show the potential bias in the estimators. Based on simulations, we illustrate the good behaviour of the newly introduced estimators. Finally, we conclude with results on human data and rat data.

## 2. Proposed estimators of correlation

### 2.1. Definition of the proposed spatial model for fMRI data

Let $\mathcal{C} \subset \mathbb{Z}^d, d \in \mathbb{N}^*$, be a finite compact set of multi-indices. In the context of our application, $d = 3$ and $\mathcal{C}$ contains all 3-tuples indexing the voxels of a three-dimensional image of a brain. Each brain is virtually partitioned into $J$ regions of interest which are represented through their subsets of voxels $\mathcal{R}_j$ of cardinality $\#\mathcal{R}_j = N_j, j = 1, \dots, J$. We thus have

$$ \mathcal{C} = \cup_{j=1}^{J} \mathcal{R}_j \text{ and } \#\mathcal{C} = \sum_{j=1}^{J} N_j. $$

For any $d$-tuple $i \in \mathcal{C}$, a signal $Y_i(\cdot)$ sampled at times $t = 1, \dots, T$ is observed and we assume that it can be decomposed as follows

$$ Y_i(t) = X_i(t) + \varepsilon_i(t) + e(t), \qquad (1) $$

where $X_i(\cdot)$ is an unobserved signal of interest, $\varepsilon_i(\cdot)$ is a local noise contaminating locally the signal $X_i(\cdot)$, and $e(\cdot)$ is a global noise corrupting in the same way the signal measured in each voxel indexed by an element of $\mathcal{C}$. This can be, e.g., a consequence of thermal or background noise (Lazar, 2008; Greve et al., 2013), which at high field strength has been shown to explain a high proportion of noise variance (Greve et al., 2011).
We make a few assumptions on these different components. First, we assume that all random variables are centered. We also assume that the signals $X_i(\cdot)$, $\varepsilon_i(\cdot)$ and $e(\cdot)$ are mutually independent and independent in time. This is not an overly restrictive assumption as for the applications presented in Section 3 we preprocess the data by applying a wavelet transform. And it is now well-known (Moulines et al., 2007), that if a random time series has short or long memory characteristics, after a wavelet decomposition, this signal can be approximated to be decorrelated in time for large wavelet scales. In addition, assuming that the $X_i$'s are centered is coherent as it is a well-known fact that a wavelet decomposition based on a wavelet mother with $K$ vanishing moments cancels out every polynomial trend with degree $K-1$. Assuming that the local and global noises are homoskedastic with a variance denoted respectively by $\sigma_{\varepsilon}$ and $\sigma_e$ is also not restrictive. Finally, we will be assuming (again following the literature Lazar (2008); Greve et al. (2011, 2013)) that the local noise $\varepsilon$ is not too strongly spatially dependent and more precisely that two voxels far away have uncorrelated local noise. This is made more precise in Section 2.2.

## 2.2. Spatial correlation structure induced by model (1)

Let $i, i' \in \mathcal{C}$, $j, j' = 1, \dots, J$ ($j \neq j'$) and for all $t = 1, \dots, T$, we assume that there exists $\sigma_j > 0$, $\sigma_{\varepsilon} \ge 0$, $r_{jj'} \in [-1, 1]$, $\rho_{ii'} \in (0, 1]$, $\eta_{ii'} \in [-1, 1]$ such that

$$ E[X_i(t)X_{i'}(t')] = \begin{cases} \sigma_j\sigma_{j'}r_{jj'} & \text{if } i \in \mathcal{R}_j, i' \in \mathcal{R}_{j'}, j \neq j', \\ \sigma_j^2\rho_{ii'} & \text{if } i, i' \in \mathcal{R}_j \end{cases} $$

and

$$ E[\varepsilon_i(t)\varepsilon_{i'}(t)] = \sigma_{\varepsilon}^2\eta_{ii'} \text{ if } i, i' \in \mathcal{R}_j. $$
The parameter $r_{jj'}$ represents the correlation between two (unobserved) signals of two different regions $\mathcal{R}_j$ and $\mathcal{R}_{j'}$ and is called inter-correlation between regions $\mathcal{R}_j$ and $\mathcal{R}_{j'}$ in the following. This is the parameter of interest we focus on in the rest of the paper. The parameter $\rho_{ii'}$ represents the intra-correlation between two (unobserved) signals inside a common region. Moreover, the parameter $\eta_{ii'}$ represents the spatial correlation between two local noises inside a common region. We assume that inside each region, the signals of interest have positive intra-correlation. This is supported by literature (Uddin et al., 2008; Tomasi and Volkow, 2010; Jiang and Zuo, 2016). We also assume that for each time $t$ and for $j = 1, \dots, J$, $\{X_i(t), i \in \mathcal{R}_j\}$ (resp. $\{\varepsilon_i(t), i \in \mathcal{C}\}$) is a second-order stationary and isotropic (with respect to the uniform norm) random field defined over $\mathcal{R}_j$ (resp. $\mathcal{C}$). This means in particular that both the correlations $\rho_{ii'}$ (for any $i, i' \in \mathcal{R}_j$ for some $j$) and $\eta_{ii'}$ (for $i, i' \in \mathcal{C}$) depend only on the (uniform) distance (that is the usual distance on the lattice, e.g. Gaetan et al. (2010)) between the two voxels $i$ and $i'$. For brevity, we still denote $\rho_{|i'-i|}$ by $\rho_{ii'}$ and $\eta_{|i'-i|}$ by $\eta_{ii'}$ where for $x \in \mathbb{Z}^d$, the notation $|x|$ stands for the uniform norm. Our a priori hypothesis is that the intra-correlation $\rho_\delta$ is close to 1 for moderate distances $\delta$, meaning that close neighbours are strongly (positively) correlated. Finally, we assume that the local noise is $p$-dependent, i.e., $\eta_\delta = 0$ for any $\delta > p$. Without loss of generality, we also intrinsically assume that for all $i \in \mathcal{R}_j$ and $i' \in \mathcal{R}_{j'}$, $\varepsilon_i(t)$ and $\varepsilon_{i'}(t)$ are uncorrelated. This simplifies the presentation in the next sections. Furthermore, in the sequel, employing a slight abuse of language, we refer to the correlation between two voxels instead of the correlation between the signals originating from those voxels.
Hence, this results in the following (spatial) correlation structure for the
signals $Y_i$ and $Y'_i$ at time $t$:

$$
\mathrm{E}[Y_i(t)Y_{i'}(t)] = \begin{cases}
\sigma_j \sigma_{j'} r_{jj'} + \sigma_e^2 & \text{if } i \in \mathcal{R}_j, i' \in \mathcal{R}_{j'}, j \neq j', \\
\sigma_j^2 \rho_{|i-i'|} + \sigma_e^2 \eta_{|i-i'|} + \sigma_e^2 & \text{if } i, i' \in \mathcal{R}_j \text{ and } |i - i'| \le p \\
\sigma_j^2 \rho_{|i-i'|} + \sigma_e^2 & \text{if } i, i' \in \mathcal{R}_j \text{ and } |i - i'| > p.
\end{cases}
$$

Given a parcellation of the brain, the objective is to estimate inter-
correlations $r_{jj'}$ for each pair of regions of interest, independently of the
parameters $\sigma_j, \sigma_j', \sigma_\epsilon, \sigma_e, \rho_{ii'}, \eta_{ii'}$ which are viewed as nuisance parameters.
We do not model the distribution of $Y_i$ but only its second-order properties
(through $X_i, \epsilon_i, e$). As said before, we consider the intra-correlations, the
correlation of the local noise and the different variances as nuisance parame-
ters that we do not want to estimate. In the next section, we present various
estimators of $r_{jj'}$ built in order to address one (or several) of the following
cases: (1) the regions of interest $\mathcal{R}_j$ and $\mathcal{R}_{j'}$ may contain a different number
of voxels; (2) the intra-correlation may deviate strongly from 1 (especially
for large regions); (3) there may be a non negligible local noise $\epsilon_i$ affecting
the signal in each region; (4) there may be a global noise affecting all regions.

2.3. Inter-correlation: notation and properties

Let $\mathbf{Y}_1 = (Y_1(1), \dots, Y_1(T))$ and $\mathbf{Y}_2 = (Y_2(1), \dots, Y_2(T))$ denote two
voxel time-series of length $T$. The notation $\widehat{\text{Cov}}(\mathbf{Y}_1, \mathbf{Y}_2)$, $\widehat{\text{Cor}}(\mathbf{Y}_1, \mathbf{Y}_2)$ and
$\widehat{\sigma}^2(\mathbf{Y}_1)$ stand for the sample covariance between $\mathbf{Y}_1$ and $\mathbf{Y}_2$, the sample
correlation between $\mathbf{Y}_1$ and $\mathbf{Y}_2$ and the sample variance of $\mathbf{Y}_1$, respectively.
For any $j = 1, \dots, J$, we define a $\nu$-neighbourhood and denote it by $\mathcal{V}$ as a
subset of $n_\nu := (2\nu + 1)^d$ indices, all of which are at a distance less than or
equal to $\nu$ from the center $j$ of the neighbourhood. For any set of indices $E$ (which could be a $\nu$-neighborhood or a region of interest) and any spatio-temporal field $Z_i(t)$ (which could be $Y_i, X_i, \varepsilon_i, \dots$) we denote by $\bar{Z}_E(t)$ for $t = 1, \dots, T$ the time series spatially averaged over $E$, that is

$$ \bar{Z}_E(t) = \frac{1}{\#E} \sum_{i \in E} Z_i(t). $$

To sum up, we reserve the bold notation to mainly denote a vector of length $T$, the notation $\hat{\cdot}$ to denote an average over time while $\cdot$ will denote an average over space. Hence, for instance $\hat{\sigma}^2(\bar{\mathbf{Y}}_E)$ denotes the sample variance of the vector with components $(\#E)^{-1} \sum_{i \in E} Y_i(t)$ for $t = 1, \dots, T$. We also let

$$ \bar{\rho}_E = \frac{1}{(\#E)^2} \sum_{i,i' \in E} \rho_{ii'} \quad \text{and} \quad \bar{\eta}_E = \frac{1}{(\#E)^2} \sum_{i,i' \in E} \eta_{ii'}. \qquad (2) $$

The quantity $\bar{\rho}_E$ represents the (spatial) average intra-correlation inside the set $E$. If $E$ corresponds to a $\nu$-neighborhood with moderate $\nu$, we may expect $\bar{\rho}_\nu$ to be close to 1. Such an observation is probably less realistic when $E = \mathcal{R}_j$ especially for large regions. The quantity $\bar{\eta}_E$ is related to the (spatial) correlation structure of the local noise. By assuming this noise to be $p$-dependent (that is $\eta_\delta = 0$ when $\delta \ge p$), it is clear that the larger $\#E$ the smaller $\bar{\eta}_E$.

Using the assumption given in Section 2.1, for any $E \subseteq \mathcal{R}_j$ and $E' \subseteq \mathcal{R}_{j'}$, we deduce

$$ \mathrm{Cov}[\bar{Y}_E(t), \bar{Y}_{E'}(t)] = \begin{cases} \sigma_j \sigma_{j'} r_{jj'} + \sigma_e^2, & \text{if } j \neq j', \\ \sigma_j^2 \bar{\rho}_{E,E'} + \sigma_e^2, & \text{if } j = j', \end{cases} $$
where

$$ \bar{\rho}_{E,E'} = \frac{1}{(\#E)(\#E')} \sum_{i \in E, i' \in E'} \rho_{|i-i'|}. $$

The variance can also be deduced as follows:

$$ \mathrm{Var}[\bar{Y}_E(t)] = \sigma_j^2 \bar{\rho}_E + \sigma_\epsilon^2 \bar{\eta}_E + \sigma_e^2. $$

The detail of this result is given in Proposition Appendix C.1.

To lighten the expression of estimators proposed in the next sections, we define for $j, j' \in \{1, \dots, J\}$

$$ \sigma_{\epsilon,j} = \frac{\sigma_{\epsilon}}{\sigma_j}, \quad \sigma_{e,j} = \frac{\sigma_e}{\sigma_j}, \quad \text{and} \quad \sigma_{e,jj'} = \frac{\sigma_e}{\sqrt{\sigma_j \sigma_{j'}}}. \qquad (3) $$

In the next sections, we set $j, j'$ and thus aim to estimate $r_{jj'}$ independently of the other parameters. The definition of standard estimators as well as novel estimators may look complicated due to the large amount of notation induced by the spatio-temporal correlation structure of $Y$ and the methods themselves. However, we have postponed as much as possible theoretical contents to Appendix and present the estimation methods from an intuitive point of view in Figures 1-2 in order to make next sections readable and reproducible.

## 2.4. Existing inter-correlation estimators

We first review existing inter-regional correlation estimators using a unified notation throughout<fnref n="1" />. Results on consistency of the estimators are provided in Appendix C-Appendix G.4.

<footnote n="1">¹In a previous study (Achard et al., 2011), we already described three of the estimators discussed here (CA, AC, $\ell_{CA}$), but not with a unified notation, as well as a fourth estimator which is only discussed in the appendix of the present paper</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1: Graphical overview of the inter-regional correlation estimators CA, AC, $\ell_{CA}$, and R discussed in this paper. Gray dots represents voxels. Dashed black lines represent brain regions. Edges between voxels represent voxel-voxel temporal correlations. Blue rectangles show the region and level of aggregation (voxels or correlations). $r_{jj'}$ shows quantities involved in the computation of the inter-regional correlations. Illustrations are approximate, please refer to the relevant equations for the exact definition. Neighbourhood versions of the estimators (starting with $\ell$) use the same principles but involve aggregating in small neighbourhoods within regions.</figcaption></figure>
### 2.4.1. Correlation of averages (method CA)

In order to increase the signal-to-noise ratio, the most standard method in fMRI is to average (or sometimes convolve with a Gaussian kernel) the signal in space (in each region of interest). The aggregated correlation estimator corresponds to the standard estimator (see Section 1) considered for example in Achard et al. (2006), Bolt et al. (2017) or Ogawa (2021):

$$ \hat{r}_{jj'}^{\text{CA}} = \frac{\widehat{\text{Cov}}(\bar{\mathbf{Y}}_{\mathcal{R}_j}, \bar{\mathbf{Y}}_{\mathcal{R}_{j'}})}{\hat{\sigma}(\bar{\mathbf{Y}}_{\mathcal{R}_j})\hat{\sigma}(\bar{\mathbf{Y}}_{\mathcal{R}_{j'}})}. \quad (4) $$

This estimator, illustrated in Figure 1 was designed to reduce the local noise. Indeed, in the absence of global noise ($\sigma_e^2 = 0$), this estimator tends to $r_{jj'}/\sqrt{(\bar{\rho}_{\mathcal{R}_j} + \sigma_{\varepsilon,j}^2\bar{\eta}_{\mathcal{R}_j})(\bar{\rho}_{\mathcal{R}_{j'}} + \sigma_{\varepsilon,j'}^2\bar{\eta}_{\mathcal{R}_{j'}})}$. The interest of averaging before computing correlations is clear: the local noise is smoothed, thus $\bar{\eta}_{\mathcal{R}_j} = O(1/N_j)$ is probably very small. However, even in absence of noise ($\sigma_\varepsilon = \sigma_e = 0$), $\hat{r}_{jj'}^{\text{CA}}$ has a serious drawback since it estimates $r_{jj'}/\sqrt{\bar{\rho}_{\mathcal{R}_j}\bar{\rho}_{\mathcal{R}_{j'}}}$ which is highly dependent on intra-correlation. Just to give an example, assume $r_{jj'} = 1/2$, $N_j = N_{j'} = 2$, $\rho_1 = 0$ then $\bar{\rho}_{\mathcal{R}_j} = \bar{\rho}_{\mathcal{R}_{j'}} = 1/2$ then $\hat{r}_{jj'}^{\text{CA}}$ will converge towards 1 and not 1/2. This is a caricature but illustrates what may happen for large regions when some of the signals $X_i$ are not enough positively intra-correlated. That fact was already pointed out by Achard et al. (2011).

### 2.4.2. Average of correlations (method AC)

Instead of evaluating correlation of spatial averages, it is natural to perform the (spatial) average of correlations. This estimator, illustrated in Fig-
306 ure 1, is given by:

$$ \hat{r}_{jj'}^{AC} = \frac{1}{N_j N_{j'}} \sum_{i \in \mathcal{R}_j, i' \in \mathcal{R}_{j'}} \widehat{\mathrm{Cor}}(\mathbf{Y}_i, \mathbf{Y}_{i'}) . \quad (5) $$

As seen from Table 1, in absence of global noise ($\sigma_e^2 = 0$) and when the variances are equal to 1, this estimator estimates the quantity $r_{jj'}/(1 + \sigma_e^2)$ which makes this estimator robust to large regions (for which $\bar{\rho}_{R_j}$ may be far from 1) but more sensitive to local noise than the estimator $\hat{r}_{jj'}^{CA}$.

### 2.4.3. Replicates for correlations (method R)

In order to cancel out the effect of local noise, we introduce to fMRI a slight adaptation of the estimator introduced by Bergholm et al. (2010), in the context of image analysis. This is based on the concept of replicates within the same region, and denoted by $\hat{r}^R$ (R for replicates). The idea is to take two samples within each region, called replicates, to be able to compute correlation using these replicates to cancel out the local noise. These replicates can be chosen randomly a certain number of times denoted B. This estimator, illustrated in Figure 1, is then obtained as a Monte-Carlo mean (or bootstrap) over different random replicates

$$ \hat{r}_{jj'}^{R} = \frac{1}{B} \sum_{b=1}^{B} \frac{\frac{1}{4} \sum_{\alpha, \beta=1}^{2} \widehat{\mathrm{Cor}}(\mathbf{Y}_{i_{\alpha}^{(b)}}, \mathbf{Y}_{i_{\beta}^{(b)}})}{\sqrt{|\widehat{\mathrm{Cor}}(\mathbf{Y}_{i_{1}^{(b)}}, \mathbf{Y}_{i_{2}^{(b)}}) \widehat{\mathrm{Cor}}(\mathbf{Y}_{i_{1}^{(b)}}, \mathbf{Y}_{i_{2}^{(b)}})|}} \quad (6) $$

where for $b = 1, \dots, B$, $i_1^{(b)}, i_2^{(b)} \in \mathcal{R}_j$ are such that $|i_2^{(b)} - i_1^{(b)}| = \delta \ge p$. In the same way, $i_1'^{(b)}, i_2'^{(b)} \in \mathcal{R}_{j'}$ are such that $|i_2'^{(b)} - i_1'^{(b)}| = \delta \ge p$. Under equal variances and absence of global noise, $\hat{r}_{jj'}^R$ estimates $1/|\rho_\delta|$ which is clearly independent of $\sigma_e^2$ and may be expected to be close to one if $\delta$ is small.
2.5. Novel estimators: discarding the effect of global and/or local noise

2.5.1. Use of a priori uncorrelated regions (method D based on differences)

We now present an estimator which handles the problem of global noise. To achieve this task, we assume that among the regions where the signal is recorded there are at least two regions say $\mathcal{R}_k$ and $\mathcal{R}_{k'}$ which are uncorrelated between themselves and from all the other ones. With a slight abuse of notation, $k, k'$ will be used for the indices of these two regions, while $j, j'$ will be used when we are interested in the inter-correlation between regions $\mathcal{R}_j$ and $\mathcal{R}_{j'}$ (hence $r_{jk} = r_{jk'} = r_{j'k} = r_{j'k'} = 0$). This assumption is realistic in the context of fMRI data where we are interested in the correlations between cortical regions. Indeed, the field of view is typically larger than the brain itself, and the definition of extra regions is possible, for instance using air voxels or muscle voxels. The estimator is illustrated in Figure 2.

We propose the following strategy: for $b = 1, \dots, B$ let $i^{(b)}, i'^{(b)}, k^{(b)}$ and $k'^{(b)}$ be voxels of $\mathcal{R}_j, \mathcal{R}_{j'}, \mathcal{R}_k$ and $\mathcal{R}_{k'}$.

$$ \hat{r}_{jj'}^{\text{D}} = \frac{1}{B} \sum_{b=1}^{B} \widetilde{\text{Cor}}(\mathbf{Y}_{i^{(b)}}, \mathbf{Y}_{i'^{(b)}}, \mathbf{Y}_{k^{(b)}}, \mathbf{Y}_{k'^{(b)}}), \quad (7) $$

where for four vectors $\mathbf{Y}_1, \mathbf{Y}_2, \mathbf{Y}_3$ and $\mathbf{Y}_4$ (with same length)

$$ \widetilde{\text{Cor}}(\mathbf{Y}_1, \mathbf{Y}_2; \mathbf{Y}_3, \mathbf{Y}_4) = \frac{\widehat{\text{Cov}}(\mathbf{Y}_1 - \mathbf{Y}_3, \mathbf{Y}_2 - \mathbf{Y}_4)}{\widehat{s}(\mathbf{Y}_1, \mathbf{Y}_3, \mathbf{Y}_4) \widehat{s}(\mathbf{Y}_2, \mathbf{Y}_3, \mathbf{Y}_4)} \quad (8) $$

and where for three vectors $\mathbf{U}, \mathbf{V}$ and $\mathbf{W}$ with same length

$$ \widehat{s}^2(\mathbf{U}, \mathbf{V}, \mathbf{W}) = (\widehat{\sigma}^2(\mathbf{U} - \mathbf{V}) + \widehat{\sigma}^2(\mathbf{U} - \mathbf{W}) - \widehat{\sigma}^2(\mathbf{V} - \mathbf{W})) / 2. $$

The intuition of this estimator is quite simple. Assume that the local noise has zero variance. Since the noise $e(\cdot)$ is global, subtracting from $Y_{i^{(b)}}(t)$ the
value $Y_{k(b)}(t)$ and from $Y_{i'(b)}(t)$ the value $Y_{k'(b)}(t)$ discards the global noise. And since the regions $\mathcal{R}_k$ and $\mathcal{R}_{k'}$ are not correlated and not correlated to the other ones, the numerator (for each $b$) is an estimate of $\sigma_j\sigma_{j'}r_{jj'}$. Then, we just have to divide by estimates of $\sigma_j$ (and $\sigma_{j'}$). We observe that this cannot be done using simply $\hat{\sigma}^2(\mathbf{Y}_{i(b)} - \mathbf{Y}_{k(b)})$ which estimates $\sigma_j^2 + \sigma_k^2 + 2\sigma_{jk}\sigma_{kj'}$. This justifies the introduction of $\hat{s}^2$.

Note that $\hat{r}_{jj'}^D$ is still biased with respect to local noise (see Table 1).

An illustration of estimator D is provided in Figure 2, and a more formal proposition and proof for this estimator are provided in Appendix E.

### method D

<figure><img src="image_3.png" /><figcaption>Figure 2: Main steps involved in computing the inter-regional correlation estimator D. Dashed black lines represent brain regions of interest $\mathcal{R}_j$ and $\mathcal{R}_{j'}$. Dotted black lines represent a priori uncorrelated region (e.g. air or muscle voxels) $\mathcal{R}_k$ and $\mathcal{R}_{k'}$. Colored rectangles show the voxels involved in the computation. The final inter-regional correlation estimator $r_{jj'}$ is defined in terms of the intermediate quantities computed in these three steps. See Figure 1 for more details and other estimators.</figcaption></figure>
### 2.5.2. Replicates and use of a priori disconnected regions: method RD

Combining replicates and the idea based on differences motivates us to propose the following estimator (see Sections 2.4.3 and 2.5.1 for notation)

$$ \hat{r}_{jj'}^{\text{RD}} = \frac{1}{B} \sum_{b=1}^{B} \frac{\frac{1}{4} \sum_{\alpha, \beta=1}^{2} \widetilde{\text{Cor}}(\mathbf{Y}_{i_{\alpha}^{(b)}}, \mathbf{Y}_{i'_{\beta}^{(b)}}; \mathbf{Y}_{k^{(b)}}, \mathbf{Y}_{k'_{(b)}})}{\sqrt{|\widetilde{\text{Cor}}(\mathbf{Y}_{i_{1}^{(b)}}, \mathbf{Y}_{i_{2}^{(b)}}; \mathbf{Y}_{k^{(b)}}, \mathbf{Y}_{k'_{(b)}}) \widetilde{\text{Cor}}(\mathbf{Y}_{i'_{1}^{(b)}}, \mathbf{Y}_{i'_{2}^{(b)}}; \mathbf{Y}_{k^{(b)}}, \mathbf{Y}_{k'_{(b)}})|}} \quad (9) $$

It is worth pointing out that $r_{jj'}^{\text{RD}}$ is independent of $\sigma_{\varepsilon}$ and $\sigma_{e}$ and equals the unknown $r_{jj'}$ if $\rho_{\delta}$ is close to 1. A more formal proposition and proof for this estimator are provided in Appendix F.

### 2.6. Localized versions of inter-correlation estimators

As mentioned previously, when noisy signals are averaged, the signal to noise ratio increases. A very popular method in neuroimaging analyses is to apply a Gaussian smoothing on the fMRI volumes (Worsley et al., 1992, 1996; Poline et al., 1997). However, applying a large kernel width may have dramatic effect on brain connectivity (Triana et al., 2020). Some earlier work on PET connectivity used a local neighbourhood centered around voxels of interest to smooth the signal in each region prior to connectivity estimation (Köhler et al., 1998). We introduce in this section estimators using local neighborhoods to control the smoothing effect on correlation estimations.

#### 2.6.1. Local correlation of averages (method $\ell_{CA}$)

Motivated by the first two estimators, we propose to estimate $r_{jj'}$ using an empirical average of local spatial averages. For $b = 1, \dots, B$, let $\mathcal{V}_{j}^{(b)}$ (resp.
372 $\mathcal{V}_{j'}^{(b)}$ be a $\nu$-neighborhood of $\mathcal{R}_j$ (resp. $\mathcal{R}_{j'}$). We define

$$ \hat{r}_{jj'}^{\ell_{CA}} = \frac{1}{B} \sum_{b=1}^{B} \widehat{\operatorname{Cor}}(\bar{\mathbf{Y}}_{\mathcal{V}_{j}^{(b)}}, \bar{\mathbf{Y}}_{\mathcal{V}_{j'}^{(b)}}). \quad (10) $$

### 2.6.2. Local average of replicates (method $\ell_R$)

374 This estimator consists in replacing single indices by neighborhoods in (6).

For $b = 1, \dots, B$, let $\mathcal{V}_{j_1}^{(b)}$ and $\mathcal{V}_{j_2}^{(b)}$ (resp. $\mathcal{V}_{j_1'}^{(b)}$ and $\mathcal{V}_{j_2'}^{(b)}$) be two $\nu$-neighborhoods in $\mathcal{R}_j$ (resp. $\mathcal{R}_{j'}$) such that for any $i_1^{(b)} \in \mathcal{V}_{j_1}^{(b)}$, $i_2^{(b)} \in \mathcal{V}_{j_2}^{(b)}$, $|i_1^{(b)} - i_2^{(b)}| = \delta \ge p$ (resp. $|i_1'^{(b)} - i_2'^{(b)}| = \delta \ge p$ for any $i_1'^{(b)} \in \mathcal{V}_{j_1}^{(b)}$, $i_2'^{(b)} \in \mathcal{V}_{j_2}^{(b)}$). The local average of replicates based estimator is defined by

$$ \hat{r}_{jj'}^{\ell_R} = \frac{1}{B} \sum_{b=1}^{B} \frac{\frac{1}{4} \sum_{\alpha, \beta=1}^{2} \widehat{\operatorname{Cor}}(\bar{\mathbf{Y}}_{\mathcal{V}_{j_\alpha}^{(b)}}, \bar{\mathbf{Y}}_{\mathcal{V}_{j_\beta}^{(b)}})}{\sqrt{|\widehat{\operatorname{Cor}}(\bar{\mathbf{Y}}_{\mathcal{V}_{j_1}^{(b)}}, \bar{\mathbf{Y}}_{\mathcal{V}_{j_2}^{(b)}}) \widehat{\operatorname{Cor}}(\bar{\mathbf{Y}}_{\mathcal{V}_{j_1}^{(b)}}, \bar{\mathbf{Y}}_{\mathcal{V}_{j_2}^{(b)}})|}}. \quad (11) $$

### 2.6.3. Local averages and use of disconnected regions (method $\ell_D$)

380 We use in particular notation introduced in Sections 2.6.2 and 2.5.1 to propose the estimator $\hat{r}_{jj'}^{\ell_D}$ given by

$$ \hat{r}_{jj'}^{\ell_D} = \frac{1}{B} \sum_{b=1}^{B} \widetilde{\operatorname{Cor}}(\bar{\mathbf{Y}}_{\mathcal{V}_{j}^{(b)}}, \bar{\mathbf{Y}}_{\mathcal{V}_{j'}^{(b)}}; \bar{\mathbf{Y}}_{\mathcal{V}_{k}^{(b)}}, \bar{\mathbf{Y}}_{\mathcal{V}_{k'}^{(b)}}). \quad (12) $$

### 2.6.4. Replicates, local averages and use of a priori disconnected regions (method $\ell_{RD}$)

384 This estimator is a local version of $\hat{r}_{jj'}^{\ell_{RD}}$ and is defined by

$$ \hat{r}_{jj'}^{\ell_{RD}} = \frac{1}{B} \sum_{b=1}^{B} \frac{\frac{1}{4} \sum_{\alpha, \beta=1}^{2} \widetilde{\operatorname{Cor}}(\bar{\mathbf{Y}}_{\mathcal{V}_{j_\alpha}^{(b)}}, \bar{\mathbf{Y}}_{\mathcal{V}_{j_\beta}^{(b)}}; \bar{\mathbf{Y}}_{\mathcal{V}_{k}^{(b)}}, \bar{\mathbf{Y}}_{\mathcal{V}_{k'}^{(b)}})}{\sqrt{\widetilde{\operatorname{Cor}}(\bar{\mathbf{Y}}_{\mathcal{V}_{j_1}^{(b)}}, \bar{\mathbf{Y}}_{\mathcal{V}_{j_2}^{(b)}}; \bar{\mathbf{Y}}_{\mathcal{V}_{k}^{(b)}}, \bar{\mathbf{Y}}_{\mathcal{V}_{k'}^{(b)}}) \widetilde{\operatorname{Cor}}(\bar{\mathbf{Y}}_{\mathcal{V}_{j_1}^{(b)}}, \bar{\mathbf{Y}}_{\mathcal{V}_{j_2}^{(b)}}; \bar{\mathbf{Y}}_{\mathcal{V}_{k}^{(b)}}, \bar{\mathbf{Y}}_{\mathcal{V}_{k'}^{(b)}})}}. \quad (13) $$
## 2.7. Summary of estimators

We have formalised 9 estimators for inter-region correlation in fMRI, 6 of which are novel to the best of our knowledge. They vary in terms of their theoretical sensitivity to three factors: differences in region sizes and region intra-correlations ($\bar{\rho}_{R_j} \ll 1$), local noise ($\sigma_\epsilon$), and global noise ($\sigma_e$). Table 1 summarises estimators properties qualitatively using — for estimators that are sensitive to these factors, + for estimators that are insensitive, and ± for those that are in-between. The CA, $\ell_{CA}$, AC, and R estimators are sketched in Figure 1 and the D estimator is illustrated in Figure 2.

As an example, let us interpret the properties of CA shown in Table 1 in terms of these factors. First, we observe that the limit of $\hat{r}_{jj'}^{CA}$ strongly depends on the region size. Indeed, even in absence of noise this limit is $r_{jj'}/\sqrt{\bar{\rho}_{R_j}\bar{\rho}_{R_{j'}}}$, which can be quite far from $r_{jj'}$ especially for very large regions (so the estimator is sensitive to local noise and denoted — in the corresponding column). Now imagine that $\bar{\rho}_{R_j}\bar{\rho}_{R_{j'}} = 1$ and that $\sigma_e^2 = 0$ then the limit becomes $r_{jj'}/\sqrt{(1+\sigma_{\epsilon,j}^2\bar{\eta}_{R_j})(1+\sigma_{\epsilon,j'}^2\bar{\eta}_{R_{j'}})}$. Since it is expected that $\bar{\eta}_E$ is small (see (C.3)), especially for large sets $E$, this limit should be quite close to $r_{jj'}$ in this situation (+). Finally, if $\bar{\rho}_{R_j}\bar{\rho}_{R_{j'}} = 1$ and $\sigma_\epsilon^2 = 0$, and assume for simplicity that $\sigma_j = \sigma_{j'} = 1$, then $\hat{r}_{jj'}^{CA}$ would converge towards $(r_{jj'} + \sigma_e^2)/(1+\sigma_e^2)$ which can significantly deviate from $r_{jj'}$ when the global noise is strong (—).

This does not describe at all finite sample properties of the different estimators. Obviously, we could be tempted to always use the last two estimators (methods RD ad $\ell_{RD}$) which seem to be the most robust to additional noises. However, these last estimators will be less robust to small sample size. We
410 propose to investigate these finite sample properties in a simulation study
(Section 3.1) and real datasets (Sections 3.2 and 3.3 ).

412 We note that evaluating asymptotic variances of the different estimators
413 would add too much notation, assumptions and technicalities, and is left for
414 future work.
<table><caption>Table 1: Expected limits and properties for existing and novel estimators of inter-correlation $r_{jj'}$, under the model (1). We refer the reader to Section 2.3 for details on notation. In particular $\sigma_{e,j}^2$, $\sigma_{e,j}^2$, and $\sigma_{e,jj'}^2$ are given by (3) while $\bar{\rho}_E$, $\bar{\eta}_E$ and $\bar{\rho}_{E,E'}$ for two sets of indices $E, E'$ are given by (2) and (C.7). Sensitivity of estimators to three factors are reported: differences in region sizes and region intra-correlations ($\bar{\rho}_{R_j} \ll 1$), local noise ($\sigma_e$), and global noise ($\sigma_e$). Estimators that are sensitive to these factors are denoted $+$ and those in-between are denoted $\pm$.</caption><thead><tr><th rowspan="2">Estimator</th><th rowspan="2">Limit</th><th rowspan="2">r<sub>jj'</sub><sup>•</sup></th><th colspan="3">Sensitivity to</th></tr><tr><th>r<sub>jj'</sub> &plusmn; &sigma;<sup>2</sup><sub>e,jj'</sub></th><th>&rho;<sub>R<sub>j</sub></sub> &ll; 1</th><th>&sigma;<sub>e</sub></th><th>&sigma;<sub>e</sub></th></tr></thead><tbody><tr><td>&#x0102;<sup>CA</sup> (see (4))</td><td><sup>r<sub>jj'</sub> + &sigma;<sup>2</sup><sub>e,jj'</sub></sup>&frasl;<sub>&radic;(r<sub>R<sub>j</sub></sub> + &sigma;<sup>2</sup><sub>e,j</sub> &eta;&#x0102;<sub>R<sub>j</sub></sub> + &sigma;<sup>2</sup><sub>e,j</sub>)(r&#x0102;<sub>R<sub>j</sub></sub> + &sigma;<sup>2</sup><sub>e,jj'</sub> + &sigma;<sup>2</sup><sub>e,jj'</sub>)</sup></td><td>-</td><td>+</td><td>-</td></tr><tr><td>&#x0102;<sup>AC</sup> (see (5))</td><td><sup>r<sub>jj'</sub> + &sigma;<sup>2</sup><sub>e,jj'</sub></sup>&frasl;<sub>&radic;(1 + &sigma;<sup>2</sup><sub>e,j</sub> + &sigma;<sup>2</sup><sub>e,j</sub>)(1 + &sigma;<sup>2</sup><sub>e,jj'</sub> + &sigma;<sup>2</sup><sub>e,jj'</sub>)</sup></td><td>+</td><td>-</td><td>-</td></tr><tr><td>&#x0102;<sup>CA</sup> (see (10))</td><td><sup>r<sub>jj'</sub> + &sigma;<sup>2</sup><sub>e,jj'</sub></sup>&frasl;<sub>&radic;((r&#x0102;<sub>V</sub> + &sigma;<sup>2</sup><sub>e,j</sub> &eta;&#x0102;<sub>V</sub> + &sigma;<sup>2</sup><sub>e,j</sub>)(r&#x0102;<sub>V</sub> + &sigma;<sup>2</sup><sub>e,jj'</sub> &eta;&#x0102;<sub>V</sub> + &sigma;<sup>2</sup><sub>e,jj'</sub>))</sup></td><td>+</td><td>&plusmn;</td><td>-</td></tr><tr><td>&#x0102;<sup>R</sup> (see (6))</td><td><sup>r<sub>jj'</sub> + &sigma;<sup>2</sup><sub>e,jj'</sub></sup>&frasl;<sub>&radic;(|(r<sub>&delta;</sub> + &sigma;<sup>2</sup><sub>e,j</sub>)(r<sub>&delta;</sub> + &sigma;<sup>2</sup><sub>e,jj'</sub>)|)</sup></td><td>+</td><td>+</td><td>-</td></tr><tr><td>&#x0102;<sup>R</sup> (see (11))</td><td><sup>r<sub>jj'</sub> + &sigma;<sup>2</sup><sub>e,jj'</sub></sup>&frasl;<sub>&radic;(|(r&#x0102;<sub>V</sub>, v', &delta; + &sigma;<sup>2</sup><sub>e,j</sub>)(r&#x0102;<sub>V</sub>, v', &delta; + &sigma;<sup>2</sup><sub>e,jj'</sub>)|)</sup></td><td>+</td><td>+</td><td>-</td></tr><tr><td>&#x0102;<sup>D</sup> (see (7))</td><td><sup>r<sub>jj'</sub></sup>&frasl;<sub>&radic;(1 + &sigma;<sup>2</sup><sub>e,j</sub>)(1 + &sigma;<sup>2</sup><sub>e,jj'</sub>)</sup></td><td>+</td><td>-</td><td>+</td></tr><tr><td>&#x0102;<sup>RD</sup> (see (12))</td><td><sup>r<sub>jj'</sub></sup>&frasl;<sub>&radic;((r&#x0102;<sub>V</sub> + &sigma;<sup>2</sup><sub>e,j</sub> &eta;&#x0102;<sub>V</sub>)(r&#x0102;<sub>V</sub> + &sigma;<sup>2</sup><sub>e,jj'</sub> &eta;&#x0102;<sub>V</sub>))</sup></td><td>+</td><td>&plusmn;</td><td>+</td></tr><tr><td>&#x0102;<sup>RD</sup> (see (9))</td><td><sup>1</sup>&frasl;<sub>|r<sub>&delta;</sub>|</sub></td><td>+</td><td>+</td><td>+</td></tr><tr><td>&#x0102;<sup>RD</sup> (see (13))</td><td><sup>1</sup>&frasl;<sub>|(r&#x0102;<sub>V</sub>, v', &delta;)|</sub></td><td>+</td><td>+</td><td>+</td></tr></tbody></table>
### 3. Description of simulated and real datasets

We employed three distinct datasets to assess the performance of our estimators. These datasets encompassed a simulated dataset, a dataset involving rats that comprised both deceased and living animals, and a dataset from a healthy human subject, which included test-retest data.

#### 3.1. Simulated data

The paper being focused on pairwise spatial (auto)correlation estimation, it is sufficient to investigate the finite sample properties of our estimators for just two regions, say $\mathcal{R}_j$ and $\mathcal{R}_{j'}$ (whose sizes are set here to 20 and 40 voxels, respectively). Also, to save time and memory, we restrict ourselves, w.l.o.g., to one-dimensional regions ($d = 1$, regions are simply intervals so they are simply made of 'voxels' along a line). For the estimators based on differences (methods D, $\ell_D$, $\ell_{RD}$), we consider two extra regions, say $\mathcal{R}_k$ and $\mathcal{R}_{k'}$, that are disconnected (i.e., $r_{jk} = r_{jk'} = r_{kk'} = r_{j'k} = r_{j'k'} = 0$). We consider two scenarios: the “relatively strong inter-correlation case” ($r_{jj'} = 0.6$) and the “no inter-correlation case” ($r_{jj'} = 0$). The intra-correlation for any given region is modelled (alike within any region $\mathcal{R}_j, \mathcal{R}_{j'}, \mathcal{R}_k, \mathcal{R}_{k'}$) using the following spatial model

$$ \rho_{ii'} = 1 - (1 - \nu) \left( 1 - \frac{|i - i'|}{40} \right) \qquad (14) $$

which only depends on the distance $|i - i'|$ between two voxels, say $i$ and $i'$, belonging to the same region. We selected either the value $\nu = 0.8$ or $\nu = 0$, and designated accordingly the region as strongly or weakly intra-correlated. Hence, when $\nu = 0$, the two voxels furthest apart in the region
of size 40 are uncorrelated, i.e., $\rho_{ii'} = 0$ when $|i - i'| = 40$. When $\nu = 0.8$, they are highly correlated with $\rho_{ii'} = 0.8$. Figure 3 represents the correlation structure of (14), as well as the quantity $\bar{\rho}_E$ for sets of indices $E \subseteq \{1, \dots, 40\}$ with increasing size $\#E$, for both intra-correlation models. We generated independently 500 series of length $T = 1000$ according to model (1), and used 500 replicates for the method R and also 500 (Monte-Carlo or bootstrap) replications of choices of neighborhoods for methods $\ell_{CA}$, $\ell_{R}$, $\ell_{D}$, $\ell_{RD}$ (where we set the length of the neighborhood to 3). The local noises $\varepsilon_i(t)$ and $\varepsilon_{i'}(t)$ are assumed to be uncorrelated, so for some set of indices $E$, $\bar{\eta}_E = \sigma_e/\#E$. This is also represented in Figure 3.

Finally, we chose two values for the variance of the global noise, $\sigma_e^2 = 0$ and $\sigma_\varepsilon^2 = 0.1$, and two values for the variance of the local noise, $\sigma_\varepsilon^2 = 0$ and $\sigma_\varepsilon^2 = 0.1$. Results consist of 500 estimates for 9 methods, 2 intra-correlation models, 2 values for $\sigma_e^2$ and 2 values for $\sigma_\varepsilon^2$, that is 16 different scenarios (involving each time the 9 methods). They are presented and discussed in Section 6.1.
<figure><img src="image_4.png" /><figcaption>Figure 3: Simulation setup and results. (A) The two simulated one-dimensional regions (one with 40 “voxels”, the other with 20 “voxels”, shown as an inset) and their intra-correlation structure. Intra-correlation $ρ_{ii'}$ is given by (14) (with $ν = 0.8$ for the strong intra-correlation and $ν = 0$ for the weak intra-correlation) and decays with distance. (B) Intra-correlation (vertical axis) as a function of the size of region $E$ (horizontal axis). From top to bottom: $\bar{\rho}_E$ (orange): average intra-correlation of signal in the strong case; $\bar{\rho}_E$ (blue): average intra-correlation of signal in the weak case; $\bar{\eta}_E$ (black): average intra-correlation of noise. Average noise intra-correlation decays sharply with region size.</figcaption></figure>

## 3.2. Rats data

Using a 9.4T machine (Paravision 6.0.1, Bruker, Ettlingen, Germany), fMRI data were acquired for both dead and alive rats in Pawela et al. (2008). Twenty-five rats were scanned and identified in 4 different groups: DEAD, ETO-L (Etomidate), ISO-W (Isoflurane) and MED-L (Medetomidine). The first group contains dead rats and the three last groups correspond to different anesthetics. In this paper, we show results with data from three rats, one dead and two alive with different anesthetics (ETO-L, ISO-W).

The duration of the scanning was 30 minutes, using single-shot echo-planar imaging with TR / TE = 500 / 20 ms, so that 3600 time points were
available at the end of experiment. The resolution was 0.47 × 0.47 × 1.00 mm, slice gap 0.1 mm, 9 slices. After preprocessing as explained in Becq et al. (2020b), 51 brain regions for each rat were extracted using an in-house atlas. Sufficiently large regions are needed to be able to use the R estimator. We hence discarded regions that contained fewer than 40 voxels, and were left with 18 brain regions: The anterior cingulate cortex (ACC), bilateral Insular cortex (Ins\_r and Ins\_l), bilateral primary motor cortex (M1\_r and M1\_l), bilateral somatosensory 1 (S1\_r and S1\_l), bilateral somatosensory 1 barrel field (S1BF\_r and S1BF\_l), bilateral auditory cortex (AU\_r and AU\_l), bilateral caudate-putamen (striatum) (CPu\_r and CPu\_l), bilateral thalamus (Th\_r and Th\_l), bilateral basal forebrain region (BF\_r and BF\_l), bilateral hippocampus (HIP\_r and HIP\_l).

Voxel time series were waveled-filtered using Daubechies orthonormal compactly supported wavelet of length 8.

### 3.3. Human Connectome Project data

We also evaluated our estimators on a subset of the Human Connectome Project (HCP) Young Adult 1200 Subjects release, WU-Minn Consortium pre-processed (Glasser et al., 2013) (connectome db data package *Resting State fMRI 1/2 Preprocessed*). We selected 100 subjects with two rs-fMRI acquisitions on different days. The TR was 720 ms and the duration of acquisition was 14 min and 24s.

The preprocessed fMRI data was segmented into 89 regions with SPM *New Segment* using a modified AAL template: merging some of the regions, reducing the parcellation to 89 regions. Merged regions are: frontal medial orbital and rectus (one region for left and one for right hemisphere); oc-
cipital superior, middle and inferior (one region for left and one for right hemisphere); temporal pole superior and medial (one region for left and one for right hemisphere); the cerebral crus (one region for left and one for right hemisphere); areas III, IV, V and VI of cerebellum (one region for left and one for right hemisphere); areas VII, VIII, IX, X of cerebellum (one region for left and one for right hemisphere) and finally, the vermis (one single region for both hemispheres). Other details are available in Termenon et al. (2016).

Voxel time series were wavelet filtered using Daubechies orthonormal compactly supported wavelet of length 8.

# 4. Evaluation and metrics

First, on *simulated data*, we qualitatively inspected the bias and variance of the distribution of correlation values with respect to known ground truth for various levels of global and local noise.

Then, using *rat data*, we performed a face validity analysis of the estimators, with the premise that dead rats should show no functional connectivity (the correlation distribution should be centered at zero). In order to quantify the differences between correlation values obtained for dead and live rats, we computed the Wasserstein distance between the correlation distributions of each anesthetized rat in comparison to that of a dead rat. A low value of the Wasserstein distance indicates that correlations values of live and dead rats are comparable and counts negatively in the evaluation of an estimator.

To evaluate the repeatability of the proposed estimators on the rat dataset, we split the time series in two equal parts. We computed the correlations on each part using the whole range of proposed estimators, and computed
the Concordance Correlation Coefficient (CCC) (Lin, 1989) between splits to provide a scaled measure of agreement, where 1 is perfect agreement and 0 is no agreement. A preferable estimator should be more repeatable and have higher CCC.

To quantify the similarity of connectivity graphs between estimators, we computed the number of common edges between graphs obtained from each estimator. To this end we used a sparsity threshold equal to 20% of the total number of edges (i.e., 27 edges in our case with 18 regions).

For *human data*, we used rs-fMRI sessions from different days, in order to evaluate reproducibility of correlation coefficients. This was analysed using CCC, again with a preferable estimator being more reproducible and having higher CCC.

We also evaluated the reproducibility of graph metrics between sessions. To this end we used a sparsity threshold equal to 20% of the total number of edges, keeping only edges with the highest correlation (i.e., 783 edges in our case with 89 regions), and binarized the edges. In order to compute graph metrics, we forced the graph to be connected by applying a minimum spanning tree (Alexander-Bloch et al., 2010). Then we computed classical graph metrics: betweenness centrality, transitivity, global and local efficiencies using package iGraph. Reproducibility was evaluated using the CCC.

We also summarized the differences of connectivity graphs between estimators, by computing the number of common edges between graphs obtained from CA and *l*CA using thresholding at the 20th percentile (i.e., 783 edges with 89 regions), and visualized the difference qualitatively by taking absolute values of correlation values for each estimator, rank-transforming, and
computing median difference in ranks across all subjects,

Additionally, we also evaluated discriminative power of the various estimators via three metrics: inter vs intra-subject graph distance, a non-parametric test of the same, and identification rate using functional connectome fingerprinting (Finn et al., 2015). A desirable estimator should provide estimates that preserve inter-individual differences.

We defined the intra-subject distance as the distance between the graph representing the first rs-fMRI session and the graph representing the second rs-fMRI session. The inter-subject distance was computed between each subject's first session and all other subject's first sessions. Separation between the intra-subject distances and the inter-subject distances was quantified by mean and standard deviation of the distributions, and by a Wilcoxon rank-sum test on multiple random splits of subject data, to avoid having multiple measurements of the same subjects. Here, we repeated 10 times the following procedure for each estimator of interest: first, split the subjects into two disjoint sets - one used to compute intra-distances (50 subjects), and one to compute inter-distances (50 subjects). Within the inter-distances set, 25 subject pairs were formed randomly. We tested the null hypothesis of no difference between inter- and intra-distances, against the alternative hypothesis that intra-subject distance < inter-subject distance, based on the assumption that subjects are more similar to themselves than to other subjects. Given the relatively narrow age range of our sample of HCP subjects (all 22-35 except one 36+), and given that our goal was to compare estimators using fixed splits, we did not adjust for covariates or match samples across splits. We used a one-sided Wilcoxon rank-sum test, yielding a W statistic and a *p*-
value for each of the 10 runs. We then computed the average W value across
runs, as well as the harmonic mean *p*-value (Wilson, 2019) across runs, a pro-
cedure with strong family-wise error rate (FWER) control even for positively
dependent tests.

To compute identification rate, functional connectome fingerprinting represent each subject's graph *g* as a vectorized version **a** of the upper-triangular (or lower-triangular) part of the full inter-region correlation matrix (whose entries are *r*<sub>ii'</sub>), and computes the fingerprinting distance between graphs as *d*(*g*<sub>1</sub>, *g*<sub>2</sub>) = 1 - Cor(<b>a</b><sub>1</sub>, <b>a</b><sub>2</sub>), where Cor denotes Pearson correlation. From the (intra, inter) fingerprinting distance distributions, the identification counts as correct if the intra-subject distance is lower than all inter-subject distances. This is equivalent to a top-1 recognition rate. We note there are many other possibilities to compute distances between such brain graphs (Richiardi et al., 2013; Ng et al., 2016; Dadi et al., 2019), including computing distances between graph embeddings, which could substantially alter results.

Finally, we evaluated the dependence on region size by computing Spear-
man correlations between atlas region size and the average of correlations in
which the region is involved (itself averaged across subjects). A preferable
estimator should minimize dependence to region size, and show lower Spear-
man correlation. We tested differences between estimators using a paired
t-test between these Spearman correlations.

5. Data and Code Availability

R and Python code implementing all estimators, to generate simulated
data, and to extract the time-series from the preprocessed HCP data is avail-
586 able at https://gitlab.inria.fr/q-func/ireco4fmri.

587 The pre-extracted, wavelet-filtered time series for the rat data are available at https://dx.doi.org/10.5281/zenodo.7254133. Human Connectome Project data is available at https://www.humanconnectome.org/

## 6. Results

### 6.1. Evaluation on simulated data

591 Simulation setup is described in Section 3.1. Figure 4 shows boxplots of estimates of $r_{jj'} = 0.6$ for all methods and different intra-correlation models, and different levels of local and global noise. In terms of bias, overall, the method $\ell_{RD}$ is the best, but it is also the one with the highest variance. In the strong intra-correlation case, and when $\sigma_e = 0$, all methods are almost unbiased. When $\sigma_e$ is increased to 0.1, the estimators $\ell_{AC}$, $\ell_{CA}$, $\ell_D$, and $\ell_D$ clearly lose this property. In the weak intra-correlation case, only the estimators $\ell_{AC}$, $\ell_R$ and $\ell_{RD}$ are unbiased, or close to, with $\ell_{RD}$ being the best overall for this criterion, while still being the more variable. Figure 5 shows boxplots of estimates of $r_{jj'} = 0$ for all methods and different intra-correlation models, and different levels of local and global noise. When $\sigma_e = 0$ all estimators are unbiased, both in the strong and weak intra-correlation case. This property remains true when $\sigma_e$ is increased to 0.1 only for the estimators $\ell_D$, $\ell_D$ and $\ell_{RD}$. Here again, the estimator $\ell_{RD}$ is the more variable. We can also notice that when $\sigma_e = 0.1$, the estimator $\ell_{AC}$ exhibits very good properties, while $\ell_{RD}$ is the worst.
<figure><img src="image_5.png" /><figcaption>Figure 4: Estimates of the inter-correlation parameter $r_{jj'} = 0.6$ between two regions, based on 500 simulation runs of the general model (1). Situations for two intra-correlation models and situations with no noise, local noise and/or global noise are considered. The true inter-correlation is depicted by a red dashed line.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5: Estimates of the inter-correlation parameter $r_{jj'} = 0$ between two regions, based on 500 simulation runs of the general model (1). Situations for two intra-correlation models and situations with no noise, local noise and/or global noise are considered. The true inter-correlation is depicted by a red dashed line.</figcaption></figure>

## 6.2. Evaluation on rat data

Figure 6(A) shows the correlation values obtained on rats for all pairs of brain regions, 153 in our case. For this data set, we know that for the dead
rat we are under the full null hypothesis as no legitimate functional activity
should be detected. Thus the estimated correlations should be close to zero.
This is the case for estimators AC, R, $l_{CA}$, D and $l_{D}$. However, the other
estimators showcase a clear bias towards positive values. The method CA
namely yields unexpectedly high values of correlations. These correlations
correspond to regions that are close together (Becq et al., 2020a). In order
to validate these methods, we also apply our estimators to live rats. The
results of two live rats is shown in Figure 6 (A, right). As expected, due to
the local noise, the methods AC and D do not provide satisfactory results as
the correlation values are very close to zero. One of the best method in this
case is $l_{CA}$, where sufficient non-zero correlations are obtained. Wasserstein
distance computations (Figure 6 (C)) show that AC, D, and RD have the low-
est Wasserstein distance values, indicating that the correlation distribution
of the live rats resemble that of a dead rat.

Figure 6(B) shows Concordance Correlation Coefficient results. Consistent with the all-noise nature of the data, the dead rat exhibited very low repeatability, with $l_{CA}$ providing the highest at 0.22. On the live Eto-L rat, estimators had approximately the same repeatability, with RD showing the lowest CCC at 0.62 and AC tied with $l_{CA}$ for highest at 0.87. For the Iso-W rat, $l_{R}$ had the lowest CCC at 0.46, CA the second lowest at 0.54, and $l_{D}$ the highest at 0.73.

Combining all of these results, $l_{CA}$, R and $l_{D}$ hence seem to be the most
adequate correlation estimators. However, as shown in formula (G.4), the
estimator $l_{D}$ is difficult to implement. Indeed, it requires the definition of
two other regions uncorrelated with the main brain regions of the parcellation
<figure><img src="image_7.png" /><figcaption>Figure 6: Rat data results. **A.** Empirical distribution of the correlation estimators for all pairs of brain regions for a dead and two anesthetized rats, for all proposed estimators. In the dead rat, the correlation of averages (CA) estimator is providing high values where null correlations should be observed. For the live rat the average of correlation estimator (AC) is providing very low values where non null correlations should be observed. **B.** The Concordance Correlation Coefficient (CCC) for the repeatability of the different estimators for all rats, calculated between the first and second half of the BOLD time series. Higher CCC corresponds to a more repeatable estimator. **C.** Wasserstein distances between the correlation distribution of each anesthetized rat and that of the dead rat, for all estimators. AC, D, RD have a very low distance, indicating that correlation values are similar between dead and live rats for these estimators.</figcaption></figure>
and uncorrelated with themselves. Moreover, R cannot be estimated when regions are too small, which is often the case in rat data. From now on, we will hence focus on estimator $\ell_{CA}$.

We then quantified the edges in common between the networks obtained via the two estimators CA (which is currently the most widely used estimator) and $\ell_{CA}$. For the dead rat, 67% of edges are in common between the two estimators. Additionally, 60% and 77% of edges are similar for the live rats. https://www.overleaf.com/project/619e15e123dc6c1b1e091576

### 6.3. Evaluation on human data

Based on our findings on the rats datasets, we evaluate the performances of the three estimators CA (most common estimator, highest dead-live rat distance), AC (low dead-live rat distance) and $\ell_{CA}$ (high dead-live rat distance) for 100 subjects of the HCP dataset.

Figure 7(A) reports the correlation values among all pairs of regions for four randomly selected HCP subjects. Consistent with the rat results, the estimator CA yields the largest values of correlations, estimator AC yields very low values, while $\ell_{CA}$ values are different from zero, but smaller that CA values.

Reproducibility results for correlation estimates are shown in Figure 7 (B,C). The Concordance Correlation Coefficient was similar between estimators (average (sd) across 100 subjects for CA: 0.64 (0.13), AC: 0.66 (0.20), $\ell_{CA}$: 0.62 (0.17), R: 0.56 (0.17), $\ell_{R}$: 0.52 (0.14)), with variations in reproducibility reflecting inter-subject variability more than differences between estimators. For graph metrics reproducibility, we report only the results with betweenness in Figure 7 (D,E), since similar results are obtained with
<figure><img src="image_8.png" /><figcaption>Figure 7: Human data results. **A.** Empirical distribution of inter-regional correlations for three selected estimators for all pairs of brain regions for four human subjects. Each subject was scanned twice, on different days. **B. and C.** Concordance correlation coefficient (CCC) for the reproducibility of the inter-regional correlation values obtained by different estimators for all human subjects, computed between the two examinations. Higher CCC indicates a more repeatable estimator. All estimators have broadly similar reproducibility. **D. and E.** Reproducibility of a topological graph metric (betweenness). Again all estimators give broadly similar results, with slightly higher reproducibility for AC</figcaption></figure>
other metrics. Here, the methods differed more, with average (sd) across 661 100 subjects for CA: 0.29 (0.14), AC: 0.55 (0.17), $\ell_{CA}$: 0.4 (0.16), R: 0.37 662 (0.14), $\ell_{R}$: 0.26 (0.15). $\ell_{CA}$ had significantly lower CCC than AC (T=-6.8, 663 $p$ = 1.6e-10), However, $\ell_{CA}$ has significantly higher CCC than CA (T=5.1, 664 $p$ = 8e-07) and $\ell_{R}$ (T=6.3, $p$ = 2e-9). Finally $\ell_{CA}$ and R are not significantly 665 different (T=1.35, $p$ = 0.18). These differences are robust to the choice of 666 threshold (cf, Appendix Appendix G.5). 667

In the thresholded graphs, the percentage of edges in common between 668 estimators CA and $\ell_{CA}$ was on average equal to 70% for the one hundred 669 subjects used in this analysis for both sessions. Figure 8 shows median 670 differences between the estimators in brain space across the HCP subjects. 671

<figure><img src="image_9.png" /><figcaption>Figure 8: Largest differences between the CA and $\ell_{CA}$ estimators, median over 100 HCP subjects. Only the top 20% differences are shown. Inter-regional correlations are taken in absolute value and rank-transformed prior to computing differences (rank 1 for the strongest correlation, rank 2 for the second-strongest, and so on). Red indicates absolute correlations that are higher for the $\ell_{CA}$ than the CA estimator, while blue indicates the reverse. Node size is proportional to region size in the atlas. Estimator CA on average shows hyperconnectivity in occipital and generally dorsal posterior regions, and hypoconnectivity in frontal, temporal, and general ventral anterior regions.</figcaption></figure>

Looking at dependence on region size, the CA estimator showed signifi- 672
cantly more correlation with region size than the $\ell_{CA}$ estimator (average (sd) across 100 subjects 0.55 (0.10) vs 0.40 (0.09), T=27.2, $p = 1.1e^{-47}$.

In terms of discriminative power between subjects, for connectome fingerprinting, CA and $\ell_{CA}$ achieved the same performance (72% correct identification), while AC had slightly lower performance (68% correct identification). Group differences were also similar between estimators. Table 2 provides details.

<table><caption>Table 2: Discriminative power of estimators on the human dataset. intra: Within-subject average and standard deviation of graph distances between first and second imaging session across 100 subjects; inter: same for between-subjects, using only the first session. W: average one-sided Wilcoxon rank-sum test value on 10 random splits, with corresponding harmonic mean p-value</caption><thead><tr><th>estimator</th><th>intra (sd)</th><th>inter (sd)</th><th>W (p-value)</th><th>identification rate</th></tr></thead><tbody><tr><td>CA</td><td>0.29 (0.10)</td><td>0.49 (0.10)</td><td>-5.82 (p<sub>hmp</sub> = 1.5e-10)</td><td>72%</td></tr><tr><td>AC</td><td>0.19 (0.10)</td><td>0.32 (0.10)</td><td>-4.88 (p<sub>hmp</sub> = 3.6e-9)</td><td>69%</td></tr><tr><td>l<sub>CA</sub></td><td>0.26 (0.10)</td><td>0.43 (0.10)</td><td>-3.315.58 (p<sub>hmp</sub> = 1.2e-9)</td><td>72%</td></tr></tbody></table>

## 7. Discussion

In this paper we illustrate the effect of averaged data on estimators of correlation when two types of noises are present, local and global noise. The use of the classical correlation of averages is hindered by the presence of these noises in addition to the presence of intra-correlations. We proposed alternative estimators including correction terms to compensate the intra-correlations, local and global noises. The performance of these estimators
687 was evaluated on simulations, rats data, and human data, yielding several
688 observations.

7.1. *The correlation of averages estimator is highly biased*

690 The CA estimator tends to be highly biased, as illustrated on synthetic
data where the ground truth is known, but also compared to other estimators,
691 as shown on live rats and human data, where the mode of the distribution of
692 correlation values is systematically among the highest found. We hypothesise
693 that this is driven by a combination of low intra-correlation and large region
694 sizes, which further lowers intra-correlation. This can be seen from the es-
695 timator definition in Eq. 4. We also note that the *l*CA estimator effectively
696 reduces this influence of region size.
697

In addition, Figure 8 revealed a systematic spatial bias between the CA
and ℓCA estimator, exhibiting dorsal posterior hyper-connectivity for CA,
and corresponding ventral anterior hypo-connectivity. The figure also sug-
gests that the largest differences between the two estimators appear between
regions that are the largest, further highlighting the reduced dependency to
region size for the ℓCA estimator. The spatial distribution of these differences
suggests that caution is in order when examining large-scale resting-state net-
works derived from the CA estimator, as some apparent topological properties
of brain networks, such as modularity, could be driven in part by region size
and region intra-correlation. Indeed, in our experiments, thresholded graphs
differed in a large proportion of edges, both in rats (around 30%-50% edge
differences) and humans (around 30% edge differences). Thus, it is probable
that both edge-level and graph-level metrics obtained from the CA estima-
tor are biased due to their over- or under-estimation of actual functional
connectivity, in a spatially-dependent manner. For clinical applications, this phenomenon could either emphasize or reduce differences between patients and controls. Since we have no ground truth available for in-vivo functional connectivity, in practical situations, we therefore recommend that results obtained with the CA estimator be re-run at least with the ℓ<sub>CA</sub> estimator as a sensitivity analysis. The computational cost is not excessive, and differences in results could indicate that estimator-induced bias was at play.

7.2. local noise and intra-correlation link to long-range correlation

In this paper, we explain the bias observed in CA estimator by introduc-
ing hypotheses on both intra-correlation and noise. Indeed, previous studies
on regional homogeneity (Zang et al., 2004) showed relevant results on clas-
sification of pathologies based only on intra-regional properties. This was
confirmed by a recent work on classification of intra-correlation (Petersen
et al., 2016) using Wasserstein distances. Based on these findings, we hy-
pothesize that bias observed on inter-correlation is driven by intra-correlation
and noise. Our simple simulation model illustrates the effect of local noise
and intra-correlation. This is clearly displayed in Figures 4-5, where the box-
plots for the various estimators are plotted. However, it is important to note
that under local noise, in this framework with controlled intra-correlation,
estimator CA is relatively close to the exact value. This may be explained
by a trade-off in the denominator of the limit as expressed in Table 1. In
our simulation, we also observed that the CA estimators bias depends on
the intra-correlation and local noise. Indeed, high values of CA tends to be
observed when low values of intra-correlation are observed. These low values
of intra-correlation have already been mentioned in the study of dynamics of
neural networks (Deco et al., 2014) where local decorrelation was reported in
real datasets. In our paper, for the first time, we proved a statistical expla-
nation of the link between local decorrelations and long-range correlations
using aggregated time series.

The model chosen in this paper for intra-correlation and local noise was
driven by statistical motivations to be able to write explicit formulas for
the limit of the estimators. However, as observed in Jiang and Zuo (2016);
Deco et al. (2014), these hypotheses are realistic for resting-state fMRI data,
where local decorrelations are observed. These local decorrelations can come
from two factors: a low intra-correlation (as modeled by the choice of the
intra-correlation coefficients of the matrix ), or a strong local noise. The
stationarity assumption may be not adequate based on raw data. However,
as mentioned in Section 2.1, it becomes very reasonable after performing a
wavelet transform of each time series voxelwise. This preprocessing signifi-
cantly reduces non stationary artefacts.

7.3. *Repeatability and reproducibility*

Repeatability of correlation values in dead rats was very low for all es-
timators, consistent with the random nature of the data. For live rats, the
CCC ranged from 0.46 to 0.87 depending on specimen and estimator. For
humans, CA and *l*CA showed approximately the same reproducibility (0.63
average (0.2)), and AC was slightly superior (0.66 average (0.2)). But repro-
ducibility differences between estimators were much less pronounced than
reproducibility differences between individual subjects.

As a representative for the reproducibility of graph metrics, we inves-
tigated betweenness. Here, AC offered the highest reproducibility (average
(sd): 0.55 (0.17)) and $\ell_{CA}$ improved markedly over CA (0.4 (0.16) vs 0.29 (0.14)). This is contrast to another study that found no effect of aggregation method (region mean time series versus region median versus 1st eigenvariate of the region) on the reproducibility of graph metrics (Braun et al., 2012) (although in that study sessions were weeks apart).

### 7.4. Discriminability

Estimators CA, AC, $\ell_{CA}$ showed similar values for discriminability, with slightly lower identification rate and intra-subject to inter-subject distribution separation for AC than the two others, and slightly lower intra-inter separation for $\ell_{CA}$ than CA. This suggests that the improved robustness to region size and intra-correlation effects of $\ell_{CA}$ does not result in a sizeable impact on discriminative ability, although this warrants further evaluation.

### 7.5. Limitations

Our signal model, and therefore the derived estimators, is a trade-off between model realism and tractability of the analysis of estimator properties. This comes with important limitations.

First, assuming stationarity and additivity of the local noise fails to capture effects like system instability due to B0 inhomogeneity, RF power variations, or gradient fluctuations (Lazar, 2008; Greve et al., 2013; Liu, 2016). Independently of the model, note that effects such as drift are mitigated by using wavelet coefficient time series as we did in this study, and that such instabilities explain proportionally less of the noise variance than thermal noise at high field (Greve et al., 2011).
Second, motion effects, and in particular differential long-vs short-range effects on correlations (Van Dijk et al., 2010; Yan et al., 2013), were not studied, and their interplay with the spatial bias exhibited by estimator CA in Figure 8 was not examined.

Third, our new estimators come with the added burden of choosing hyper-parameters such as neighbourhood size. These are currently selected empirically, and no systematic sensitivity analysis has been performed. However, our proposed approach may be used to redefine the brain regions by grouping voxels with high intra-correlation. This would allow to define new brain regions using intra-correlation in addition to anatomical criterion.

Despite these limitations, we believe our empirical tests served to bridge the gap towards applicability, since our model yielded at least an estimator, $ℓ_{CA}$, with useful properties for use in neuroimaging - namely, reduced dependency to region size and low intra-correlation, and improved reproducibility of graph metrics.

# 8. Acknowledgements

This work was partly funded by French National Research Agency project ANR-20-NEUC-0003-02. Rat data was acquired at the IRMaGe MRI facility, partly funded by the French program Investissement d'Avenir run by the French National Research Agency, grant Infrastructure d'avenir en Biologie Santé ANR-11-INBS-0006. The work was also partly funded by the transnational Alliance Campus Rhodanien programme. Human data were provided by the Human Connectome Project, WU-Minn Consortium (Principal Investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657) funded by
the 16 NIH Institutes and Centers that support the NIH Blueprint for Neuroscience Research; and by the McDonnell Center for Systems Neuroscience at Washington University. We also thank Ayoub Cherkaoui and Ossama Ziri for coding the Python implementation of the estimators.

References

Achard, S., Coeurjolly, J.F., Marcillaud, R., Richiardi, J., 2011. fMRI functional connectivity estimators robust to region size bias, in: IEEE Workshop on Statistical Signal Processing, SSP2011, Nice, France. pp. 813–816.

Achard, S., Salvador, R., Whitcher, B., Suckling, J., Bullmore, E., 2006. A resilient, low-frequency, small-world human brain functional network with highly connected association cortical hubs. Journal of Neuroscience 26, 63–72.

Afyouni, S., Smith, S.M., Nichols, T.E., 2019. Effective degrees of freedom of the Pearson’s correlation coefficient under autocorrelation. NeuroImage 199, 609–625.

Alexander-Bloch, A.F., Gogtay, N., Meunier, D., Birn, R., Clasen, L., Lalonde, F., Lenroot, R., Giedd, J., Bullmore, E.T., 2010. Disrupted Modularity and Local Connectivity of Brain Functional Networks in Childhood-Onset Schizophrenia. Frontiers in Systems Neuroscience 4.

Alexander-Bloch, A.F., Shou, H., Liu, S., Satterthwaite, T.D., Glahn, D.C., Shinohara, R.T., Vandekar, S.N., Raznahan, A., 2018. On testing for spatial correspondence between maps of human brain structure and function. NeuroImage 178, 540–551.
Becq, G.G.J.P.C., 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, G.J.P., Habet, T., Collomb, N., Faucher, M., Delon-Martin, C., Coizet, V., Achard, S., Barbier, E.L., 2020b. Functional connectivity is preserved but reorganized across several anesthetic regimes. NeuroImage 219, 116945.

Bergholm, F., Adler, J., Parmryd, I., 2010. Analysis of bias in the apparent correlation coefficient between image pairs corrupted by severe noise. J Math Imaging Vis 37, 204-219.

Bolt, T., Nomi, J.S., Rubinov, M., Uddin, L.Q., 2017. Correspondence between evoked and intrinsic functional brain network configurations. Human Brain Mapping 38, 1992-2007. https://onlinelibrary.wiley.com/doi/pdf/10.1002/hbm.23500.

Braun, U., Plichta, M.M., Esslinger, C., Sauer, C., Haddad, L., Grimm, O., Mier, D., Mohnke, S., Heinz, A., Erk, S., Walter, H., Seiferth, N., Kirsch, P., Meyer-Lindenberg, A., 2012. Test-retest reliability of resting-state connectivity network characteristics using fMRI and graph theoretical measures. NeuroImage 59, 1404-1412.

Burt, J.B., Helmer, M., Shinn, M., Anticevic, A., Murray, J.D., 2020. Generative modeling of brain maps with spatial autocorrelation. NeuroImage 220, 117038.
Büchel, C., Friston, K.J., 1997. Modulation of connectivity in visual pathways by attention: cortical interactions evaluated with structural equation modelling and fMRI. Cerebral Cortex 7, 768-778.

Caballero-Gaudes, C., Reynolds, R.C., 2017. Methods for cleaning the BOLD fMRI signal. NeuroImage 154, 128-149.

Cao, H., McEwen, S.C., Forsyth, J.K., Gee, D.G., Bearden, C.E., Addington, J., Goodyear, B., Cadenhead, K.S., Mirzakhanian, H., Cornblatt, B.A., Carrión, R.E., Mathalon, D.H., McGlashan, T.H., Perkins, D.O., Belger, A., Seidman, L.J., Thermenos, H., Tsuang, M.T., van Erp, T.G.M., Walker, E.F., Hamann, S., Anticevic, A., Woods, S.W., Cannon, T.D., 2019. Toward Leveraging Human Connectomic Data in Large Consortia: Generalizability of fMRI-Based Brain Graphs Across Sites, Sessions, and Paradigms. Cerebral Cortex 29, 1263-1279.

Cao, J., Worsley, K., 1999. The geometry of correlation fields with an application to functional connectivity of the brain. The Annals of Applied Probability 9, 1021-1057. Publisher: Institute of Mathematical Statistics.

Castellanos, F.X., Di Martino, A., Craddock, R.C., Mehta, A.D., Milham, M.P., 2013. Clinical applications of the functional connectome. NeuroImage 80, 527-540.

Castrillon, J.G., Ahmadi, A., Navab, N., Richiardi, J., 2014. Learning with multi-site fMRI graph data, in: 2014 48th Asilomar Conference on Signals, Systems and Computers, pp. 608-612. ISSN: 1058-6393.
876 Chen, J.E., Glover, G.H., 2015. BOLD fractional contribution to resting-state functional connectivity above 0.1 Hz. *NeuroImage* 107, 207–218.

878 Clifford, P., Richardson, S., Hemon, D., 1989. Assessing the significance of the correlation between two spatial processes. *Biometrics* 45, 123–134.

880 Cordes, D., Haughton, V.M., Arfanakis, K., Carew, J.D., Turski, P.A., Moritz, C.H., Quigley, M.A., Meyerand, M.E., 2001. Frequencies Contributing to Functional Connectivity in the Cerebral Cortex in “Resting-state” Data. *American Journal of Neuroradiology* 22, 1326–1333. Publisher: *American Journal of Neuroradiology Section: BRAIN*.

885 Dadi, K., Rahim, M., Abraham, A., Chyzhyk, D., Milham, M., Thirion, B., Varoquaux, G., 2019. Benchmarking functional connectome-based predictive models for resting-state fMRI. *NeuroImage* 192, 115–134.

888 Deco, G., Ponce-Alvarez, A., Hagmann, P., Romani, G.L., Mantini, D., Corbetta, M., 2014. How Local Excitation-Inhibition Ratio Impacts the Whole Brain Dynamics. *Journal of Neuroscience* 34, 7886–7898.

891 Descombes, X., Kruggel, F., Von Cramon, D., 1998. fMRI Signal Restoration Using a Spatio-Temporal Markov Random Field Preserving Transitions. *NeuroImage* 8, 340–349.

894 Eickhoff, S.B., Thirion, B., Varoquaux, G., Bzdok, D., 2015. Connectivity-based parcellation: Critique and implications. *Human Brain Mapping* 36, 4771–4792.

897 Fan, F., Liao, X., Lei, T., Zhao, T., Xia, M., Men, W., Wang, Y., Hu, M.,
898

Liu, J., Qin, S., Tan, S., Gao, J.H., Dong, Q., Tao, S., He, Y., 2021. Development of the default-mode network during childhood and adolescence: A longitudinal resting-state fMRI study. *NeuroImage* 226, 117581.

899

900

901

Figueroa-Jimenez, M.D., Cañete-Massé, C., Carbó-Carreté, M., Zarabozo-Hurtado, D., Peró-Cebollero, M., Salazar-Estrada, J.G., Guàrdia-Olmos, J., 2021. Resting-state default mode network connectivity in young individuals with down syndrome. *Brain and Behavior* 11, e01905. https://onlinelibrary.wiley.com/doi/pdf/10.1002/brb3.1905.

902

903

904

905

906

907

908

909

Finn, E.S., Shen, X., Scheinost, D., Rosenberg, M.D., Huang, J., Chun, M.M., Papademetris, X., Constable, R.T., 2015. Functional connectome fingerprinting: identifying individuals using patterns of brain connectivity. *Nature Neuroscience* 18, 1664–1671.

910

911

912

Fornito, A., Bullmore, E.T., Zalesky, A., 2017. Opportunities and Challenges for Psychiatry in the Connectomic Era. *Biological Psychiatry: Cognitive Neuroscience and Neuroimaging* 2, 9–19.

913

914

Fornito, A., Zalesky, A., Breakspear, M., 2015. The connectomics of brain disorders. *Nature Reviews Neuroscience* 16, 159–172.

915

916

Gaetan, C., Guyon, X., et al., 2010. Spatial statistics and modeling. volume 90. Springer.

917

918

919

920

Glasser, M.F., Sotiropoulos, S.N., Wilson, J.A., Coalson, T.S., Fischl, B., Andersson, J.L., Xu, J., Jbabdi, S., Webster, M., Polimeni, J.R., Van EsSEN, D.C., Jenkinson, M., 2013. The minimal preprocessing pipelines for the Human Connectome Project. *NeuroImage* 80, 105–124.
Greve, D.N., Brown, G.G., Mueller, B.A., Glover, G., Liu, T.T., Function 921 Biomedical Research Network, 2013. A Survey of the Sources of Noise in 922 fMRI. Psychometrika 78, 396–416. 923

Greve, D.N., Mueller, B.A., Liu, T., Turner, J.A., Voyvodic, J., Yetter, E., 924 Diaz, M., McCarthy, G., Wallace, S., Roach, B.J., Ford, J.M., Mathalon, 925 D.H., Calhoun, V.D., Wible, C.G., Brown, G.G., Potkin, S.G., Glover, 926 G., 2011. A novel method for quantifying scanner instability in fMRI: 927 Quantifying Scanner Instability in fMRI. Magnetic Resonance in Medicine 928 65, 1053–1061. 929

Gunst, R.F., 1995. Estimating spatial correlations from spatial-temporal 930 meteorological data. Journal of climate 8, 2454–69. 931

Hartvig, N.V., Jensen, J.L., 2000. Spatial mixture modeling of fMRI data. 932 Human Brain Mapping 11, 233–248. 933

Jiang, L., Zuo, X.N., 2016. Regional Homogeneity: A Multimodal, Multiscale 934 Neuroimaging Marker of the Human Connectome. The Neuroscientist 22, 935 486–505. 936

Jo, H.J., Saad, Z.S., Simmons, W.K., Milbury, L.A., Cox, R.W., 2010. Map- 937 ping sources of correlation in resting state fMRI, with artifact detection 938 and removal. NeuroImage 52, 571–582. 939

Köhler, S., McIntosh, A.R., Moscovitch, M., Winocur, G., 1998. Functional 940 interactions between the medial temporal lobes and posterior neocortex 941 related to episodic memory retrieval. Cerebral Cortex 8, 451–461. 942
943

Lazar, N.A., 2008. Noise and Data Preprocessing, in: The Statistical Analysis of Functional MRI Data. Springer New York, New York, NY, pp. 37–51.
ISSN: 1431-8776 Series Title: Statistics for Biology and Health.

944

945

946

Lin, L.I.K., 1989. A Concordance Correlation Coefficient to Evaluate Re-
producibility. Biometrics 45, 255–268. Publisher: [Wiley, International
Biometric Society].

947

948

949

Liu, T.T., 2016. Noise contributions to the fMRI signal: An overview. NeuroImage 143, 141–151.

950

951

Markello, R.D., Misic, B., 2021. Comparing spatial null models for brain maps. NeuroImage 236, 118052.

952

953

Moulines, E., Roueff, F., Taqqu, M.S., 2007. On the spectral density of the wavelet coefficients of long-memory time series with application to the log-regression estimation of the memory parameter. Journal of Time Series Analysis 28, 155–187.

954

955

956

957

Murphy, K., Birn, R.M., Bandettini, P.A., 2013. Resting-state fMRI con-
stands and cleanup. NeuroImage 80, 349–359.

958

959

Ng, B., Varoquaux, G., Poline, J.B., Greicius, M., Thirion, B., 2016. Transport on Riemannian Manifold for Connectivity-Based Brain Decoding. IEEE Transactions on Medical Imaging 35, 208–216.

960

961

962

Ogawa, A., 2021. Time-varying measures of cerebral network centrality correlate with visual saliency during movie watching. Brain and Behavior 11, e2334. https://onlinelibrary.wiley.com/doi/pdf/10.1002/brb3.2334.

963

964

965
966

967

Ostroff, C., 1993. Comparing correlations based on individual-level and aggregated data. Journal of Applied Psychology 78, 569.

968

969

970

971

Pawela, C.P., Biswal, B.B., Cho, Y.R., Kao, D.S., Li, R., Jones, S.R., Schulte, M.L., Matloub, H.S., Hudetz, A.G., Hyde, J.S., 2008. Resting-state functional connectivity of the rat brain. Magnetic Resonance in Medicine 59, 1021-1029.

972

973

974

Petersen, A., Zhao, J., Carmichael, O., Müller, H.G., 2016. Quantifying Individual Brain Connectivity with Functional Principal Component Analysis for Networks. Brain Connectivity 6, 540-547.

975

976

977

Poline, J.B., Worsley, K.J., Evans, A.C., Friston, K.J., 1997. Combining spatial extent and peak intensity to test for activations in functional imaging. NeuroImage 5, 83-96.

978

979

980

Poline, J.B., Worsley, K.J., Holmes, A.P., Frackowiak, R.S.J., Friston, K.J., 1995. Estimating Smoothness in Statistical Parametric Maps: Variability of p Values. Journal of Computer Assisted Tomography 19, 788-796.

981

982

983

984

Richiardi, J., Achard, S., Bunke, H., Van De Ville, D., 2013. Machine Learning with Brain Graphs: Predictive Modeling Approaches for Functional Imaging in Systems Neuroscience. IEEE Signal Processing Magazine 30, 58-70.

985

986

Rosner, B., Donner, A., Hennekens, C., 1977. Estimation of interclass correlation from familial data. Applied statistics 26, 179-187.

987

Salvador, R., Suckling, J., Schwarzbauer, C., Bullmore, E., 2005. Undirected
graphs of frequency-dependent functional connectivity in whole brain networks. Philosophical Transactions of the Royal Society B: Biological Sciences 360, 937-946.

Stanley, M.L., Moussa, M.N., Paolini, B.M., Lyday, R.G., Burdette, J.H., Laurienti, P.J., 2013. Defining nodes in complex brain networks. Frontiers in Computational Neuroscience 7.

Student, 1914. The elimination of spurious correlation due to position in time and space. Biometrika 10, 179-181.

Termenon, M., Jaillard, A., Delon-Martin, C., Achard, S., 2016. Reliability of graph analysis of resting state fMRI using test-retest dataset from the human connectome project. Neuroimage 142, 172-187.

Tomasi, D., Volkow, N.D., 2010. Functional connectivity density mapping. Proceedings of the National Academy of Sciences 107, 9885-9890.

Tooley, U.A., Bassett, D.S., Mackey, A.P., 2021. Environmental influences on the pace of brain development. Nature Reviews Neuroscience 22, 372-384.

Triana, A.M., Glerean, E., Saramäki, J., Korhonen, O., 2020. Effects of spatial smoothing on group-level differences in functional brain networks. Network Neuroscience 4, 556-574.

Uddin, L.Q., Kelly, A.C., Biswal, B.B., Margulies, D.S., Shehzad, Z., Shaw, D., Ghaffari, M., Rotrosen, J., Adler, L.A., Castellanos, F.X., Milham, M.P., 2008. Network homogeneity reveals decreased integrity of default-mode network in ADHD. Journal of Neuroscience Methods 169, 249-254.
1010 Van Dijk, K.R.A., Hedden, T., Venkataraman, A., Evans, K.C., Lazar, S.W.,
1011 Buckner, R.L., 2010. Intrinsic Functional Connectivity As a Tool For
1012 Human Connectomics: Theory, Properties, and Optimization. Journal of
1013 Neurophysiology 103, 297-321.

1014 Whitlow, C.T., Casanova, R., Maldjian, J.A., 2011. Effect of Resting-State
1015 Functional MR Imaging Duration on Stability of Graph Theory Metrics of
1016 Brain Network Connectivity. Radiology 259, 516-524.

1017 Wilson, D.J., 2019. The harmonic mean *p*-value for combining dependent
1018 tests. Proceedings of the National Academy of Sciences 116, 1195-1200.

1019 Worsley, K.J., Evans, A.C., Marrett, S., Neelin, P., 1992. A Three-Dimensional Statistical Analysis for CBF Activation Studies in Human Brain. Journal of Cerebral Blood Flow & Metabolism 12, 900-918. Publisher: SAGE Publications Ltd STM.

1023 Worsley, K.J., Marrett, S., Neelin, P., Vandal, A.C., Fris-
ton, K.J., Evans, A.C., 1996. A unified statistical ap-
proach for determining significant signals in images of cere-
bral activation. Human Brain Mapping 4, 58-73. _eprint:
https://onlinelibrary.wiley.com/doi/pdf/10.1002/%28SICI%291097-
0193%281996%294%3A1%3C58%3A%3AAID-HBM4%3E3.0.CO%3B2-O.

1029 Yan, C.G., Cheung, B., Kelly, C., Colcombe, S., Craddock, R.C., Di Martino, A., Li, Q., Zuo, X.N., Castellanos, F.X., Milham, M.P., 2013. A comprehensive assessment of regional variation in the impact of head micromovements on functional connectomics. NeuroImage 76, 183-201.
Ye, J., Lazar, N., Li, Y., 2011. Sparse geostatistical analysis in clustering fMRI time series. Journal of Neuroscience Methods 199, 336-345.

Ye, J., Lazar, N.A., Li, Y., 2009. Geostatistical analysis in clustering fMRI time series. Statistics in Medicine 28, 2490-2508.

Zang, Y., Jiang, T., Lu, Y., He, Y., Tian, L., 2004. Regional homogeneity approach to fMRI data analysis. NeuroImage 22, 394-400.

Zhang, C., Cahill, N., Arbabshirani, M., White, T., Baum, S., Michael, A., 2016. Sex and age effects of functional connectivity in early adulthood. Brain Connectivity 6.

## 1042 Appendix A. Brain functional connectivity review

The literature review was conducted on PubMed using the keywords "brain connectivity graph resting state 'human connectome project'" on September 30, 2021. The search returned 32 papers written between 2014 and 2021. Out of those papers, 5 were not open access and 2 papers were literature reviews, and were not considered further. 3 papers were either using seed-based or voxel-to-voxel correlation. Out of the remaining 24 papers 71% (17/24) first averaged voxels before computing the inter-regional correlations and 88% (21/24) employed some kind of spatial aggregation method, including but not limited to averaging over voxels, ICA or dictionary learning.
## 1052 Appendix B. Hypotheses for the spatio-temporal model

The assumptions on the model can be written as follows. For any $i, i' \in \mathcal{C}$ and $s, t = 1, \dots, T$,

$$E[X_i(t)] = E[\varepsilon_i(t)] = E[e(t)] = 0,$$

$$E[X_i(s)X_i(t)] = E[\varepsilon_i(s)\varepsilon_i(t)] = E[e(s)e(t)] = 0,$$

$$E[X_i(s)\varepsilon_{i'}(t)] = E[X_i(s)e(t)] = E[\varepsilon_i(s)e(t)] = 0,$$

$$E[e(t)^2] = \sigma_e^2.$$

Let $\Sigma$ be the covariance matrix of the vector $(Y_i(t))_{i \in \mathcal{C}, t=1,\dots,T}$. In this paper, we assume without referring specifically to this assumption that the parameters $\sigma_j^2$, $\sigma_\varepsilon^2$, $\sigma_e^2$, $\rho_{ii'}$, $\eta_{ii'}$, $r_{jj'}$ are such that $\Sigma$ is a positive definite matrix.

We also assume that the random variables are independent in time. This is not overly restrictive: in particular, if the random variables have long memory, after a wavelet decomposition, the random variables can be approximated to be decorrelated in time for large wavelet scales (Moulines et al., 2007). In addition, assuming that the $X_i$'s are centered is coherent as it is a well-known fact that a wavelet decomposition based on a wavelet mother with $K$ vanishing moments cancels out every polynomial trend with degree $K-1$.

Finally, to apply the law of large numbers, we also assume that all random variables are absolutely integrable, that is $E[|Z_i(t)|] < \infty$ for $Z = X, \varepsilon, e$, $i \in \mathcal{C}$ and $t = 1, \dots, T$.
## 1068 Appendix C. Properties of the estimators of interest

For any set of indices $E$ with cardinality $\#E$, we let

$$ \bar{\rho}_E = \frac{1}{(\#E)^2} \sum_{i,i' \in E} \rho_{ii'} \quad \text{and} \quad \bar{\eta}_E = \frac{1}{(\#E)^2} \sum_{i,i' \in E} \eta_{ii'}. \qquad (C.1) $$

The results of the paper are based on this proposition:

1069 **Proposition Appendix C.1.** Consider the notation of Section 2.1 and
1070 assumptions described in Appendix B. Let $j, j' \in \{1, \dots, J\}$.

1071
(i) Let $E \subseteq \mathcal{R}_j$, then for any $t = 1, \dots, T$

$$ \mathrm{Var}[\bar{X}_E(t)] = \sigma_j^2 \bar{\rho}_E \qquad (C.2) $$

$$ \mathrm{Var}[\bar{\varepsilon}_E(t)] = \sigma_{\varepsilon}^2 \bar{\eta}_E = \mathcal{O}(1/(\#E)) \qquad (C.3) $$

$$ \mathrm{Var}[\bar{e}_E(t)] = \mathrm{Var}(e(t)) = \sigma_e^2 \qquad (C.4) $$

$$ \mathrm{Var}[\bar{Y}_E(t)] = \sigma_j^2 \bar{\rho}_E + \sigma_{\varepsilon}^2 \bar{\eta}_E + \sigma_e^2. \qquad (C.5) $$

1073 (ii) Let $E \subseteq \mathcal{R}_j$ and $E' \subseteq \mathcal{R}_{j'}$, then

$$ \mathrm{Cov}[\bar{Y}_E(t), \bar{Y}_{E'}(t)] = \begin{cases} \sigma_j \sigma_{j'} r_{jj'} + \sigma_e^2 & \text{if } j \neq j' \\ \sigma_j^2 \bar{\rho}_{E,E'} + \sigma_e^2 & \text{if } j = j' \end{cases} \qquad (C.6) $$

1074 where

$$ \bar{\rho}_{E,E'} = \frac{1}{(\#E)(\#E')} \sum_{i \in E, i' \in E'} \rho_{|i-i'|}. \qquad (C.7) $$

(iii) Let $i \in E \subseteq \mathcal{R}_j$ and $i' \in E' \subseteq \mathcal{R}_{j'}$ and assume for any $i \in E$ and $i' \in E'$
$|i - i'| \ge p$ (in the case $j = j'$). Then as $T \to \infty$, the following statements
hold almost surely.

$$ \hat{\sigma}^2(\mathbf{Y}_i) \xrightarrow{a.s.} \mathrm{Var}[Y_i(1)] \quad \text{and} \quad \widehat{\mathrm{Cov}}[\mathbf{Y}_i, \mathbf{Y}_{i'}] \xrightarrow{a.s.} \mathrm{Cov}[Y_i(1), Y_{i'}(1)] \quad (C.8) $$

$$ \hat{\sigma}^2(\bar{\mathbf{Y}}_E) \xrightarrow{a.s.} \mathrm{Var}[\bar{\mathbf{Y}}_E(1)] \quad \text{and} \quad \widehat{\mathrm{Cov}}[\bar{\mathbf{Y}}_E, \bar{\mathbf{Y}}_{E'}] \xrightarrow{a.s.} \mathrm{Cov}[\bar{\mathbf{Y}}_E(1), \bar{\mathbf{Y}}_{E'}(1)]. \quad (C.9) $$

Proposition Appendix C.1 is given without proof. (i)-(ii) ensue from the model (1) while (iii) is quite straightforward since we have assumed independence in time.

As seen from Proposition Appendix C.1, the quantity $\bar{\eta}_E$ is related to the correlation structure of the local noise. By assuming this noise to be p-dependent (that is $\eta_\delta = 0$ when $\delta \ge p$), it is clear that the larger $\#E$ the smaller $\bar{\eta}_E$.

## Appendix D. Consistency results for the existing estimators

### Appendix D.1. Consistency of $\hat{r}_{jj'}^{CA}$

Proposition Appendix C.1 shows $\hat{r}_{jj'}^{CA}$ is a strongly consistent estimator of $r_{jj'}^{CA}$ as $T \to \infty$ where

$$ r_{jj'}^{CA} = r_{jj'} \frac{1 + \sigma_{e,jj'}^2 / r_{jj'}}{\sqrt{(\bar{\rho}_{R_j} + \sigma_{\varepsilon,j}^2 \bar{\eta}_{R_j} + \sigma_{e,j}^2)(\bar{\rho}_{R_{j'}} + \sigma_{\varepsilon,j'}^2 \bar{\eta}_{R_{j'}} + \sigma_{e,j'}^2)}}. \quad (D.1) $$

Another way to correct the size effect is to compensate the inter-correlation by the intra-correlation. This would lead to the following estimator:

$$ \hat{r}_{jj'}^{\widetilde{\Lambda}^{\mathrm{C}}} = \frac{1}{N_j N_{j'}} \left( \sum_{i, i' \in \mathcal{R}_j} \widehat{\mathrm{Cor}}(\mathbf{Y}_i, \mathbf{Y}_{i'}) \sum_{i, i' \in \mathcal{R}_{j'}} \widehat{\mathrm{Cor}}(\mathbf{Y}_i, \mathbf{Y}_{i'}) \right)^{1/2} \hat{r}^{\Lambda^{\mathrm{C}}}. \quad (D.2) $$

The two estimators (5) and (D.2) have the important property to remove the size effect (since when $\sigma_\varepsilon = \sigma_e = 0$, $r_{jj'}^{AC} = r_{jj'}$). Note that both estimators tend to the same limit.
## Appendix D.2. Consistency of $\hat{r}_{jj'}^{AC}$

Proposition Appendix C.1 shows that $\hat{r}_{jj'}^{AC}$ is a strongly consistent estimator of $r_{jj'}^{AC}$ given by

As revealed by (D.1) and (D.3), $\hat{r}_{jj'}^{CA}$ and $\hat{r}_{jj'}^{AC}$ do not converge toward $r_{jj'}$ when a local noise or global noise is present. We could ask why $\hat{r}_{jj'}^{CA}$ is interesting. Actually, a first spatial averaging tends to decrease the effect of the local noise. Indeed, when $\sigma_e^2 = 0$ (and with equal unit variances to simplify), we have

Hence, if we expect that $\bar{\rho}_{R_j} \approx \bar{\rho}_{R_{j'}} \approx 1$, $\hat{r}_{jj'}^{CA}$ will be a better estimator since $\bar{\eta}_{R_j} = \mathcal{O}(1/N_j)$. A natural compromise between $\hat{r}_{jj'}^{CA}$ and $\hat{r}_{jj'}^{AC}$ can be defined using local neighborhood as defined by $\ell_{CA}$.

## Appendix D.3. Consistency of $\hat{r}_{jj'}^{R}$

From Proposition Appendix C.1, as $T \to \infty$

$$ \frac{1}{4} \sum_{\alpha, \beta=1}^{2} \widehat{\text{Cor}}(\mathbf{Y}_{i_{\alpha}^{(b)}}, \mathbf{Y}_{i_{\beta}^{(b)}}) \xrightarrow{\text{a.s.}} \frac{\sigma_j \sigma_{j'} r_{jj'} + \sigma_e^2}{\sqrt{(\sigma_j^2 + \sigma_e^2 + \sigma_e^2)(\sigma_{j'}^2 + \sigma_e^2 + \sigma_e^2)}} $$
and

$$ \widehat{\text{Cor}}(\mathbf{Y}_{i_1^{(b)}}, \mathbf{Y}_{i_2^{(b)}}) \xrightarrow{\text{a.s.}} \frac{\sigma_j^2 \rho_{\delta} + \sigma_e^2}{\sigma_j^2 + \sigma_e^2 + \sigma_e^2}, $$
whereby we deduce that $\hat{r}^R$ is a strongly consistent estimator of

$$ r_{jj'}^{R} = r_{jj'} \frac{1 + \sigma_{e,jj'}^{2}/r_{jj'}}{\sqrt{|(\rho_{\delta} + \sigma_{e,j}^{2})(\rho_{\delta} + \sigma_{e,j'}^{2})|}}. \quad (D.4) $$
From (D.4), we observe that when $\sigma_e = 0$ then for any unknown value of $\sigma_\epsilon$, $\hat{r}_{jj'}^R$ estimates consistently $r_{jj'}/|\rho_\delta|$ which should be close to $r_{jj'}$ if we take $\delta = p$ and expect that $\rho_p$ is close to 1. In other words, the estimator $\hat{r}_{jj'}^R$ is robust to the size of the regions and robust to a possible local noise.

To reduce the assumption that $\rho_p$ is close to 1, we can combine this idea of replicates with local averaging. This is the topic of the next section.

## 1112 Appendix E. Consistency of $\hat{r}_{jj'}^D$

The following result is the key ingredient:

**Proposition Appendix E.1.** Under the notation of this section, as $T \to \infty$, the following statements hold almost surely.

$$ (i) \qquad \widehat{\mathrm{Cov}}(\mathbf{Y}_{i^{(b)}}, \mathbf{Y}_{k^{(b)}}, \mathbf{Y}_{i'^{(b)}} - \mathbf{Y}_{k'^{(b)}}) \xrightarrow{a.s.} \sigma_j \sigma_{j'} r_{jj'}. $$

$$ (ii) \qquad 2\hat{s}^2(\mathbf{Y}_{i^{(b)}}, \mathbf{Y}_{k^{(b)}}, \mathbf{Y}_{k'^{(b)}}) \xrightarrow{a.s.} 2\sigma_j^2 + 2\sigma_\epsilon^2. \qquad (E.1) $$

1115 *Proof.* (i) Using the independence in time, it is clear that the left-hand side converges almost surely to $\mathrm{Cov}(Y_{i^{(b)}}(1) - Y_{k^{(b)}}(1), Y_{i'^{(b)}}(1) - Y_{k'^{(b)}}(1)) = \sigma_j \sigma_{j'} r_{jj'} + \sigma_e^2 - 2\sigma_e^2 + \sigma_e^2$ since the two regions $\mathcal{R}_k$ and $\mathcal{R}_{k'}$ are disconnected, which leads to the result.

1116
1117
1118

(ii) In the same way, the left-hand side tends to $\mathrm{Var}(Y_{i^{(b)}}(1) - Y_{k^{(b)}}(1)) + \mathrm{Var}(Y_{i^{(b)}}(1) - Y_{k'^{(b)}}(1)) - \mathrm{Var}(Y_{k^{(b)}}(1) - Y_{k'^{(b)}}(1)) = \sigma_j^2 + \sigma_k^2 + \sigma_{k'}^2 + 4\sigma_\epsilon^2 - \sigma_k^2 - \sigma_{k'}^2 - 2\sigma_\epsilon^2$ which yields the stated limit. $\square$

1122 In other words, Proposition Appendix E.1 shows that $\hat{r}_{jj'}^D$ is a strongly
1123 consistent estimator of $r_{jj'}^D$ given by

$$r^D = r_{jj'} \frac{1}{\sqrt{(1 + \sigma_{\varepsilon,j}^2)(1 + \sigma_{\varepsilon,j'}^2)}} \quad (E.2)$$

1124 which, in the situation where $\sigma_\varepsilon = 0$, is nothing else than $r_{jj'}$.

## 1125 Appendix F. Consistency of $\hat{r}_{jj'}^{RD}$

1126 The following result is a consequence of Propositions Appendix C.1- Appendix G.1.

1127 **Proposition Appendix F.1.** As $T \to \infty$, the following statements hold almost surely.

1128 (i) For any $i_1, i_2 \in \mathcal{R}_j, i'_1, i'_2 \in \mathcal{R}_j, i_k \in \mathcal{R}_k$ and $i_{k'} \in \mathcal{R}_{k'}$, such that $|i_2 - i_1| = |i'_2 - i'_1| = \delta \ge p$

$$\sqrt{|\widetilde{\operatorname{Cor}}(\mathbf{Y}_{i_1}, \mathbf{Y}_{i_2}; \mathbf{Y}_{i_k}, \mathbf{Y}_{i_{k'}})\widetilde{\operatorname{Cor}}(\mathbf{Y}_{i'_1}, \mathbf{Y}_{i'_2}; \mathbf{Y}_{i_k}, \mathbf{Y}_{i_{k'}})|} \xrightarrow{\text{a.s.}} \frac{\sigma_j \sigma_{j'} |\rho_\delta|}{\sqrt{(\sigma_j^2 + \sigma_\varepsilon^2)(\sigma_{j'}^2 + \sigma_\varepsilon^2)}} \quad (\text{F.1})$$

1129 (ii) For any $\nu$-neighborhoods $\mathcal{V}_{j_1}, \mathcal{V}_{j_2} \in \mathcal{R}_j, \mathcal{V}_{j'_1}, \mathcal{V}_{j'_2} \in \mathcal{R}_{j'}, \mathcal{V}_k \in \mathcal{R}_k, \mathcal{V}_{k'} \in \mathcal{R}_{k'}$, such that for any $i_1 \in \mathcal{V}_{j_1}, i_2 \in \mathcal{V}_{j_2}, i'_1 \in \mathcal{V}_{j'_1}, i'_2 \in \mathcal{V}_{j'_2}, |i_1 - i_2| = |i'_1 - i'_2| = \delta \ge p$

$$|\widetilde{\operatorname{Cor}}(\bar{\mathbf{Y}}_{\mathcal{V}_{j_1}}, \bar{\mathbf{Y}}_{\mathcal{V}_{j_2}}; \bar{\mathbf{Y}}_{\mathcal{V}_k}, \bar{\mathbf{Y}}_{\mathcal{V}_{k'}})\widetilde{\operatorname{Cor}}(\bar{\mathbf{Y}}_{\mathcal{V}_{j'_1}}, \bar{\mathbf{Y}}_{\mathcal{V}_{j'_2}}; \bar{\mathbf{Y}}_{\mathcal{V}_k}, \bar{\mathbf{Y}}_{\mathcal{V}_{k'}})| \xrightarrow{\text{a.s.}} \frac{\sigma_j^2 \sigma_{j'}^2 \rho_{\mathcal{V},\mathcal{V}'}^2 \delta}{(\sigma_j^2 \bar{\rho}_{\mathcal{V}} + \sigma_\varepsilon^2 \bar{\eta}_{\mathcal{V}})(\sigma_{j'}^2 \bar{\rho}_{\mathcal{V}} + \sigma_\varepsilon^2 \bar{\eta}_{\mathcal{V}})} \quad (\text{F.2})$$

1132 where $\mathcal{V}, \mathcal{V}'$ are two $\nu$-neighborhoods at distance $\delta$.

1133 Propositions Appendix G.1-Appendix F.1 show that $\hat{r}_{jj'}^{RD}$ is a strongly consistent estimator of $r_{jj'}^{RD}$ given by

$$r_{jj'}^{RD} = \frac{r_{jj'}}{|\rho_{\delta}|} \quad (\text{F.3})$$
## Appendix G. Consistency of localized versions of estimators

### Appendix G.1. Consistency of $\hat{r}_{jj'}^{\ell_{CA}}$

We can apply Proposition Appendix C.1 to show that $r_{jj'}^{\ell_{CA}}$ is a strongly consistent estimator of

$$r_{jj'}^{\ell_{CA}} = r_{jj'} \frac{1 + \sigma_{e,jj'}^2 / r_{jj'}}{\sqrt{(\bar{\rho}_{\nu} + \sigma_{\varepsilon,j}^2 \bar{\eta}_{\nu} + \sigma_{e,j}^2)(\bar{\rho}_{\nu} + \sigma_{\varepsilon,j'}^2 \bar{\eta}_{\nu} + \sigma_{e,j'}^2)}} \quad (G.1)$$

where $\mathcal{V}$ is any $\nu$-neighborhood. When there is no global noise ($\sigma_e = 0$) and for moderate $\nu$, it may be expected than the denominator of $r_{jj'}^{\ell_{CA}}$ is closer to 1 than the ones of $r_{jj'}^{\text{CA}}$ and $r_{jj'}^{\text{AC}}$.

### Appendix G.2. Consistency of $\hat{r}_{jj'}^{\ell_{R}}$

Proposition Appendix C.1 shows that $\hat{r}_{jj'}^{\ell_{R}}$ is a strongly consistent estimator of $r_{jj'}^{\ell_{R}}$ defined by

$$r_{jj'}^{\ell_{R}} = r_{jj'} \frac{1 + \sigma_{e,jj'}^2 / r_{jj'}}{\sqrt{|(\bar{\rho}_{\nu,\nu',\delta} + \sigma_{e,j}^2)(\bar{\rho}_{\nu,\nu',\delta} + \sigma_{e,j'}^2)|}} \quad (G.2)$$

where $\bar{\rho}_{\nu,\nu',\delta}$ is defined by (C.7) with $\mathcal{V}$ and $\mathcal{V}'$ two $\nu$-neighborhoods at distance $\delta$. Similarly to the estimator $\hat{r}^R$, when $\sigma_e = 0$, the previous expression reduces to $r_{jj'}^{\ell_{R}} = r_{jj'} / |\bar{\rho}_{\nu,\nu',\delta}|$ and again it is not unreasonable to think that $\bar{\rho}_{\nu,\nu',\delta}$ is close to 1.

### Appendix G.3. Consistency of $\hat{r}_{jj'}^{\ell_{D}}$

The following result is an adaptation of Proposition Appendix E.1 to local averages.
**Proposition Appendix G.1.** As $T \to \infty$, the following statements hold almost surely.

(i)

$$ \widehat{\mathrm{Cov}}(\bar{\mathbf{Y}}_{\mathcal{V}_j^{(b)}}, \bar{\mathbf{Y}}_{\mathcal{V}_k^{(b)}}, \bar{\mathbf{Y}}_{\mathcal{V}_{j'}^{(b)}} - \bar{\mathbf{Y}}_{\mathcal{V}_{k'}^{(b)}}) \xrightarrow{\text{a.s.}} \sigma_j \sigma_{j'} r_{jj'}. $$

1152

(ii)

$$ 2\hat{s}^2(\bar{\mathbf{Y}}_{\mathcal{V}_j^{(b)}}, \bar{\mathbf{Y}}_{\mathcal{V}_k^{(b)}}, \bar{\mathbf{Y}}_{\mathcal{V}_{k'}^{(b)}}) \xrightarrow{\text{a.s.}} 2\sigma_j^2 \bar{\rho}_\nu + 2\sigma_\varepsilon^2 \bar{\eta}_\nu. \quad (G.3) $$

1153

1154

1155

Using Proposition Appendix G.1 (for which proof follows along similar lines as Proposition Appendix E.1), we deduce that $\hat{r}_{jj'}^{\ell D}$ is a strongly consistent estimator of $r_{jj'}^{\ell D}$ given by

$$ r_{jj'}^{\ell D} = r_{jj'} \frac{1}{\sqrt{(\bar{\rho}_{\nu} + \sigma_{\varepsilon,j}^2 \bar{\eta}_{\nu})(\bar{\rho}_{\nu} + \sigma_{\varepsilon,j'}^2 \bar{\eta}_{\nu})}} \quad (G.4) $$

1156

where $\nu$ is any $\nu$-neighborhood.

1157

Appendix G.4. Consistency of $\hat{r}_{jj'}^{\ell RD}$

1158

1159

Propositions Appendix G.1-Appendix F.1 show that $\hat{r}_{jj'}^{\ell RD}$ is a strongly consistent estimator of $r_{jj'}^{\ell RD}$ given by

$$ r_{jj'}^{\ell RD} = \frac{r_{jj'}}{|\bar{\rho}_{\nu,\nu',\delta}|} \quad (G.5) $$

1160

1161

where $\nu$ and $\nu'$ are two $\nu$-neighborhoods at distance $\delta$. Similarly to the previous estimator, $\hat{r}_{jj'}^{\ell RD}$ is robust to an additive global and local noise.

1162

Appendix G.5. Robustness of CCC differences in terms of threshold

1163

1164

1165

In the main text, we presented the difference of CCC for the different estimators based on a single threshold corresponding to 20% of the edges of the graph. Figure G.9 displays the variability of CCC according to different
number of edges selected to construct the graph. This shows the robustness
of our findings where *l*CA has always higher CCC than CA for all possible
thresholds.

<figure><img src="image_10.png" /><figcaption>Figure G.9: CCC of betweenness according to different choices of threshold and for the different estimators.</figcaption></figure>