#  Observation of Electronic Structure Modification
in the Hidden Order Phase of CeCoSi

Shin-ichi Kimura1,2,3<EMAIL_ADDRESS>Hiroshi Watanabe1,2
Shingo Tatsukawa2 and Hiroshi Tanida4 1Graduate School of Frontier
Biosciences1Graduate School of Frontier Biosciences Osaka University Osaka
University Suita Suita Osaka 565-0871 Osaka 565-0871 Japan
2Department of Physics Japan
2Department of Physics Graduate School of Science Graduate School of Science
Osaka University Osaka University Toyonaka Toyonaka Osaka 560-0043 Osaka
560-0043 Japan
3Institute for Molecular Science Japan
3Institute for Molecular Science Okazaki Okazaki Aichi 444-8585 Aichi 444-8585
Japan
4Liberal Arts and Sciences Japan
4Liberal Arts and Sciences Toyama Prefectural University Toyama Prefectural
University Imizu Imizu Toyama 939-0398 Toyama 939-0398 Japan Japan

###### Abstract

CeCoSi with no local inversion symmetric crystal structure ($P4/nmm$) exhibits
a phase transition of unknown origin (Hidden Order: HO) at about 12 K
($T_{0}$) above the antiferromagnetic transition temperature ($T_{\rm N}$
$=9.4$ K). The electronic structure change across $T_{0}$ was investigated
with high-precision optical reflection spectroscopy. The optical spectrum
changed from a typical metallic behavior above $T_{0}$ to a gap-like structure
at around 15 meV below $T_{0}$. The gap-like structure was unchanged across
$T_{\rm N}$ except for the narrowing of the Drude component of carriers due to
the suppression of magnetic fluctuations. This result suggests a slight change
from the typical metallic electronic structure above $T_{0}$ to that with an
energy gap near the Fermi level in the HO phase. The change in electronic
structure in the HO phase was concluded to be due to electron/valence
instability.

Materials without inversion symmetry have recently attracted much attention
due to their topological phenomena [1]. The electric field generated by the
non-inversion symmetry produces Weyl semimetallic electronic structures with
strong spin-orbit interactions in the bulk [2, 3] and giant Rashba effects on
the surface of solids [4]. Such novel characteristic physical properties and
electronic structures can be realized in materials with non-inversion
symmetric crystal structures.

In rare-earth compounds with non-inversion symmetric / non-centrosymmetric
crystal structures, fascinating physical properties have been observed mainly
in superconducting properties: A Cooper pair with a helical spin structure
produces a superconducting state with a high upper critical field in global
non-centrosymmetric materials CePt3Si, CeCoGe3, and CeRhIn3. [5, 6, 7]. Also,
a locally non-centrosymmetric material, CeRh2As2, shows a superconducting
transition with changing the parity by applying magnetic field [8]. Locally
non-centrosymmetric Kondo semiconductors, Ce$M_{2}$Al10 ($M=$ Ru, Os), show an
antiferromagnetic transition with a tiny magnetic moment, namely hidden order
(HO), due to an electronic instability [9, 10, 11].

CeCoSi, the title compound, is one of the locally non-centrosymmetric rare-
earth compounds. It has a CeFeSi-type crystal structure ($P4/nmm$,
$D^{7}_{4h}$, No. 129), in which there is no local inversion symmetry on the
Ce site [12]. Since the electronic specific-heat coefficient $\gamma$ is
almost consistent with that of LaCoSi without $4f$ electrons, the $4f$
electrons are almost localized at low temperatures [13]. The Néel temperature
($T_{\rm N}$) is 9.4 K and the first excited state of the $4f$ crystal-field
splitting is about 100 K ($\sim 10$ meV) evaluated by the specific heat [13]
and inelastic neutron scattering (INS) measurements [14]. Recently, HO, in
which the order parameter has not been revealed yet, has been observed at
$12-13$ K ($T_{0}$) higher than $T_{\rm N}$ [13, 15]. The origin is still
unclear even though a theoretical work on the presence of odd-parity
multipoles, originating from staggered antiferromagnetic and antiferro-
quadrupole orderings, has been proposed [16]. Experimentally, $T_{0}$ shifts
to the high-temperature side, and the anomaly at $T_{0}$ becomes more
pronounced by applying pressure, i.e., 38 K at 1.5 GPa [17, 18] and also by
using magnetic field [15]. With these results, some anomalies in the HO phase
have been observed: The symmetry of the Co site becomes lowering observed by
NQR [19], the crystal symmetry becomes lowering slightly detected by a
synchrotron x-ray diffraction [20], and the sign of the magnetic resistance
changes at $T_{0}$ [15]. These results suggest that the electronic structure
changes at $T_{0}$, but the direct evidence of the electronic structure
modification in the HO phase has not been observed.

Optical measurements such as optical conductivity ($\sigma_{1}(\omega)$) and
angle-resolved photoelectron spectroscopy (ARPES) are essential in studying
HO. In URu2Si2 [21], the most famous heavy-fermion material with HO, a
pseudogap was reported in the infrared $\sigma_{1}(\omega)$ spectra in the
early stages of a series of studies [22]. In addition, in the HO phase of
Ce$M_{2}$Al10, a pseudogap was observed in the infrared $\sigma_{1}(\omega)$
spectra only in a certain direction [10, 11]. These examples show that
infrared spectroscopy is useful for observing slight changes in electronic
structures.

In this Letter, we performed high-precision infrared spectroscopy to
investigate whether a change in the electronic structure of CeCoSi occurs at
$T_{0}$. We observed the change in the reflectivity ($R(\omega)$) spectrum and
the $\sigma_{1}(\omega)$ spectrum at $T_{0}$ and also at $T_{\rm N}$. Below
the photon energy ($\hbar\omega$) of about 300 meV, in addition to two optical
phonon peaks at $\hbar\omega$ $\sim 20$ and 40 meV, a monotonic increase
toward lower energies due to a Drude absorption of carriers and interband
transitions of low-energy excitations appeared. High-precision temperature
dependence measurements revealed that the electronic structure near the Fermi
level ($E_{\rm F}$) slightly changes to a gap-like form below $T_{0}$. At
$T_{\rm N}$, the relaxation time of the Drude parameter was discontinuously
changed, implying the suppression of magnetic fluctuations below the
antiferromagnetic transition.

Single-crystalline CeCoSi samples were synthesized by a Ce/Co eutectic flux
method [13]. The $R(\omega)$ spectra have been measured using the as-grown
(001) plane, where the electric vector of light was set perpendicular to the
$c$-axis ($E$ $\perp c$). Near-normal-incident $R(\omega)$ spectra were
acquired in a photon-energy range of 5 meV – 1.5 eV at 6 – 300 K to obtain
information on the temperature-dependent $\sigma_{1}(\omega)$ spectra via the
Kramers-Kronig analysis (KKA) [23]. To get the highly accurate $R(\omega)$
spectra less than $\pm 0.3$ %, we adopted a feed-back positioning system
combined with an in-situ gold evaporation method [24]. For extrapolation of
the low and high energy sides in KKA, we used the conventional methods of the
Hagen-Rubens function for metallic $R(\omega)$ spectra below the lowest
accessible energy and measuring $R(\omega)$ spectra up to 30 eV only at 300 K
[25] connected with a free-electron approximation
$R(\omega)\propto\omega^{-4}$ above the highest accessible energy [26].
Additionally, in order to clarify small changes in the $R(\omega)$ spectrum at
metallic reflectivity higher than 90 % for the electronic structure
modification at $T_{0}$, we measured the relative $R(\omega)$ change by
temperature from 6 K (below $T_{\rm N}$) to 20 K (above $T_{0}$) with fixing
the sample position. With this method, we were able to capture changes in the
$R(\omega)$ spectrum with an accuracy better than 0.1 %.

Figure 1:  Temperature-dependent reflectivity ($R(\omega)$) spectra of the as-
grown $(001)$ surfaces of CeCoSi in the photon energy $\hbar\omega$ range of
$5-1000$ meV. Dispersive peaks at the photon energies $\hbar\omega$ of about
20 and 42 meV originate from TO phonons. (Inset) Temperature-dependent
$R(\omega)$ spectra of CeCoSi in the $\hbar\omega$ range of $10-60$ meV at
lower temperatures than 80 K. The spectra at 6 and 10 K show a steep rise
structure below 15 meV.

Obtained $R(\omega)$ spectra of CeCoSi in a wide temperature range are shown
in Fig. 1. Although the overall structure is metallic, dispersive peaks due to
optical phonons were observed at about 20 and 42 meV. As the temperature is
lowered, $R(\omega)$ at $\sim 10$ meV increases by about 2 %, while it
decreases by about 2 % at $\sim 300$ meV, indicating that the relaxation time
increases at low temperatures, behaving like a typical metal. In ordinary
heavy-fermion systems, a double peak with a splitting energy of 0.25 eV
corresponding to the Ce $4f$ spin-orbit splitting width appears in this
region, and a significant temperature dependence associated with the
development of the Kondo effect is observed [27, 28]. However, in CeCoSi,
there is no such change, suggesting that the $4f$ electrons are very
localized, which is consistent with the speculation from the low electronic
specific-heat coefficient [13]. The inset shows the temperature variation of
the far-infrared $R(\omega)$ spectrum below 80 K. At 80 and 20 K, except for
two phonon peaks, the backgrounds of the $R(\omega)$ spectra monotonically
increase toward lower energies. However, the $R(\omega)$ spectra at 10 and 6 K
show a rapid increase with decreasing $\hbar\omega$ below 15 meV. This change
suggests that a variation in the electronic structure appears below $\sim 15$
meV.

Figure 2:  (a) Temperature-dependent $R(\omega)$ spectra divided by that at 20
K ($R(\omega,T)/R(\omega,{\rm 20~{}K})$) of CeCoSi. (b) Temperature dependence
of the spectral integrations [$\Xi(\Omega,T)=\int_{\Omega}R(\omega,T)d\omega$]
of the range $\Omega$ of 9–13 meV (open square) and 20–40 meV (solid circle)
in (a), which are representative regions of the rapid increase below 15 meV
and the gradual decrease above 15 meV, respectively, with decreasing
temperature. $\Xi(\Omega,T)$ is normalized by that value at $T=20$ K. Solid
lines are guides for the eye. $T_{0}$ and $T_{\rm N}$ are shown by vertical
thick-solid and dashed lines, respectively.

To examine the temperature variation of the $R(\omega)$ spectrum below 20 K in
detail, we measured the relative deviation to the $R(\omega)$ spectrum at 20 K
as shown in Fig. 2(a). The data show that the relative spectral change is
about 0.1 % and less, which is hardly detected using conventional $R(\omega)$
measurements. As observed in the inset of Fig. 1, the intensity above 15 meV
decreases and conversely increases at lower energies with decreasing
temperature. To better clarify this change, we plotted the integration
[$\Xi(\Omega,T)=\int_{\Omega}R(\omega,T)d\omega$] normalized by the value at
20 K over the ranges $\Omega$ of 9–13 meV and 20–40 meV as a function of
temperature in Fig. 2(b).

$\Xi(\Omega,T)$ of $\Omega=9-13$ meV shows only a change of $\pm 0.1$ % or
less above $T_{\rm N}$, and is almost constant within the measurement accuracy
of about $\pm 0.1$ %. On the other hand, the value discontinuously changes at
$T_{\rm N}$ and rapidly increases with decreasing temperature. $\Xi(\Omega,T)$
of $\Omega=20-40$ meV also has a step-like change at $T_{\rm N}$. The spectral
variation of $\Xi(\Omega,T)$ below $T_{\rm N}$ corresponds to an increase in
the relaxation time of the Drude component and is consistent with the fact
that the electrical resistivity drops below $T_{\rm N}$ [13] due to a
suppression of magnetic fluctuations by the magnetic ordering.

At around $T_{0}$, $\Xi(\Omega,T)$ of $\Omega=20-40$ meV shows the increasing
ratio with increasing temperature [$d\Xi(\Omega,T)/dT/\Xi(\Omega,{\rm
20~{}K})$] changes from $\sim 1\times 10^{-3}$ K-1 below $T_{0}$ to $\sim
5\times 10^{-4}$ K-1 above $T_{0}$. This result suggests that a modification
of the electronic structure near $E_{\rm F}$ occurs at $T_{0}$. To clarify the
electronic structure change in more detail, $\sigma_{1}(\omega)$ spectra,
which correspond to the joint density of states, were derived from KKA of
$R(\omega)$ spectra.

Figure 3:  Optical conductivity ($\sigma_{1}(\omega)$) spectra of CeCoSi at
representative temperatures of 20 K ($>$ $T_{0}$), 10 K ($T_{\rm N}$ $<T<$
$T_{0}$), and 6 K ($<$ $T_{\rm N}$) (solid lines). Dot-dashed lines are Drude
curves obtained from the values of electrical resistivity and the
$\sigma_{1}(\omega)$ values at $\hbar\omega\sim 5$ meV. The interband
components evaluated by subtracting the Drude components from the
$\sigma_{1}(\omega)$ spectra are shown by dashed lines.

Figure 3 shows $\sigma_{1}(\omega)$ spectra at representative temperatures
below $T_{\rm N}$, between $T_{\rm N}$ and $T_{0}$, and above $T_{0}$ to
investigate the change of electronic states at the different phases. Here, we
obtained the $R(\omega)$ spectra below 20 K by multiplying the relative
$R(\omega)$ spectra [$R(\omega,T)/R(\omega,{\rm 20~{}K})$] in Fig. 2(a) to the
$R(\omega)$ spectrum at 20 K shown in Fig. 1, and derived $\sigma_{1}(\omega)$
spectra from KKA of the $R(\omega)$ spectra. In this figure, the sharp peaks
at 20 and 42 meV originate from optical phonons, and no evident change was
observed in this temperature range. This result suggests that there is no
significant change in crystal structure or symmetry at $T_{0}$ and $T_{\rm
N}$, which is consistent with only a slight structural change observed with a
synchrotron x-ray diffraction [20]. Besides these sharp phonon peaks, the
spectra have several broad peaks, possibly due to interband transitions.
However, they are not treated in this study because they hardly change in
temperature. At 20 K, the spectral intensity monotonically increases toward
lower energies, which qualitatively represents the Drude component leading to
direct current conductivity. However, with decreasing temperature, it can be
seen that the intensity around 15 meV decreases more rapidly than that above
20 meV. This temperature dependence suggests a modulation in the electronic
structure at low temperatures from a simple conduction band at 20 K to that
with a gap-like structure. The upward structure toward low energy below 15 meV
becomes sharper at low temperatures, which reduces the intensity at around 15
meV. Since the upward structure can be regarded as the Drude component
originating from the carriers, the Drude components expected from the
electrical resistivity and the slope at the lowest accessible photon energy
(dot-dashed lines) are subtracted from the $\sigma_{1}(\omega)$ spectra as
shown by dashed lines [29]. It should be noted that the Drude components in
Fig. 3 were evaluated as no interband transition component at the lowest
energy in $\sigma_{1}(\omega)$ spectra. This assumption may be unrealistic
because interband transitions can exist even at low energy. However, since the
same evaluation method for the Drude components is used for all
$\sigma_{1}(\omega)$ spectra, the temperature dependence of the
$\sigma_{1}(\omega)$ spectrum after the Drude component subtraction can be
compared with each other.

Compared to the Drude-subtracted $\sigma_{1}(\omega)$ spectrum at 20 K ($>$
$T_{0}$), the intensity below 20 meV is suppressed at lower temperatures than
$T_{0}$, indicating the appearance of an energy-gap-like structure. The Drude-
subtracted spectrum was almost unchanged across $T_{\rm N}$ from 10 K to 6 K.
This result suggests that the electronic structure near $E_{\rm F}$ is
modified at $T_{0}$, where an energy-gap-like form appears, but does not
change at $T_{\rm N}$.

Figure 4:  (a) Temperature-dependent $\sigma_{1}(\omega)$ spectra normalized
by that at 20 K of CeCoSi. The marks are the same as those in (b). (b)
Temperature dependence of the spectral integrations
[$\Sigma(\Omega,T)=\int_{\Omega}\sigma_{1}(\omega,T)d\omega$] of the range
$\Omega$ of 10–15 meV (solid circle), 15–20 meV (open circle), 20–25 meV (open
square), and 20–40 meV (solid square) in (a), which are representative regions
of relative $\sigma_{1}(\omega)$ spectra in (a). The $\Sigma(\Omega,T)$
intensity at each $\Omega$ is normalized by the value at $T=20$ K. Solid lines
are guides for the eye. $T_{0}$ and $T_{\rm N}$ are shown by vertical solid
and dashed lines, respectively.

To clarify the temperature dependence of $\sigma_{1}(\omega)$ in more detail,
the $\sigma_{1}(\omega)$ spectra below 20 K are divided by that at 20 K and
plotted as functions of temperature and photon energy in Fig. 4(a). At low
temperatures, the gap formation below 15 meV is clearly visible. To clarify
the spectral change at the ordering temperatures, the temperature dependence
of the integration of the $\sigma_{1}(\omega)$ spectrum
[$\Sigma(\Omega,T)=\int_{\Omega}\sigma_{1}(\omega,T)d\omega$] in the four
energy regions is shown in Fig. 4(b). In $\Sigma(\Omega,T)$ of $\Omega=10-15$
meV, a clear jump appears at $T_{\rm N}$, which is the same as the jump in
$R(\omega)$ at $T_{\rm N}$ in Fig. 2(b). This evidence can be attributed to
the increase in the relaxation time of carriers caused by the suppression of
magnetic fluctuations below $T_{\rm N}$. In addition, the increasing ratio of
$\Sigma(\Omega,T)$ ($d\Sigma(\Omega,T)/dT$) below $T_{0}$ is steeper than that
above $T_{0}$. This fact suggests the development of the gap-like structure
below $T_{0}$ shown in Fig. 3. The intensity of $\Sigma(\Omega,T)$ at higher
energies than 15 meV hardly decreases with decreasing temperature, confirming
that the size of the energy gap is about 15 meV. This gap-like structure does
not change at $T_{\rm N}$, which corresponds to the fact that
$d\Sigma(\Omega,T)/dT$ is almost constant before and after $T_{\rm N}$, even
though a step-like change appears at $T_{\rm N}$. In other words, the gap-like
structure observed at low temperatures in Fig. 3 starts to develop at $T_{0}$.
This fact suggests that the HO of CeCoSi is accompanied by a modification in
the electronic state near $E_{\rm F}$.

Here, we discuss the origin of the electronic structure modification below
$T_{0}$. An energy-gap-like structure with a gap size ($\Delta E$) of $\sim
15$ meV slightly appears below $T_{0}$. The amount of the $\sigma_{1}(\omega)$
intensity change at $\sim 10$ meV below $\Delta E$ is less than 10 % at 10 K
($T_{\rm N}$ $<T<$ $T_{0}$) compared to that at 14 K ($>$ $T_{0}$), and about
20 % in maximum at 6 K ($<$ $T_{\rm N}$) as shown in Fig. 4. The change is
much smaller than that of other HO materials; for instance, more than 50 % in
URu2Si2 [30] and CeOs2Al10 [10]. The reason is considered to be due to the
remaining large background of the Drude component and low-energy interband
transition that does not change with temperature. This energy gap-like
structure formation suggests that some bands close to $E_{\rm F}$, including
conduction bands, are modulated at $T_{0}$. The size of $\Delta E$ is similar
to a broad peak observed by INS [14] and the $4f$ crystal-field splitting
energy [13]. The optical gap observed here below $T_{0}$ may be related to
both the INS peak and the crystal field, but we must consider that optical
measurements detect the excitation at the momentum transfer of $q$ $\sim 0$,
in contrast to a finite $q$ detected by INS. The origin of the optical gap and
INS peak may be the $c$-$f$ hybridization band because the evidence of a weak
$c$-$f$ hybridization appears in the electrical resistivity data [13, 31].
Then, the valence of Ce-ions is slightly shifted from trivalent ($3+\delta$)
above $T_{\rm N}$. On the other hand, below $T_{\rm N}$, the Ce valence
becomes trivalent because the $4f$ electrons in Ce are fully localized [13].
Therefore, at temperatures slightly above $T_{\rm N}$, around $T_{0}$, an
electronic/valence instability can appear. By applying pressure, the anomaly
in the electrical resistivity at $T_{0}$ becomes more visible and moves to the
high-temperature side, suggesting that the electronic/valence instability is
grown with increasing pressure because the $c$-$f$ hybridization intensity
increases with applying pressure. The pressure dependence typically appears in
Ce-based heavy fermions [32]. The INS peak also shifts to the high-energy side
with increasing pressure, which is the same as the behavior of the anomaly at
$T_{0}$, suggesting the enlargement of the $c$-$f$ hybridization gap size.
Another possibility is that the electronic/valence instability may originate
from the characteristic degeneracy specific to the nonsymmorphic crystal
structure discussed in LaMnSi [33]. Optical measurements under pressure may
provide conclusive information for the origin of HO and its relation to the
INS peak.

The observed electronic/valence instability at $T_{0}$ is similar to that of
Ce$M_{2}$Al10. The crystal structure of Ce$M_{2}$Al10 is also a locally non-
centrosymmetric one ($Cmcm$), and the pressure-temperature phase diagram seems
similar to that of CeCoSi [34]. However, in Ce$M_{2}$Al10, the lattice
distortion at $T_{0}$ has yet to be observed, inconsistent with CeCoSi. Like
Ce$M_{2}$Al10, in URu2Si2, the lattice distortion in the HO phase has not been
observed [35, 36]. Since all of these phase transitions are considered to
originate from electronic/valence instability, optical measurements, including
ARPES, with high accuracy and resolution, will be essential in investigating
the origin of HO.

To summarize, we investigated the origin of the HO phase of CeCoSi using
temperature-dependent high-precision relative infrared $R(\omega)$
spectroscopy. The $R(\omega)$ spectra and the $\sigma_{1}(\omega)$ spectra
after the Kramer-Kronig analysis clearly show non-monotonic temperature
dependence at $T_{0}$ and $T_{\rm N}$. The anomaly at $T_{\rm N}$ originates
from the elongation of the carrier relaxation time due to the suppression of
magnetic fluctuations in the antiferromagnetic phase. On the other hand, the
change at $T_{0}$ is derived from the deviation of the low-energy interband
transition component, with the emergence of an energy-gap-like structure with
the size of $\sim 15$ meV. This result indicates that the electronic structure
modification due to electronic/valence instability appears in the HO phase of
CeCoSi.

We would like to thank Mr. Kiichi Ikai for his support of the IR experiments,
and UVSOR Synchrotron staff members for their support during synchrotron
radiation experiments. Part of this work was performed under the Use-of-UVSOR
Synchrotron Facility Program (Proposals No. 22IMS6015) of the Institute for
Molecular Science, National Institutes of Natural Sciences. This work was
partly supported by JSPS KAKENHI (Grant Nos. 20H04453 and 18H03683).

## References

  * [1] R. Moessner and J. E. Moore, Topological Phases of Matter (Cambridge University Press, 2021).
  * [2] B. A. Bernevig, Nat. Phys. 11, 698 (2015).
  * [3] N. P. Armitage, E. J. Mele, and A. Vishwanath, Rev. Mod. Phys. 90, 015001 (2018).
  * [4] A. Manchon, H. C. Koo, J. Nitta, S. M. Frolov, and R. A. Duine, Nat. Mater. 14, 871 (2015).
  * [5] E. Bauer, G. Hilscher, H. Michor, C. Paul, E. W. Scheidt, A. Gribanov, Y. Seropegin, H. Noël, M. Sigrist, and P. Rogl, Phys. Rev. Lett. 92, 027003 (2004).
  * [6] M.-a. Méasson, H. Muranaka, T. Kawai, Y. Ota, K. Sugiyama, M. Hagiwara, K. Kindo, T. Takeuchi, K. Shimizu, F. Honda, R. Settai, and Y. Ônuki, J. Phys. Soc. Jpn. 78, 124713 (2009).
  * [7] N. Kimura, K. Ito, H. Aoki, S. Uji, and T. Terashima, Phys. Rev. Lett. 98, 197001 (2007).
  * [8] S. Khim, J. F. Landaeta, J. Banda, N. Bannor, M. Brando, P. M. R. Brydon, D. Hafner, R. Küchler, R. Cardoso-Gil, U. Stockert, A. P. Mackenzie, D. F. Agterberg, C. Geibel, and E. Hassinger, Science 373, 1012 (2021).
  * [9] T. Nishioka, Y. Kawamura, T. Takesaka, R. Kobayashi, H. Kato, M. Matsumura, K. Kodama, K. Matsubayashi, and Y. Uwatoko, J. Phys. Soc. Jpn. 78, 123705 (2009).
  * [10] S. Kimura, T. Iizuka, H. Miyazaki, A. Irizawa, Y. Muro, and T. Takabatake, Phys. Rev. Lett. 106, 056404 (2011).
  * [11] S. Kimura, T. Iizuka, H. Miyazaki, T. Hajiri, M. Matsunami, T. Mori, A. Irizawa, Y. Muro, J. Kajino, and T. Takabatake, Phys. Rev. B 84, 165125 (2011).
  * [12] O. I. Bodak, E. I. Gladyshevskii, and P. I. Kripyakevich, J. Struct. Chem. 11, 283 (1970).
  * [13] H. Tanida, K. Mitsumoto, Y. Muro, T. Fukuhara, Y. Kawamura, A. Kondo, K. Kindo, Y. Matsumoto, T. Namiki, T. Kuwai, and T. Matsumura, J. Phys. Soc. Jpn. 88, 054716 (2019).
  * [14] S. E. Nikitin, D. G. Franco, J. Kwon, R. Bewley, A. Podlesnyak, A. Hoser, M. M. Koza, C. Geibel, and O. Stockert, Phys. Rev. B 101, 214426 (2020).
  * [15] H. Hidaka, S. Yanagiya, E. Hayasaka, Y. Kaneko, T. Yanagisawa, H. Tanida, and H. Amitsuka, J. Phys. Soc. Jpn. 91, 094701 (2022).
  * [16] M. Yatsushiro and S. Hayami, Phys. Rev. B 102, 195147 (2020).
  * [17] E. Lengyel, M. Nicklas, N. Caroca-Canales, and C. Geibel, Phys. Rev. B 88 155137 (2013).
  * [18] H. Tanida, Y. Muro, and T. Matsumura, J. Phys. Soc. Jpn. 87, 023705 (2018).
  * [19] M. Manago, H. Kotegawa, H. Tou, H. Harima, and H. Tanida, J. Phys. Soc. Jpn. 90, 023702 (2021).
  * [20] T. Matsumura, S. Kishida, M. Tsukagoshi, Y. Kawamura, H. Nakao, and H. Tanida, J. Phys. Soc. Jpn. 91, 064704 (2022).
  * [21] J. A. Mydosh and P. M. Oppeneer, Rev. Mod. Phys. 83, 1301 (2011).
  * [22] D. A. Bonn, J. D. Garrett, and T. Timusk, Phys. Rev. Lett. 61, 1305 (1988).
  * [23] S. Kimura and H. Okamura, J. Phys. Soc. Jpn. 82, 021004 (2013).
  * [24] S. Kimura, Jasco Rep. 50, 6 (2008).
  * [25] K. Fukui, R.-i. Ikematsu, Y. Imoto, M. Kitaura, K. Nakagawa, T. Ejima, E. Nakamura, M. Sakai, M. Hasumoto, and S. Kimura, J. Synchrotron Radiat. 21, 452 (2014).
  * [26] M. Dressel and G. Grüner, Electrodynamics of Solids (Cambridge University Press, 2002).
  * [27] S. Kimura, Y. S. Kwon, Y. Matsumoto, H. Aoki, and O. Sakai, J. Phys. Soc. Jpn. 85, 083702 (2016).
  * [28] S. Kimura, J. Sichelschmidt, and S. Khim, Phys. Rev. B 104, 245116 (2021).
  * [29] The expected Drude curves do not follow the experimental $\sigma_{1}(\omega)$ spectra because the $\sigma_{1}(\omega=0)$ values evaluated from the electrical resistivity ($\sim 2\times 10^{4}~{}\Omega^{-1}{\rm cm}^{-1}$ at 20 K, $\sim 2.5\times 10^{4}~{}\Omega^{-1}{\rm cm}^{-1}$ at 10 K, $\sim 1\times 10^{5}~{}\Omega^{-1}{\rm cm}^{-1}$ at 6 K) are much higher than the $\sigma_{1}(\omega)$ intensity at the lowest accessible energy. The same is usually observed in heavy fermion compounds at low temperatures due to the strong renormalization effect. In the Drude parameters [26], the dumping constant (corresponding to the inverse relaxation time) at 6 K ($<$$T_{\rm N}$) was evaluated as 2 meV, which is reduced from 5 meV for 20 and 10 K ($>$$T_{\rm N}$). The fact also suggests the suppression of magnetic fluctuations below $T_{\rm N}$.
  * [30] J. Levallois, F. Lévy-Bertrand, M. K. Tran, D. Stricker, J. A. Mydosh, Y.-K. Huang, and D. van der Marel, Phys. Rev. B 84, 184420 (2011).
  * [31] Y. Kawamura, K. Ikeda, A. N. B. A. Dalan, J. Hayashi, K. Takeda, C. Sekine, T. Matsumura, J. Gouchi, Y. Uwatoko, T. Tomita, H. Takahashi, and H. Tanida, J. Phys. Soc. Jpn. 91, 064714 (2022).
  * [32] H. Yamaoka, Y. Yamamoto, E. F. Schwier, F. Honda, Y. Zekko, Y. Ohta, J.-F. Lin, M. Nakatake, H. Iwasawa, M. Arita, K. Shimada, N. Hiraoka, H. Ishii, K.-D. Tsuei, and J. Mizuki, Phys. Rev. B 92, 235110 (2015).
  * [33] H. Tanida, H. Matsuoka, K. Mitsumoto, Y. Muro, T. Fukuhara, and H. Harima, J. Phys. Soc. Jpn. 91, 013704 (2022).
  * [34] K. Umeo, T. Ohsuka, Y. Muro, J. Kajino, and T. Takabatake, J. Phys. Soc. Jpn. 80, 064709 (2011).
  * [35] A. de Visser, F. E. Kayzel, A. A. Menovsky, J. J. M. Franse, J. van den Berg, and G. J. Nieuwenhuys, Phys. Rev. B 34, 8168 (1986).
  * [36] P. G. Niklowitz, C. Pfleiderer, T. Keller, M. Vojta, Y.-K. Huang, and J. A. Mydosh, Phys. Rev. Lett. 104, 106406 (2010).


# Robust correlation for aggregated data with spatial characteristics

Sophie Acharda, Jean-François Coeurjollya,
Pierre Lafaye de Micheauxb and Jonas Richiardic
a Univ. Grenoble Alpes, CNRS, Inria, Grenoble-INP, LJK, 38000 Grenoble,
France;
b School of Mathematics and Statistics, UNSW Sydney, NSW, Australia;
c Department of Radiology, Lausanne University Hospital and University of
Lausanne,
Lausanne, Switzerland

###### Abstract

The objective of this paper is to study the robustness of computation of
correlations under spatial constraints. The motivation of our paper is the
specific case of functional magnetic resonance (fMRI) brain imaging data,
where voxels are aggregated to compute correlations. In this paper we show
that the way the spatial components are aggregating to compute correlation may
have a strong influence on the resulting estimations. We then propose various
estimators which take into account this spatial structure.

keywords correlation; aggregated data; familial correlations; serial
correlations

##  1 Introduction

The use of aggregated data is particularly common in various field of
sciences. In social sciences, both individual and organizational data are
collected (Ostroff 1993). The objective is to evaluate the relationships among
variables at different level of analysis. For example, in (Ostroff 1993), the
author took the example of studying correlations between satisfaction and
technology at two levels, individual scores, denoted individual correlation,
and when individuals are grouped into organizations, also called
organizational correlation. Depending on three main factors mainly the
measurement errors, the variance within group and the variance-within ratios,
the ratio of organizational correlation and individual correlation was shown
to take values between -1 and 2. This means that it is crucial to identify the
way the data are generated.

In the studies of familial data (Rosner et al. 1977), specific characteristics
are obtained for different families with different sizes. The quantification
of family resemblance is studied in this context, with the estimation of
sibling correlations and/or parent-offspring correlations, see for a review
(Donner and Eliasziw 1991). Each data is composed with characteristics extract
for families with a mother and/or father and the children. The difficulties
arise because of the dependence between the children and the different number
of children per family. It was noted that correlation between the average of
time series and the average of correlations between all possible pairs of time
series are not equal in the majority of cases.

In psychology, (Vul et al. 2009) shows that the strength of relationship
depends on the reliability of the measurements. Finally, the problem of
correlation computation is also a common problem in ecology, geography,
climate studies …. The data collected in this field are attached to a spatial
position and usually with spatial correlation. This 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 ecological data (Liebhold and Sharov
1998), in fMRI data (Ye et al. 2011). The main difficulty in these analyses is
to take into account the spatial correlations in order to construct estimator
of correlations and testing procedures when the averaged variables are not
independent.

We are facing analogous difficulties in brain imaging, in particular when
studying so-called functional connectivity. Functional connectivity of the
brain is estimated from observations using non invasive techniques such as
electroencephalography (EEG), magnetoencephalography (MEG) or 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 between
the time series. It has already been shown that computation of connectomes is
affected by three main parameters: the length of the acquisition, the number
of regions and the chosen frequency band. In addition, the number of subjects
available in the sample will also play a role in terms of group comparisons
(Termenon et al. 2016). Using fMRI data, we record thousand of voxels. Each
region of the brain is then associated to a given set of voxels. The idea is
then to extract a representative of the set of voxels to attach one time
series to each region. The most common approach is to take the average of the
voxels. In this paper, we question this choice by studying statistically the
computation of correlation of this average of voxels.

The paper is organized as follows. In a first part, we are motivating the
definition of a new spatial model of fMRI. Then, computations of correlations
are described. Based on simulations, we illustrate the good behavior of the
newly introduced estimators.

##  2 Definition of the proposed spatial model for fMRI data

Let $\mathcal{C}$ denote a finite compact subset of indices of
$\mathbb{Z}^{d}$. This subset $\mathcal{C}$ will play the role of the set of
indices of voxels of the brain when $d=3$. The $J$ regions of interest of the
brain are represented through their set of voxels denoted by
${\mathcal{R}}_{j}$ for $j=1,\ldots,N$. The number of voxels in each region is
denoted by $\\#{\mathcal{R}}_{j}=N_{j}$. So,

$\mathcal{C}=\cup_{j=1}^{J}{\mathcal{R}}_{j}\mbox{ and
}\\#\mathcal{C}=\sum_{j=1}^{J}N_{j}.$

For any $i\in\mathcal{C}$, we assume observing the signal $Y_{i}(\cdot)$
sampled at times $t=1,\ldots,T$ which can be decomposed as follows

$Y_{i}(t)=X_{i}(t)+\varepsilon_{i}(t)+e(t),$ (1)

where $X_{i}(\cdot)$ represents the signal of interest,
$\varepsilon_{i}(\cdot)$ represents a noise contaminating locally the voxel
$i$ and the signal $e(\cdot)$ is a noise corrupting in the same way all the
voxels $i\in\mathcal{C}$. We now make a few assumptions on these different
components. First, we assume that all the random variables are centered, the
signals $X_{i}(\cdot)$, $\varepsilon_{i}(\cdot)$ and $e(\cdot)$ are mutually
independent and independent in time and that the global noise is
homoscedastic. This implies that for any $i,i^{\prime}\in\mathcal{C}$ and
$s,t=1,\ldots,T$ ($s\neq t$)

$\displaystyle E[X_{i}(t)]$ $\displaystyle=E[\varepsilon_{i}(t)]=E[e(t)]=0$
$\displaystyle E[X_{i}(s)X_{i}(t)]$
$\displaystyle=E[\varepsilon_{i}(s)\varepsilon_{i}(t)]=E[e(s)e(t)]=0$
$\displaystyle E[X_{i}(s)\varepsilon_{i^{\prime}}(t)]$
$\displaystyle=E[X_{i}(s)e(t)]=E[\varepsilon_{i}(s)e(t)]=0$ $\displaystyle
E[e(t)^{2}]$ $\displaystyle=\sigma^{2}_{e},$

Let us now describe the spatial nature of the signal and the local noise. For
any $i,i^{\prime}\in\mathcal{C}$ $j,j^{\prime}=1,\ldots,J$ ($j\neq
j^{\prime}$) and for all $t=1,\ldots,T$, we assume that there exists
$\sigma_{j}>0$, $\sigma_{e}\geq 0$, $r_{jj^{\prime}}\in[-1,1]$,
$\rho_{ii^{\prime}}\in(0,1]$, $\eta_{ii^{\prime}}\in[-1,1]$ such that

$E[X_{i}(t)X_{i^{\prime}}(t)]=\left\\{\begin{array}[]{ll}\sigma_{j}\sigma_{j^{\prime}}r_{jj^{\prime}}&\mbox{
if }i\in{\mathcal{R}}_{j},i^{\prime}\in{\mathcal{R}}_{j^{\prime}},j\neq
j^{\prime}\\\ \sigma_{j}^{2}\rho_{ii^{\prime}}&\mbox{ if
}i,i^{\prime}\in{\mathcal{R}}_{j}\end{array}\right.$

and

$E[\varepsilon_{i}(t)\varepsilon_{i^{\prime}}(t)]=\sigma^{2}_{\varepsilon}\eta_{ii^{\prime}}.$

The parameter $r_{jj^{\prime}}$ represents the correlation between two signals
of two different regions ${\mathcal{R}}_{j}$ and ${\mathcal{R}}_{j^{\prime}}$
and is called inter-correlation between regions ${\mathcal{R}}_{j}$ and
${\mathcal{R}}_{j^{\prime}}$ in the following. The parameter
$\rho_{ii^{\prime}}$ (resp. $\eta_{ii^{\prime}}$) represents the intra-
correlation between two signals (resp. 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 and that for each time $t$
and for $j=1,\ldots,J$, $(X_{i}(t),i\in{\mathcal{R}}_{j})$ (resp.
($\varepsilon_{i}(t),i\in\mathcal{C}$)) is a stationary random field observed
in ${\mathcal{R}}_{j}$ (resp. $\mathcal{C}$). We furthermore assume that both
the correlations $\rho_{ii^{\prime}}$ (for any
$i,i^{\prime}\in{\mathcal{R}}_{j}$ for some $j$) and $\eta_{ii^{\prime}}$ (for
$i,i^{\prime}\in\mathcal{C}$) depend only on the (uniform) distance between
the two voxels $i$ and $i^{\prime}$. For brevity, we still denote
$\rho_{|i^{\prime}-i|}$ by $\rho_{ii^{\prime}}$ and $\eta_{|i^{\prime}-i|}$ by
$\eta_{ii^{\prime}}$ where for $x\in\mathbb{Z}^{d}$, the notation $|x|$ stands
for the uniform norm. Our a priori is that the intra-correlation $\rho_{d}$ is
close to 1 for moderate distances $d$, meaning that close neighbors are
strongly (positively) correlated. For the local noise, we assume that the
spatial correlation function is such that $\eta_{0}=1$ and $\eta_{d}=0$ for
$d\geq p$. When $p=1$, this of course means that for any
$i,i^{\prime}\in\mathcal{C},i\neq i^{\prime}$, $\varepsilon_{i}(t)$ and
$\varepsilon_{i^{\prime}}(t)$ are uncorrelated.

######  Remarque 2.1.

We assume that the random variables are independent in time. This is not a too
restrictive hypothesis, in particular, if the random variables have long
memory hypothesis, after a wavelet decomposition, the random variables 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$.

###  2.1 Objectives and issues

Given a parcellation of the brain, the objective is to estimate inter-
correlations for each pair of parcel. The estimator has to be non parametric,
robust and fast. In particular we do not want to focus on the intra-
correlations, nor the estimations of variances and correlations of the
additional noises.

##  3 Inter-correlation estimators

###  3.1 Notation

Let $\mathbf{Y}_{1}=(Y_{1}(1),\ldots,Y_{1}(T))$ and
$\mathbf{Y}_{2}=(Y_{2}(1),\ldots,Y_{2}(T))$ denote two samples of length $T$
and let $\widehat{\mathrm{cov}}(\mathbf{Y}_{1},\mathbf{Y}_{2})$,
$\widehat{\mathrm{cor}}(\mathbf{Y}_{1},\mathbf{Y}_{2})$ and
$\widehat{\sigma}(\mathbf{Y}_{1})$ denote respectively the sample covariance
between $\mathbf{Y}_{1}$ and $\mathbf{Y}_{2}$, the sample correlation between
$\mathbf{Y}_{1}$ and $\mathbf{Y}_{2}$ and the standard deviation of
$\mathbf{Y}_{1}$.

In this section, we recall some known estimates of the inter-correlation
parameter $r_{jj^{\prime}}$ (that we simply denote by $r$ for brevity) for two
regions of interest ${\mathcal{R}}_{j}$ and ${\mathcal{R}}_{j^{\prime}}$. We
also introduce some new ones which can handle the problems occurring when
dealing with spatio-temporal data such as fMRI data.

To understand the limit of each estimator, we need some additional notation.
For any $j=1,\ldots,J$, we first define a $\nu$-neighborhood as a subset of
$n_{\nu}:=(2\nu+1)^{d}$ indices for which all indices are at distance
$\leq\nu$ from the center.

Then, for any $\nu$-neighborhood $\mathcal{V}$ ($\nu\geq 1$) we let

$\bar{\rho}_{\nu}=\frac{1}{n_{\nu}^{2}}\sum_{i,i^{\prime}\in\mathcal{V}}\rho_{|i^{\prime}-i|}\quad\mbox{
and
}\quad\bar{\eta}_{\nu}=\frac{1}{n_{\nu}^{2}}\sum_{i,i^{\prime}\in\mathcal{V}}\eta_{|i^{\prime}-i|}=\frac{1}{n_{\nu}^{2}}\sum_{\begin{subarray}{l}i,i^{\prime}\in\mathcal{V}\\\
|i^{\prime}-i|<p\end{subarray}}\eta_{|i^{\prime}-i|}.$

We also need the related notation where $\mathcal{V}$ corresponds to the
entire region of interest ${\mathcal{R}}_{j}$

$\bar{\rho}^{(j)}=\frac{1}{N_{j}^{2}}\sum_{i,i^{\prime}\in\mathcal{{\mathcal{R}}}_{j}}\rho_{|i^{\prime}-i|}\quad\mbox{
and
}\quad\bar{\eta}^{(j)}=\frac{1}{N_{j}^{2}}\sum_{\begin{subarray}{l}i,i^{\prime}\in\mathcal{{\mathcal{R}}}_{j}\\\
|i^{\prime}-i|<p\end{subarray}}\eta_{|i^{\prime}-i|}$

If $\varepsilon_{i}$ and $\varepsilon_{i^{\prime}}$ are spatially uncorrelated
for any $i,i^{\prime}\in\mathcal{C}$ ($i\neq i^{\prime}$), then
$\bar{\eta}_{\nu}=1/n_{\nu}$ and $\bar{\eta}^{(j)}=1/N_{j}$. Otherwise since
we assume that $\eta_{d}=0$ for $d\geq p$,
$\bar{\eta}_{\nu}=\mathcal{O}(1/n_{\nu})$ and
$\bar{\eta}^{(j)}=\mathcal{O}(1/N_{j})$. For moderate $\nu$, it is expected
that the intra-correlations satisfy $\bar{\rho}_{\nu}$ is close to 1. However,
$\bar{\rho}^{(j)}$ can be significantly far from 1, especially for large
regions of interest. Finally for more complex estimators developed hereafter,
we need the additional notation

$\bar{\rho}_{\nu,\delta}=\frac{1}{n_{\nu}^{2}}\sum_{i\in\mathcal{V},i^{\prime}\in\mathcal{V}^{\prime}}\rho_{|i^{\prime}-i|}$

for two $\nu$-neighborhoods $\mathcal{V},\mathcal{V}^{\prime}$ (from the same
region of interest) distant from $\delta$ (i.e.
$d(\mathcal{V},\mathcal{V}^{\prime})=\delta$). Again, if $\nu$ and $\delta$
are not too large, we can expect that $\bar{\rho}_{\nu,\delta}$ is close to 1.

From Section 2, we may derive the following calculations. Let
$j,j^{\prime}\in\\{1,\ldots,J\\}$ two different indices. Let
$i\in{\mathcal{R}}_{j},i^{\prime}\in{\mathcal{R}}_{j}\cup{\mathcal{R}}_{j^{\prime}}$.

$\mathrm{cov}(Y_{i}(t),Y_{i^{\prime}}(t))=\left\\{\begin{array}[]{lll}\sigma_{j}\sigma_{j^{\prime}}r_{jj^{\prime}}+\sigma_{\varepsilon}^{2}\eta_{|i^{\prime}-i|}+\sigma^{2}_{e}&\mbox{
if }i\in{\mathcal{R}}_{j}\mbox{ and
}i^{\prime}\in{\mathcal{R}}_{j^{\prime}}\\\
\sigma_{j}^{2}\rho_{|i^{\prime}-i|}+\sigma_{\varepsilon}^{2}\eta_{|i^{\prime}-i|}+\sigma^{2}_{e}&\mbox{
if }i,i^{\prime}\in{\mathcal{R}}_{j}.\end{array}\right.$ (2)

If $\bar{Y}_{j}(t)$ denotes the spatial average signal at time $t$ over the
region ${\mathcal{R}}_{j}$, i.e.
$\bar{Y}_{j}(t)=N_{j}^{-1}\sum_{i\in{\mathcal{R}}_{j}}Y_{i}(t)$, then for
$t=1,\ldots,T$

$\mathrm{Var}(\bar{Y}_{j}(t))=\sigma_{j}^{2}\bar{\rho}_{j}+\sigma_{\varepsilon}^{2}\bar{\eta}_{j}+\sigma_{e}^{2}.$
(3)

Since, we assume the independence in time of the signals, the local noise and
the global noise most of the results presented in the next sections are
straightforward. Without loss of generality, we intrinsically assume that
$d({\mathcal{R}}_{j},{\mathcal{R}}_{j^{\prime}})\geq p$ which ensures that for
any $i\in{\mathcal{R}}_{j}$ and $i^{\prime}\in{\mathcal{R}}_{j^{\prime}}$,
$\varepsilon_{i}(t)$ and $\varepsilon_{i^{\prime}}(t)$ are uncorrelated. This
slightly simplifies the numerator of $r^{\mbox{\scriptsize\sc\,ca}}$ and
$r^{\mbox{\scriptsize\sc\,ac}}$ defined below by (5) and (7).

###  3.2 Classical estimator : correlation of averages (method ca)

In order to increase the signal-to-noise ratio, classical methods in fMRI or
EEG use to average or convolve with Gaussian kernel the signal in space. The
aggregated correlation estimator corresponds to the classic estimator
considered for example in Achard et al. (2006) by taking the correlation
between the aggregated variables corresponding to a pair of given groups:

$\widehat{r}^{{\mbox{\scriptsize\sc\,ca}}}=\frac{\widehat{\mathrm{cov}}(\mathbf{\bar{Y}}_{j},\mathbf{\bar{Y}}_{j^{\prime}})}{{\widehat{\sigma}(\mathbf{\bar{Y}}_{j})}{\widehat{\sigma}(\mathbf{\bar{Y}}_{j^{\prime}})}}$
(4)

Using (2) and (3), it can be shown that
$\widehat{r}^{{\mbox{\scriptsize\sc\,ca}}}$ is a strongly consistent estimator
of $r^{\mbox{\scriptsize\sc\,ca}}$ as $T\to\infty$ where

$r^{\mbox{\scriptsize\sc\,ca}}=\frac{\sigma_{j}\sigma_{j^{\prime}}r+\sigma_{e}^{2}}{(\sigma_{j}^{2}\bar{\rho}^{(j)}+\sigma_{\varepsilon}^{2}\bar{\eta}^{(j)}+\sigma_{e}^{2})^{1/2}(\sigma_{j^{\prime}}^{2}\bar{\rho}^{(j^{\prime})}+\sigma_{\varepsilon}^{2}\bar{\eta}^{(j^{\prime})}+\sigma_{e}^{2})^{1/2}}.$
(5)

Even in absence of noise ($\sigma_{\varepsilon}=\sigma_{e}=0$),
$\widehat{r}^{\mbox{\scriptsize\sc\,ca}}$ has the drawback to bias the
estimation because the average of correlations is not equal to the correlation
of the averages. In Achard et al. (2011), we demonstrated that, even without
noise, the correlation of averages is multiplied by the number of points taken
in the average. So when the group of voxels in the brain have different sizes,
the largest ones have higher correlation simply because of this bias mainly
due to the fact that for large regions $\bar{\rho}^{(j)}$ can be far from 1.

###  3.3 Average of correlations (method ac)

This estimator defined by

$\widehat{r}^{{\mbox{\scriptsize\sc\,ac}}}=\frac{1}{N_{j}N_{j^{\prime}}}\sum_{i\in{\mathcal{R}}_{j},i^{\prime}\in{\mathcal{R}}_{j^{\prime}}}\widehat{\mathrm{cor}}(\mathbf{Y}_{i},\mathbf{Y}_{i^{\prime}}).$
(6)

clearly tends (almost surely) as $T\to\infty$ toward
$r^{\mbox{\scriptsize\sc\,ac}}$ given by

$r^{\mbox{\scriptsize\sc\,ac}}=\frac{\sigma_{j}\sigma_{j^{\prime}}r+\sigma_{e}^{2}}{(\sigma_{j}^{2}+\sigma_{\varepsilon}^{2}+\sigma_{e}^{2})^{1/2}(\sigma_{j^{\prime}}^{2}+\sigma_{\varepsilon}^{2}+\sigma_{e}^{2})^{1/2}}.$
(7)

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

$\widehat{r}^{\widetilde{{\mbox{\scriptsize\sc\,ac}}}}=\frac{1}{N_{j}N_{j^{\prime}}}\left(\sum_{i,i^{\prime}\in{\mathcal{R}}_{j}}\widehat{\mathrm{cor}}(\mathbf{Y}_{i},\mathbf{Y}_{i^{\prime}})\sum_{i,i^{\prime}\in{\mathcal{R}}_{j^{\prime}}}\widehat{\mathrm{cor}}(\mathbf{Y}_{i},\mathbf{Y}_{i^{\prime}})\right)^{1/2}\widehat{r}^{\mbox{\scriptsize\sc\,ac}}.$
(8)

The two estimators (6) and (8) have the important property to remove the size
effect (since when $\sigma_{\varepsilon}=\sigma_{e}=0$,
$r^{\mbox{\scriptsize\sc\,ac}}=r$). As it can be straightforwardly shown that
these two estimators tend to the same limit, we only focus in Sections 4, we
skip the estimator $\widehat{r}^{\widetilde{{\mbox{\scriptsize\sc\,ac}}}}$.

We can underline that both estimators are very sensitive to (local and global)
noise. Indeed, with the presence of local noise, the variance of noise appears
directly in the denominator which is decreasing the values of correlation
estimations. In the sequel of this section, we introduce first estimators to
be robust to the local noise using either local average or replicated data.
Then, we introduce estimators to compensate the global noise by using 2 other
decorrelated groups of variables, and then we will combine the two to define
estimators robust to both local and global noise.

###  3.4 Local correlation of averages (method ${\ell\mbox{\sc ca}}$)

The following estimator is defined by averaging correlations between average
of variables belonging to two different regions on small neighborhood inside
the regions. Its aim is to reduce the local noise. As the choice of
neighborhood to average and the choice of pairs of neighborhood are arbitrary,
we propose to repeat draws of small neighborhood and to take the average of
all computed correlations. For $b=1,\ldots,B$, let $\mathcal{V}_{j}^{(b)}$
(resp. $\mathcal{V}_{j^{\prime}}^{(b)}$) be a $\nu$-neighborhood of
${\mathcal{R}}_{j}$ (resp. ${\mathcal{R}}_{j^{\prime}}$), then we define

$\widehat{r}^{{\,\ell\mbox{\scriptsize\sc
ca}}}=\frac{1}{B}\sum_{b=1}^{B}\widehat{\mathrm{cor}}(\mathbf{\bar{Y}}_{\mathcal{V}_{j}^{(b)}},\mathbf{\bar{Y}}_{\mathcal{V}_{j^{\prime}}^{(b)}}).$
(9)

Since $d({\mathcal{R}}_{j},r_{j^{\prime}})\geq p$ then
$d(\mathcal{V}_{j}^{(b)},\mathcal{V}_{j^{\prime}}^{(b)})\geq p$ for any $b$
and $\widehat{r}^{{\,\ell\mbox{\scriptsize\sc ca}}}$ tends (almost surely) as
$T\to\infty$ towards $r^{\,\ell\mbox{\scriptsize\sc ca}}$

$r^{\,\ell\mbox{\scriptsize\sc
ca}}=\frac{\sigma_{j}\sigma_{j^{\prime}}r+\sigma_{e}^{2}}{(\sigma_{j}^{2}\bar{\rho}_{\nu}+\sigma_{e}^{2}\bar{\eta}_{\nu}+\sigma_{e}^{2})^{1/2}(\sigma_{j^{\prime}}^{2}\bar{\rho}_{\nu}+\sigma_{e}^{2}\bar{\eta}_{\nu}+\sigma_{e}^{2})^{1/2}}.$
(10)

When there is no global noise ($\sigma_{e}=0$), we observe that
$\widehat{r}^{\,\ell\mbox{\scriptsize\sc ca}}$ is more robust to the local
noise since $\bar{\eta}_{\nu}=\mathcal{O}(n_{\nu}^{-1})$. However, it is
important to choose moderate $\nu$ otherwise $\bar{\rho}_{\nu}$ can deviate
significantly from 1.

###  3.5 Replicates for correlations (method r)

In order to suppress the effect of local noise, we introduce a new estimator
computed using replicates within the same region and the estimator
$\widehat{r}^{\mbox{\scriptsize\sc\,ac}}$ is corrected by the correlation
between the chosen replicates. This estimator, denoted by
$\widehat{r}^{\mbox{\scriptsize\sc\,R}}$ (r for replicates), was first
introduced by Bergholm et al. (2010), in the context of image analysis. It is
defined by

$\widehat{r}^{{\mbox{\scriptsize\sc\,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^{\prime}}_{\beta}^{(b)}})}{\sqrt{\widehat{\mathrm{cor}}(\mathbf{Y}_{i_{1}^{(b)}},\mathbf{Y}_{i_{2}^{(b)}})\;\widehat{\mathrm{cor}}(\mathbf{Y}_{{i^{\prime}}_{1}^{(b)}},\mathbf{Y}_{{i^{\prime}}_{2}^{(b)}})}}$
(11)

where for $b=1,\ldots,B$, $i_{1}^{(b)}$ and $i_{2}^{(b)}$ belong to
${\mathcal{R}}_{j}$, and are expected to verify
$\rho_{i_{1}^{(b)}i_{2}^{(b)}}\approx 1$. Respectively,
${i^{\prime}}_{1}^{(b)}$ and ${i^{\prime}}_{2}^{(b)}$ belong to
${\mathcal{R}}_{j^{\prime}}$ and are expected to verify
$\rho_{{i^{\prime}}_{1}^{(b)}{i^{\prime}}_{2}^{(b)}}\approx 1$.

Let us detail the convergence of $\widehat{r}^{\mbox{\scriptsize\sc\,R}}$. We
assume that any $b=1,\ldots,B$,
$|i_{2}^{(b)}-i_{1}^{(b)}|=|{i^{\prime}_{2}}^{(b)}-{i^{\prime}_{1}}^{(b)}|=\delta$
for some fixed integer $\delta\geq p$. Since

$\frac{1}{4}\sum_{\alpha,\beta=1}^{2}\widehat{\mathrm{cor}}(\mathbf{Y}_{i_{\alpha}^{(b)}},\mathbf{Y}_{{i^{\prime}}_{\beta}^{(b)}})\stackrel{{\scriptstyle
a.s.}}{{\rightarrow}}\frac{\sigma_{j}\sigma_{j^{\prime}}r_{jj^{\prime}}+\sigma_{e}^{2}}{\left(\sigma_{j}^{2}+\sigma_{\varepsilon}^{2}+\sigma_{e}^{2}\right)^{1/2}\;\left(\sigma_{j^{\prime}}^{2}+\sigma_{\varepsilon}^{2}+\sigma_{e}^{2}\right)^{1/2}}$

and

$\widehat{\mathrm{cor}}(\mathbf{Y}_{i_{1}^{(b)}},\mathbf{Y}_{i_{2}^{(b)}})\stackrel{{\scriptstyle
a.s.}}{{\rightarrow}}\frac{\sigma_{j}^{2}\rho_{\delta}+\sigma_{\varepsilon}^{2}\eta_{\delta}+\sigma_{e}^{2}}{\sigma_{j}^{2}+\sigma_{\varepsilon}^{2}+\sigma_{e}^{2}},$

and since $\eta_{\delta}=0$, we obtain that
$\widehat{r}^{\mbox{\scriptsize\sc\,R}}$ converges almost surely towards
$r^{\mbox{\scriptsize\sc\,R}}$ as $T\to\infty$ where

$r^{\mbox{\scriptsize\sc\,R}}=\frac{\sigma_{j}\sigma_{j^{\prime}}r+\sigma_{e}^{2}}{(\sigma_{j}^{2}\rho_{\delta}+\sigma_{e}^{2})^{1/2}(\sigma_{j^{\prime}}^{2}\rho_{\delta}+\sigma_{e}^{2})^{1/2}}$
(12)

When $\sigma_{e}=0$ and then choosing $\delta=p$ (so $1$ if the local noise is
spatially independent) leads to
$r^{\mbox{\scriptsize\sc\,R}}=\frac{r}{\rho_{p}}$ and we may hope that in this
situation, $\rho_{p}$ is very close to 1. In other words,
$\widehat{r}^{\mbox{\scriptsize\sc\,R}}$ is an estimator robust to the size of
the regions of interest and robust to a local noise.

In the case of low signal-to-noise ratio, the quantity in the denominator of
(11) may be negative and the square root becomes undefined. We propose then to
combine the local averaging on small neighborhood and the introduction of
replicates to compensate for local noise corruption of signals.

###  3.6 Local average of replicates (method ${\ell\mbox{\sc r}}$)

This estimator consists in combining the idea of replicates with the one
consisting in averaging locally the signals.

$\widehat{r}^{{\,\ell\mbox{\scriptsize\sc
R}}}=\frac{1}{B}\sum_{b=1}^{B}\frac{\frac{1}{4}\sum_{\alpha,\beta=1}^{2}\widehat{\mathrm{cor}}(\mathbf{\bar{Y}}_{\mathcal{V}_{j_{\alpha}}^{(b)}},\mathbf{\bar{Y}}_{\mathcal{V}_{j^{\prime}_{\beta}}^{(b)}})}{\sqrt{\widehat{\mathrm{cor}}(\mathbf{\bar{Y}}_{\mathcal{V}_{j_{1}}^{(b)}},\mathbf{\bar{Y}}_{\mathcal{V}_{j_{2}}^{(b)}})\;\widehat{\mathrm{cor}}(\mathbf{\bar{Y}}_{\mathcal{V}_{j^{\prime}_{1}}^{(b)}},\mathbf{\bar{Y}}_{\mathcal{V}_{j^{\prime}_{2}}^{(b)}})}}.$
(13)

Using arguments developed in Sections 3.4 and 3.5 then if we assume that the
distance between $\nu$-neighborhoods involved in (13) is at least $p$, then it
can be shown $\widehat{r}^{\,\ell\mbox{\scriptsize\sc R}}$ is a strongly
consistent estimator of $r^{\,\ell\mbox{\scriptsize\sc R}}$ defined by

$r^{\,\ell\mbox{\scriptsize\sc
R}}=\frac{\sigma_{j}\sigma_{j^{\prime}}r+\sigma_{e}^{2}}{(\sigma_{j}^{2}\bar{\rho}_{\nu,\delta}+\sigma_{e}^{2})^{1/2}(\sigma_{j^{\prime}}^{2}\bar{\rho}_{\nu,\delta}+\sigma_{e}^{2})^{1/2}}$
(14)

When $\sigma_{e}=0$, then letting $\delta=p$ reduces
$r^{\,\ell\mbox{\scriptsize\sc R}}$ to $\frac{r}{\bar{\rho}_{\nu,p}}$.

###  3.7 Use of a priori disconnected regions (method d based on differences)

We now present an estimator which handles the problem of global noise, based
on the selection of two extra regions that are uncorrelated. This assumption
is realistic in the context of fMRI data where we are interested in the
correlations between cortical regions. Indeed, the volume recorded by the
scanning may be bigger, and the definition of extra regions is possible. In
this case, if we want to estimate the correlation between two regions, say
${\mathcal{R}}_{j}$ and ${\mathcal{R}}_{j^{\prime}}$, the idea is to assume
the extra information that neither ${\mathcal{R}}_{j}$ nor
${\mathcal{R}}_{j^{\prime}}$ is connected to some regions ${\mathcal{R}}_{k}$
and ${\mathcal{R}}_{k^{\prime}}$ where
$k,k^{\prime}\in\\{1,\ldots,J\\}\setminus\\{j,j^{\prime}\\}$ with $k\neq
k^{\prime}$ and we assume in addition that ${\mathcal{R}}_{k}$ and
${\mathcal{R}}_{k^{\prime}}$ are also not connected. Then, we propose the
following estimator: for $b=1,\ldots,B$ let $i^{(b)}$, ${i^{\prime}}^{(b)}$,
$k^{(b)}$ and ${k^{\prime}}^{(b)}$ be voxels of ${\mathcal{R}}_{j}$,
${\mathcal{R}}_{j^{\prime}}$, ${\mathcal{R}}_{k}$ and
${\mathcal{R}}_{k^{\prime}}$.

$\widehat{r}^{\mbox{\scriptsize\sc\,D}}=\frac{1}{B}\sum_{b=1}^{B}\widetilde{\mathrm{cor}}(\mathbf{Y}_{i^{(b)}},\mathbf{Y}_{{i^{\prime}}^{(b)}};\mathbf{Y}_{k^{(b)}},\mathbf{Y}_{{k^{\prime}}^{(b)}}),$
(15)

where for four vectors $\mathbf{Y}_{1}$, $\mathbf{Y}_{2}$, $\mathbf{Y}_{3}$
and $\mathbf{Y}_{4}$ (with same length)

$\widetilde{\mathrm{cor}}(\mathbf{Y}_{1},\mathbf{Y}_{2};\mathbf{Y}_{3},\mathbf{Y}_{4})=\frac{\widehat{\mathrm{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})}$
(16)

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

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

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

From a theoretical point of view, we can show that almost surely
$\widehat{r}^{\mbox{\scriptsize\sc\,D}}$ converges almost surely as
$T\to\infty$ towards $r^{\mbox{\scriptsize\sc\,D}}$ given by

$r^{\mbox{\scriptsize\sc\,D}}=\frac{\sigma_{j}\sigma_{j^{\prime}}r+\tau_{jj^{\prime}}}{\left(\sigma_{j}^{2}+2\sigma_{\varepsilon}^{2}+\tau_{j}\right)^{1/2}\;\left(\sigma_{j^{\prime}}^{2}+2\sigma_{\varepsilon}^{2}+\tau_{j^{\prime}}\right)^{1/2}}$

where

$\tau_{jj^{\prime}}=\sigma_{j}\sigma_{k^{\prime}}r_{jk^{\prime}}+\sigma_{j^{\prime}}\sigma_{k}r_{j^{\prime}k}\quad\mbox{
and
}\quad\tau_{j}=-\sigma_{j}\sigma_{k}r_{jk}-\sigma_{j}\sigma_{k^{\prime}}r_{jk^{\prime}}+\sigma_{k}\sigma_{k^{\prime}}r_{kk^{\prime}}.$

Since we assume a priori that the regions ${\mathcal{R}}_{k}$ and
${\mathcal{R}}_{k^{\prime}}$ are not connected each other and to the other
regions then $\tau_{j}=\tau_{j^{\prime}}=\tau_{jj^{\prime}}=0$ which reduces
the previous expression to

$r^{\mbox{\scriptsize\sc\,D}}=\frac{\sigma_{j}\sigma_{j^{\prime}}r}{\left(\sigma_{j}^{2}+2\sigma_{\varepsilon}^{2}\right)^{1/2}\;\left(\sigma_{j^{\prime}}^{2}+2\sigma_{\varepsilon}^{2}\right)^{1/2}}$
(17)

which in case of absence of local noise is nothing else than $r$.

###  3.8 Combinations of previous estimators: methods ${\ell\mbox{\sc d}}$,
rd, ${\ell\mbox{\sc rd}}$

The three following estimators extend respectively the methods ${\ell\mbox{\sc
ca}}$, r, ${\ell\mbox{\sc r}}$ using the idea developed above to get rid of
the global noise. We use the notation presented in the previous sections.
Theoretical results are also derived along similar lines.

### Local averages and use of a priori disconnected regions: method
${\ell\mbox{\sc d}}$

This estimator given by

$\widehat{r}^{\,\ell\mbox{\scriptsize\sc
D}}=\frac{1}{B}\sum_{b=1}^{B}\widetilde{\mathrm{cor}}(\mathbf{\bar{Y}}_{\mathcal{V}_{j}^{(b)}},\mathbf{\bar{Y}}_{\mathcal{V}_{j^{\prime}}^{(b)}};\mathbf{\bar{Y}}_{\mathcal{V}_{k}^{(b)}},\mathbf{\bar{Y}}_{\mathcal{V}_{k^{\prime}}^{(b)}})$
(18)

is a strongly consistent estimator of $r^{\,\ell\mbox{\scriptsize\sc D}}$
given by

$r^{\,\ell\mbox{\scriptsize\sc
D}}=\frac{\sigma_{j}\sigma_{j^{\prime}}r}{(\sigma_{j}^{2}+2\sigma_{\varepsilon}^{2})^{1/2}(\sigma_{j^{\prime}}^{2}+2\sigma_{\varepsilon}^{2})^{1/2}}.$
(19)

### Replicates and use of a priori disconnected regions: method rd

This estimator given by

$\widehat{r}^{\,\mbox{\scriptsize\sc
RD}}=\frac{1}{B}\sum_{b=1}^{B}\frac{\frac{1}{4}\sum_{\alpha,\beta=1}^{2}\widetilde{\mathrm{cor}}(\mathbf{Y}_{i_{\alpha}^{(b)}},\mathbf{Y}_{{i^{\prime}}_{\beta}^{(b)}};\mathbf{Y}_{{i}_{3}^{(b)}},\mathbf{Y}_{{i}_{4}^{(b)}})}{\sqrt{\widetilde{\mathrm{cor}}(\mathbf{Y}_{i_{1}^{(b)}},\mathbf{Y}_{i_{2}^{(b)}};\mathbf{Y}_{k^{(b)}},\mathbf{Y}_{{k^{\prime}}^{(b)}})\;\widetilde{\mathrm{cor}}(\mathbf{Y}_{{i^{\prime}}_{1}^{(b)}},\mathbf{Y}_{{i^{\prime}}_{2}^{(b)}};\mathbf{Y}_{k^{(b)}},\mathbf{Y}_{{k^{\prime}}^{(b)}})}}$
(20)

is a strongly consistent estimator of $r^{\,\mbox{\scriptsize\sc RD}}$ given
by

$r^{\,\mbox{\scriptsize\sc RD}}=\frac{r}{\bar{\rho}_{\delta}}.$ (21)

### Replicates, local averages and use of a priori disconnected regions:
method ${\ell\mbox{\sc rd}}$

This estimator given by

$\widehat{r}^{\,\ell\mbox{\scriptsize\sc
RD}}=\frac{1}{B}\sum_{b=1}^{B}\frac{\frac{1}{4}\sum_{\alpha,\beta=1}^{2}\widetilde{\mathrm{cor}}(\mathbf{\bar{Y}}_{\mathcal{V}_{j_{\alpha}}^{(b)}},\mathbf{\bar{Y}}_{\mathcal{V}_{{j^{\prime}_{\beta}}}^{(b)}};\mathbf{\bar{Y}}_{\mathcal{V}_{k}^{(b)}},\mathbf{\bar{Y}}_{\mathcal{V}_{k^{\prime}}^{(b)}})}{\sqrt{\widetilde{\mathrm{cor}}(\mathbf{\bar{Y}}_{\mathcal{V}_{j_{1}}^{(b)}},\mathbf{\bar{Y}}_{\mathcal{V}_{{j}_{2}}^{(b)}};\mathbf{\bar{Y}}_{\mathcal{V}_{k}^{(b)}},\mathbf{\bar{Y}}_{\mathcal{V}_{k^{\prime}}^{(b)}})\;\widetilde{\mathrm{cor}}(\mathbf{\bar{Y}}_{\mathcal{V}_{j_{1}^{\prime}}^{(b)}},\mathbf{\bar{Y}}_{\mathcal{V}_{{j^{\prime}}_{2}}^{(b)}};\mathbf{\bar{Y}}_{\mathcal{V}_{k}^{(b)}},\mathbf{\bar{Y}}_{\mathcal{V}_{k^{\prime}}^{(b)}})}}$
(22)

is a strongly consistent estimator of $r^{\,\ell\mbox{\scriptsize\sc RD}}$
given by

$r^{\,\ell\mbox{\scriptsize\sc RD}}=\frac{r}{\bar{\rho}_{\nu,\delta}}.$ (23)

##  4 Simulations

In this section we investigate the finite sample properties of the estimators
of $r$ according to different settings. We focus on the planar case ($d=2$)
and only on two regions of interest, say ${\mathcal{R}}_{j}$ and
${\mathcal{R}}_{j^{\prime}}$ which contain respectively $20^{2}$ and $40^{2}$
voxels. To implement the estimators $\widehat{r}^{\bullet}$ for
$\bullet={\mbox{\sc d}},{\ell\mbox{\sc d}},{\mbox{\sc rd}},{\ell\mbox{\sc
rd}}$ we also assume to observe the signals in two other regions
${\mathcal{R}}_{k}$ and ${\mathcal{R}}_{k^{\prime}}$ which are disconnected
each other and to the two regions of interest. We assume that the signal, the
local noise and the global noise follow Gaussian distribution. The variances
are set to the following values: $\sigma_{j}=1$, $\sigma_{j^{\prime}}=2$,
$\sigma_{k}=\sigma_{k^{\prime}}=1$. Throughout the simulations the parameter
$r=r_{jj^{\prime}}$ is fixed to 0.6, the sample size $T$ is set to 1000 and we
set $\eta_{k}=0$ for $k\geq 1$ (i.e. $\varepsilon_{i}(t)$ and
$\varepsilon_{i^{\prime}}(t)$ are uncorrelated). Finally, we consider the
following spatial function for the intra-correlation:
$\rho_{|i^{\prime}-i|}=\max(1-|i^{\prime}-i|/K_{\max},r_{\min})$. In
particular we consider two different models for the intra-correlation:

  * •

Model 1: $K_{\max}=300$, $r_{\min}=0.9$.

  * •

Model 2: $K_{\max}=100$, $r_{\min}=0.6$.

Model 1 is model ling a region of interest where the voxels are well connected
while model 2 will be used to show the problems that could occur when two
voxels inside a common region are not as connected as they should be. Figure 1
shows the intra-correlation matrix for both these models.

(a) Model 1 (b) Model 2

Figure 1: Intra-correlation matrices for the models 1 and 2 and for the
largest region containing $40^{2}$ voxels. The bottom-left square corresponds
to the intra-correlation matrix of the smallest region containing $20^{2}$
voxels.

###  4.1 Influence of the intra-correlation in absence of noise

We first study the properties of the estimators in absence of noise
($\sigma_{\varepsilon}=\sigma_{e}=0$). Figure 2 shows boxplots based on $500$
replications of the general model (1) with parameters described previously.
Figure 2 (with $\sigma_{\varepsilon}=\sigma_{e}=0$) shows that when the intra-
correlation within each region ${\mathcal{R}}_{j}$ and
${\mathcal{R}}_{j^{\prime}}$ is high for any pair of voxels, then the
estimates of $r$ are quite satisfactory even if we can already observe a
slight bias for the method ca. However, when the model 2 is considered, the
bias of the latter method is significant. It is to be noticed that in terms of
variance the methods based on replicates (i.e. methods r, ${\ell\mbox{\sc
r}}$) are as efficient as the methods ac and ${\ell\mbox{\sc ca}}$. The
methods based on the a priori knowledge of two other disconnected regions
(methods based on “differences”, i.e. the methods d, ${\ell\mbox{\sc d}}$, rd
and ${\ell\mbox{\sc rd}}$) have higher dispersion than the other ones.
Finally, we also observe that for the “complex” methods ${\ell\mbox{\sc r}}$
and ${\ell\mbox{\sc rd}}$ based on local averages generate a more important
bias than other methods. This bias is clearly smaller than the one observed
for the method ca and illustrates the following fact: when we suspect that the
intra-correlation matrix is not very well concentrated on the diagonal then
the size of the neighborhood ($\nu$ in the definition of the estimators)
should be chosen sufficiently small.

Figure 2: Boxplots of estimates of the inter-correlation parameter $r$ based
on 500 replications of the general model (1) for the models 1 and 2 of the
intra-correlation matrix. Situations with no noise, local noise or global
noise are considered with different levels of signal to noise ratio.

###  4.2 Influence of the local noise and global noise

We now turn to study the influence of a local noise or a global noise are also
illustrated in Figure 2. The variances of the local noise or the global noise
are parameterized by fixing the signal-to-noise ratio: when for instance
$\sigma_{e}=0$, we fixed $\sigma_{\varepsilon}^{2}$ as follows

$\mathrm{SNR}_{\varepsilon}=10\log_{10}\left(\frac{\min(\sigma_{j}^{2},\sigma_{j^{\prime}}^{2})}{\sigma_{\varepsilon}^{2}}\right)\Leftrightarrow\sigma_{\varepsilon}^{2}=10^{-\mathrm{SNR}_{\varepsilon}/10}\min(\sigma_{j}^{2},\sigma_{j^{\prime}}^{2}).$

When $\sigma_{\varepsilon}\neq 0$ and $\sigma_{e}=0$, as expected, Figure 2
shows that the methods based on replicates are able to estimate correctly the
inter-correlation parameter $r$ and these methods remain efficient whatever
the value of the $\mathrm{SNR}$ and for both models of the intra-correlation
matrix. The methods ac, ${\ell\mbox{\sc ca}}$, d and ${\ell\mbox{\sc d}}$ are
strongly affected by this additional local noise and exhibit a high negative
bias. As explained in Section 3.2, the method ca which averages the signals in
the regions ${\mathcal{R}}_{j}$ and ${\mathcal{R}}_{j^{\prime}}$ is able to
reduce the effect of the local noise.

Now, focusing on Figure 2 when $\sigma_{\varepsilon}=0$ and $\sigma_{e}\neq
0$, which corresponds to situations where there is only a global noise, we
observe that the higher the $\mathrm{SNR}_{e}$ the higher the positive bias
for all the methods except the ones which were expected to handle the global
noise, namely the methods d, ${\ell\mbox{\sc d}}$, rd and ${\ell\mbox{\sc
rd}}$. It is also to be noticed that the dispersion of all the methods do not
seem to be very affected by this extra noise.

##  5 Conclusion

In this paper we illustrate the effect of averaged data on estimations of
correlation when two noises are present, local and global noise. The use of
the classical correlation of averages is perturbed by the presence of these
noises in addition to the presence of within correlations. We proposed
alternative estimators including correction terms to compensate the intra-
correlations, local and global noises. The performances of these estimators
are illustrated on simulations.

## References

  * Achard et al. (2011) Achard, S., Coeurjolly, J., Marcillaud, R., Richiardi, J., 2011\.  fMRI functional connectivity estimators robust to region size bias, in: IEEE Workshop on Statistical Signal Processing, SSP2011, Nice, France.
  * Achard et al. (2006) 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.
  * Bergholm et al. (2010) 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.
  * Clifford et al. (1989) Clifford, P., Richardson, S., Hemon, D., 1989.  Assessing the significance of the correlation between two spatial processes.  Biometrics 45, 123–134.
  * Donner and Eliasziw (1991) Donner, A., Eliasziw, M., 1991\.  Methodology for inferences concerning familial correlations: a review.  J Clin Epidemio 44, 449–455.
  * Gunst (1995) Gunst, R.F., 1995.  Estimating spatial correlations from spatial-temporal meteorological data.  Journal of climate 8, 2454–69.
  * Liebhold and Sharov (1998) Liebhold, A., Sharov, A., 1998\.  Testing for correlation in the presence of spatial autocorrelation in insect count data.  Population and Community Ecology for Insect Management and Conservation , 1–8.
  * Ostroff (1993) Ostroff, C., 1993.  Comparing correlations based on individual-level and aggregated data.  Journal of Applied Psychology 78, 569.
  * Rosner et al. (1977) Rosner, B., Donner, A., Hennekens, C., 1977.  Estimation of interclass correlation from familial data.  Applied statistics 26, 179–187.
  * Student (1914) Student, 1914.  The elimination of spurious correlation due to position in time and space.  Biometrika 10, 179–181.
  * Termenon et al. (2016) 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.
  * Vul et al. (2009) Vul, E., Harris, C., Winkielman, P., Pashler, H., 2009\.  Puzzlingly high correlations in fmri studies of emotion, personality, and social cognition.  Perspectives on Psychological Science 4, 274–290.
  * Ye et al. (2011) Ye, J., Lazar, N., Li, Y., 2011.  Sparse geostatistical analysis in clustering fMRI time series.  Journal of Neuroscience Methods 199, 336–345.