# A Mixed Model Approach for Estimating Regional Functional Connectivity from
Voxel-level BOLD Signals

Chao<EMAIL_ADDRESS>[ Chau B.
<EMAIL_ADDRESS>[ Sophie Achardlabel=e3]sophie.achard@univ-
grenoble-alpes.fr [ Wendy<EMAIL_ADDRESS>[ Alexander
Petersen<EMAIL_ADDRESS>Correspondence to: Alexander Petersen,
Department of Statistics, Brigham Young University, Provo, Utah, 84602, USA[
Department of Statistics and Applied Probability, University of California,
Santa Barbara, Santa Barbara, CA, USApresep=, ]e1,e2,e4  Univ. Grenoble Alpes,
CNRS, Inria, Grenoble INP, LJK, F-38000, Grenoble, Francepresep=, ]e3
Department of Statistics, Brigham Young University, Provo, UT, USApresep=, ]e5

###### Abstract

Resting state brain functional connectivity quantifies the similarity between
brain regions, each of which consists of voxels at which dynamic signals are
acquired via neuroimaging techniques such as blood-oxygen-level-dependent
signals in functional magnetic resonance imaging. Pearson correlation and
similar metrics have been adopted by neuroscientists to estimate inter-
regional connectivity, usually after averaging of signals within regions.
However, dependencies between signals within each region and the presence of
noise could contaminate such inter-regional correlation estimates. We propose
a mixed-effects model with a novel covariance structure that explicitly
isolates the different sources of variability in the observed BOLD signals,
including correlated regional signals, local spatiotemporal variability, and
measurement error. Methods for tackling the computational challenges
associated with restricted maximum likelihood estimation will be discussed.
Large sample properties are discussed and used for uncertainty quantification.
Simulation results demonstrate that the parameters of the proposed model
parameters can be accurately estimated and is superior to the Pearson
correlation of averages in the presence of spatiotemporal noise. The proposed
model is also applied to a real data set of BOLD signals collected from rats
to construct individual brain networks.

functional brain connectivity,

mixed-effects model,

spatiotemporal modeling,

correlation networks,

###### keywords:

, , , and

##  1 Introduction

In recent years, the rapid advancement and increasing accessibility of
neuroimaging techniques, such as electroencephalogram (EEG) scans and
functional magnetic resonance imaging (fMRI) , data sets that record dynamic
brain activity are becoming widely available to scientists and medical
practitioners. The rich volume of previously unattainable data sets has
consequently catalyzed broad interest in modeling and estimating functional
brain connectivity, which is of paramount importance in shedding light on the
evolution of pathologies such as neurodegenerative diseases or consciousness
disorders.

A brain can be considered as a network. Thus, it is not surprising that
network-based approaches are prevalent in connectivity modeling (Van Den
Heuvel and Pol, 2010; Meskaldji et al., 2011; Kaiser, 2011). The nodes and
edges in a brain network represent brain regions and connections,
respectively. An edge is identified by a high degree of some similarity metric
for the brain activity signals, which is usually quantified empirically by the
Pearson correlation (Achard et al., 2006; Zalesky, Fornito and Bullmore, 2012;
Becq et al., 2020). However, challenges arise in preprocessing and conducting
robust and reproducible analyses with the massive data volumes, particularly
in the case of fMRI data sets. Signals collected from human brains are
frequently modeled as realizations of random fields that traverse across both
the spatial and temporal domains (Achard et al., 2011; Achard and Gannaz,
2019). These signals are usually of noisy nature due to the inherent
properties of brain activity and the errors incurred by the measuring
instrument (Chaimow et al., 2018; Lohmann et al., 2018; Park and Fiecas,
2022). In addition, the spatial dimension is usually high compared to the
temporal dimension, thus posing difficulties in conducting efficient and
reproducible analyses.

Various methods have been proposed to overcome these challenges. Achard et al.
(2006) applied a wavelet transformation to obtain brain networks based on
frequency-dependent correlations between regional fMRI time series. Machine
learning techniques, in particular, feature embedding and clustering are
investigated to classify and characterize changes in brain dynamics due to
pathology or cognitive state changes (Richiardi et al., 2013). Termenon et al.
(2016) studied the relationship between the reproducibility of brain networks
and the subject count and fMRI scan length using the large test-retest (TRT)
resting-state fMRI data set from the Human Connectome Project (HCP). In
addition, Petersen et al. (2016) proposed to view the path lengths of brain
networks as functions of network density; functional principal component
analysis (fPCA) is then applied and covariates are regressed on functional
principal component scores to study the association between connectivity and
subjects’ age, episodic memory, and executive function. It should be
emphasized that, while some of these existing methods can localize
connectivity to the subject level (Petersen et al., 2016; Bowman, 2007), the
analysis phase usually requires group level information.

Functional connectivity quantification and estimation is a primary goal of
these neuroscience inquiries. Following the characterization of Van Den Heuvel
and Pol (2010), functional connectivity is the dependency of simultaneous
neuronal activation patterns of anatomically separated brain regions.
Consequently, rather than directly modeling signals at the voxel level, the
vast majority of existing studies aggregate signals within certain and use an
intuitive quantification of dependence between two regions, such as the
Pearson correlation or partial correlation. These approaches do not consider
the spatiotemporal noise or voxel level fluctuations that are inherent to the
BOLD signals (Achard and Gannaz, 2019).

In brain activation studies, as opposed to connectivity studies, spatial and
spatiotemporal models have already been proposed to quantify brain regions
activated under a certain task (Taylor and Worsley, 2007; Bowman, 2007). These
models have been applied in fMRI (McIntosh, Chau and Protzner, 2004; Woolrich
et al., 2004; Zhang et al., 2014), but are not directly applicable for
constructing a single brain network from an individual scan taken at rest,
which is the goal of this paper. As motivation, we will use an anatomical
parcellation of a rat brain that contains $51$ regions, with the aim to
quantify and parameterize regional functional connectivity within a model for
voxel-level signals to separately construct an individual brain network from
each rat’s fMRI scan. This will allow us to give an estimate of the
classically defined _inter-regional_ , or _long-range_ connectivity network.
Another key measurement of connectivity, termed the _intra-regional_
connectivity, measures the dependence within a particular brain region. While
we will emphasize estimation of inter-regional connectivity, as it is of
primary interest in a large body of neuroimaging studies (Eickhoff, Yeo and
Genon, 2018; Moghimi et al., 2022), intra-regional connectivity will play a
key role in the estimation process.

With this purpose in mind, we highlight the main contributions of this paper.
First, we propose a novel spatiotemporal statistical model for the BOLD
signals at the voxel level. In principle, the collection of BOLD signals is
considered to be a realization of a spatiotemporal process with a carefully
designed covariance structure that explicitly models the different sources of
variability leading to spatiotemporal dependency of the fMRI BOLD signals. The
inter-regional connectivity between two arbitrary regions is quantified by a
primal model parameter, which needs to be estimated along with a collection of
auxiliary parameters that characterize the contributions of other underlying
variability sources. The model will be presented in relatively general terms,
after which an inferential pipeline will be developed in the Gaussian setting.
As a second main contribution, we devise an efficient two-stage strategy to
estimate model parameters as computations associated with Gaussian processes
are usually prohibitive. In particular, the proposed approach leverages large
scale parallel computing with GPU accelerated matrix operations. In addition,
we outline a justification for the use of large sample properties of the
inter-regional connectivity estimates to provide sampling error quantification
for network construction.

The remainder of the paper is arranged as follows. In Section 2, we present
the model construction, parameterization, and estimation procedure. In
particular, we propose the two-stage estimation procedure in Section 2.3.
Assumptions and results on the asymptotic behavior of the estimators are
presented in Section 2.4. Section 3 provides a set of simulation studies that
demonstrates the empirical performance of our estimators on synthetic data. In
Section 4, our model is applied to real rat data sets and brain connectivity
is quantified using the new estimation procedure and standard errors provided
by the asymptotic theory. Section 5 provides a discussion of the results and
potential extensions.

##  2 Methodology

For an individual brain volume, let $\mathcal{B}\subset\mathbb{R}^{3}$ denote
the connected set of spatial coordinates which is to be sampled. $\mathcal{B}$
is divided into spatially disjoint and contiguous regions $\mathcal{R}_{j}$,
$j=1,\ldots,J$, that is, $\mathcal{B}=\cup_{j=1}^{J}\mathcal{R}_{j}$ and
$\mathcal{R}_{j}\cap\mathcal{R}_{j^{\prime}}=\emptyset$, $j\neq j^{\prime}$.
Typically, studies of resting-state functional connectivity take a network-
based approach to model the functional connectivity properties of the brain.
In this network, the nodes correspond to the different brain regions
$\mathcal{R}_{j},$ while connections between regions are determined by a
collection of unknown parameters $\rho_{jj^{\prime}},$ $1\leq j<j^{\prime}\leq
J,$ that are often taken to be correlation values. For now, let us assume that
$|\rho_{jj^{\prime}}|\leq 1$ represent correlations, so that we are dealing
with a correlation network; that is, the edges may be considered binary
according to whether $\rho_{jj^{\prime}}=0$, in which case there is no edge,
or $\rho_{jj^{\prime}}\neq 0,$ indicating the presence of an edge between
regions $j$ and $j^{\prime}$. Alternatively, the raw value
$\rho_{jj^{\prime}}$ can represent the strength of the edge between these
regions in a weighted network.

In principle, brain activity evolves over a continuum, though it is not
observed directly, continuously, or cleanly. Hence, to begin, denote by
$Y_{j}(v,t)$ a latent BOLD signal at location $v\in\mathcal{R}_{j}$,
$j=1,\dots,J$, and at time $t\in[0,\mathcal{T}]$, where $\mathcal{T}$ is the
duration of the fMRI scan. The observed data are

$X_{jlm}=Y_{j}(v_{jl},t_{m})+\epsilon_{jlm},$ (2.1)

being sampled on a grid of voxels $v_{jl}\in\mathcal{R}_{j},$ $1\leq l\leq
L_{j},$ and a sequence of time points
$0=t_{0}<t_{1}<\cdots<t_{M}=\mathcal{T},$ where $\epsilon_{jlm}$ are assumed
to be independent and identically distributed zero-mean additive errors with
variance $\sigma^{2}$. At this point, it should be emphasized that, without
further assumptions and modeling, a principled approach to statistical
inference to reconstruct a network from the data would be difficult to
justify, even if the latent processes $Y_{j}(v,t)$ were available to the
scientist. This is due to the disparity between the resolution of the data,
measured at unique voxels, and the parameters that are specified at the
regional level. Put another way, the role of $\rho_{jj^{\prime}},$ a regional
parameter, on the joint distribution of the spatiotemporal fields $Y_{j}$ and
$Y_{j^{\prime}}$ is not yet specified. Nevertheless, it has become commonplace
in neuroscience to construct a network using the so-called “correlation of
averages” that will be described next, along with an analysis of its
statistical properties under a basic preliminary model for the fields $Y_{j}$.

###  2.1 Correlation of Averages

For a pair of regions, the common approach to quantify an empirically-based
network is to first average voxel-level signals across space to obtain a mean
signal for each region, that is,

$\bar{X}_{jm}=L_{j}^{-1}\sum_{l=1}^{L_{j}}X_{jlm}.$ (2.2)

The use of a single signal per region simplifies not only calculations by
reducing the data dimension, but also the modeling that is required, since
voxel-specific fluctuations are not taken into account. Letting
$\tilde{\mu}_{j}$ and $\tilde{\zeta}^{2}_{j}$ be the empirical mean and
variance of $\bar{X}_{jm}$ across $m$ (time), the usual Pearson correlation
type quantification of connectivity is

$\hat{\rho}^{\textrm{CA}}_{jj^{\prime}}=\frac{M^{-1}\sum_{m=1}^{M}\left(\bar{X}_{jm}-\tilde{\mu}_{j}\right)\left(\bar{X}_{j^{\prime}m}-\tilde{\mu}_{j^{\prime}}\right)}{\tilde{\zeta}_{j}\tilde{\zeta}_{j^{\prime}}}.$
(2.3)

We term the estimator in (2.3) the “correlation of averages” estimator as it
averages signals by regions prior to the calculation of the Pearson
correlation. Although this estimator is commonly used, only recently has more
rigorous attention been given to the corresponding estimand. Ignoring the
contribution of noise, which is somewhat mitigated by the averaging operation
if the number of voxels is large, (2.3) would be a reasonable estimator if one
assumes the correlation between two average regional signals is constant in
time. Indeed, in those examples in the literature where a theoretical analysis
if $\hat{\rho}^{\textrm{CA}}$ is provided, said analysis is typically done
under the hypothetical framework of a single signal observed per region
(Afyouni, Smith and Nichols, 2019; Azevedo et al., 2022), so that averaging
across voxels is merely viewed as a preprocessing step without assessing the
additional uncertainty that it induces. However, it is not clear what, if any,
model imposed on the observed, voxel-level data would be consistent with this
assumption.

To provide an initial and rudimentary analysis of this estimator for data at
the voxel level, Achard et al. (2011) considered a model111In fact, the model
of Achard et al. (2011) was applied to the discretely observed signals, while
we generalize it to the latent fields. for the $Y_{j}$ as jointly Gaussian
random fields with
$\mathrm{Corr}(Y_{j}(v,t),Y_{j^{\prime}}(w,t))=\rho_{jj^{\prime}}^{*}$ for all
$(v,w)\in\mathcal{R}_{j}\times\mathcal{R}_{j^{\prime}}$ and
$t\in[0,\mathcal{T}].$ That is, the correlation between latent signals
observed simultaneously at any pair of spatial indices $(v,w)$ is constant for
each pair of regions. This is a natural way to generalize the role of the
region-level correlation parameter to be homogeneous across different voxels,
as it treats these as exchangeable, though not independent. Note that this
correlation $\rho_{jj^{\prime}}^{*}$ may differ from the aforementioned
parameter $\rho_{jj^{\prime}},$ as will now be elucidated.

Let us further suppose that $\xi^{2}_{j}=\mathrm{Var}(Y_{j}(v,t))$ for all
$v\in\mathcal{R}_{j}$ and $t\in[0,\mathcal{T}].$ Then the analysis of Achard
et al. (2011) shows that

$\operatorname{Corr}\left(\bar{X}_{jm},\bar{X}_{j^{\prime}m}\right)=\frac{\rho_{jj^{\prime}}^{*}}{\sqrt{(\alpha_{j}+\beta_{j})(\alpha_{j^{\prime}}+\beta_{j^{\prime}})}},$
(2.4)

where
$\alpha_{j}=L_{j}^{-2}\sum_{l,l^{\prime}=1}^{L_{j}}\operatorname{Corr}\left(Y_{j}(v_{jl},t_{m}),Y_{j}(v_{jl^{\prime}},t_{m})\right)$
is the average of all pairwise intra-regional correlations for region $j$ and
$\beta_{j}=\sigma^{2}/(L_{j}\xi_{j}^{2})$ can be considered as the noise-to-
signal ratio. Hence, under mild assumptions on the temporal dependence,
$\hat{\rho}^{\textrm{CA}}$ will converge almost surely to the right-hand side
of (2.4) as the duration $\mathcal{T}$ of the scan expands and the number $M$
of temporal observation points diverges. At least two important facts can be
gleaned from (2.4): first, that $\hat{\rho}^{\textrm{CA}}$ can be extremely
biased as an estimator of $\rho_{jj^{\prime}}^{*}$, and second, that the
correlation value targeted by $\hat{\rho}^{\textrm{CA}}$ is heavily influenced
by strength of within-region correlations and the noise level. For instance,
in an extreme scenario in which there is no noise and
$\alpha_{j}\alpha_{j^{\prime}}$ is small, one is forced to conclude that
$\rho_{jj^{\prime}}^{*}$ cannot be large.

Achard et al. (2011) proposed various corrections to
$\hat{\rho}^{\textrm{CA}}$ to appropriately target the parameter
$\rho_{jj^{\prime}}^{*}$. However, this does not eliminate the constraints
placed on $\rho_{jj^{\prime}}^{*}$ by the intra-regional correlation parameter
$\alpha_{j}$ and noise variance. These effects are somewhat unnatural, as the
connectivity parameter should intuitively be intrinsic to the brain itself and
not dependent on the sampling scheme (via the spatially averaged intra-
regional correlations) or noisy perturbations imposed by the machine.
Moreover, this basic model has other inadequacies; for example, it is not
clear what intra-regional correlation values are valid in the sense that they
lead to a spatiotemporal covariance field that is positive semidefinite. In
the next section, these technical difficulties will be resolved by formulating
a linear mixed model for the voxel-wise fMRI signals that encapsulates the
essence of this basic model while clearly distinguishing the various sources
of inter- and intra-regional dependencies.

###  2.2 A BOLD Random Effects Model

From a statistical modeling perspective, it is natural to view the collection
of signals in each region as a cluster that is, in some sense, homogeneous
compared to signals from other regions. At the same time, there must be a
common source of variability across different clusters in order to allow for
correlation between brain regions. Thus, to account for different sources of
variability that induce the characteristics observed in fMRI data of inter-
regional dependence, intra-regional dependence, and noise, we propose to
formulate a linear mixed effect model. Although a long-standing staple in
statistic modeling (not to mention specific fields such as psychometrics,
econometrics, and many others), to the knowledge of the authors, this basic
approach has not been used to develop a systematic approach to individual
brain resting-state functional connectivity estimation from voxel-level data.

Begin by assuming a basic linear model for the latent BOLD signals as

$Y_{j}(v,t)=\mu_{j}+\eta_{j}(t)+\gamma_{j}(v,t),\,v\in\mathcal{R}_{j},\,j=1,\ldots,J,\,t\in[0,\mathcal{T}].$
(2.5)

Here, $\mu_{j}$ is the deterministic mean signal, which will be assumed to be
constant in time throughout, but may vary across regions. The process
$\eta_{j}$ is a common zero-mean signal that is shared by all voxels in a
given region, and is treated as a random effect. This term represents the
desired homogeneity of signals within each region, consistent with the
characterization of connectivity found elsewhere (e.g., Van Den Heuvel and
Pol, 2010). The zero-mean spatiotemporal process $\gamma_{j}$ represents the
voxel-specific variations from the regional signals that are correlated across
both spatial and temporal domains. However, the spatial dependence of the
fields $\gamma_{j}$ is assumed to persist over a short distance; specifically,
the fields for different regions are assumed to be uncorrelated. Thus, the
dependence between random effects $\eta_{j}$ of different regions
characterizes inter-regional connectivity while also inducing intra-regional
correlation, and the random perturbations from $\gamma_{j}$ can be considered
as voxel-level spatiotemporal fluctuations resulting in intra-regional
correlations that are idiosyncratic to each brain region.

The assumed dependence structures for the latent random fields $Y_{j}$ can be
summarized as follows. First, we assume the collection of temporal regional
signals $\\{\eta_{j}:j=1,\ldots,J\\}$ is uncorrelated with the voxel-specific
spatiotemporal fluctuations $\\{\gamma_{j}:j=1,\ldots,J\\}$. Furthermore, the
random fields $\gamma_{j}$ and $\gamma_{j^{\prime}}$ are uncorrelated whenever
$j\neq j^{\prime}$. The spatiotemporal covariance structures for $\eta_{j}$
(across regions and time) and $\gamma_{j}$ (within region and across time) can
be generically modeled as follows. Let
$\bm{R}=\\{\rho_{jj^{\prime}}\\}_{j,j^{\prime}=1}^{J}$ be a $J\times J$
positive definite correlation matrix, $A_{\eta}$ and $B_{j}$ be positive
definite covariance kernels on $[0,\mathcal{T}]^{2}$, and $C_{j}$ be a
positive definite spatial correlation kernel on $\mathcal{R}_{j}$. Then the
assumed covariance model for (2.5) is

$\operatorname{Cov}\left(\eta_{j}(t),\eta_{j^{\prime}}(t^{\prime})\right)=\rho_{j{j^{\prime}}}A_{\eta}(t,t^{\prime}),\quad\operatorname{Cov}\left(\gamma_{j}(v,t),\gamma_{j}(v^{\prime},t^{\prime})\right)=B_{j}(t,t^{\prime})C_{j}(v,v^{\prime}).$
(2.6)

Due to the problem motivating the model in (2.5), the primary parameters of
interest are contained in the correlation matrix $\bm{R}$, while all other
components are viewed as auxiliary nuisance parameters.

We pause to remark that (2.5) can be viewed as a model for spatial functional
data, for which significant work exists; see, for example, Delicado et al.
(2010); Hörmann and Kokoszka (2011); Gromenko et al. (2012); French and
Kokoszka (2021); Zhang and Li (2022); Yarger, Stoev and Hsing (2022). However,
the usual treatment of the temporal aspect of spatially-indexed functional
data does not align particularly well with fMRI data. Spatial functional data
methods typically assume a fixed compact domain $[0,\mathcal{T}]$ with
increasingly dense observations $t_{m}$ while, as a practical matter, one is
forced to expand the domain $[0,\mathcal{T}]$ in order to obtain further
temporal measurements of the fMRI signal. In addition, the spatial observation
grid is essentially fixed for an fMRI scan, so that even in-fill spatial
asymptotics are not appropriate. Spatial functional data models are also
typically nonparametric to allow maximum flexibility, while the goal of
functional connectivity estimation is to infer the inter-regional correlation
parameters, so that nonparametric modeling may be unnecessary and even
inefficient. Finally, tasks such as smoothing and kriging that are staples in
functional methods do not directly address the matter at hand.

Combining (2.1) and (2.5), we obtain

$X_{jlm}=\mu_{j}+\eta_{jm}+\gamma_{jlm}+\epsilon_{jlm},$ (2.7)

where $\eta_{jm}=\eta_{j}(t_{m})$ and $\gamma_{jlm}=\gamma_{j}(v_{jl},t_{m}),$
$1\leq l\leq L_{j},$ $1\leq j\leq J,$ $1\leq m\leq M$ are the random effects,
$\mu_{j}$ is the fixed effect, and $\epsilon_{jlm}$ are errors independent
across all indices. Moreover,
$\operatorname{Cov}\left(\gamma_{jlm},\gamma_{j^{\prime}l^{\prime}m^{\prime}}\right)=0$
whenever $j\neq j^{\prime}$ and
$\operatorname{Cov}\left(\eta_{jm},\gamma_{j^{\prime}l^{\prime}m^{\prime}}\right)=0$
for all possible indices. These various effects can be practically interpreted
as follows. $\mu_{j}$ is an overall signal level related to the
neurophysiological behavior during rest, $\eta_{jm}$ are random regional
effects that characterize the long-range, inter-regional connectivity that is
the target of brain functional connectivity studies, $\gamma_{jlm}$ represent
local variations in blood flow that induce further correlated behavior between
signals that are measured closely in time and space, and $\epsilon_{jlm}$
represent additional noise inherent to the scanner or other external sources.

###  2.3 Model Estimation

The proposed model in (2.7) coherently incorporates the inter-regional
functional connectivity parameters $\rho_{jj^{\prime}}$ with the voxel-level
resolution of the observed data, all while allowing for heterogeneous intra-
regional correlation behavior as well as noise, thus providing a promising
foundation for voxel-level assessment of regional functional connectivity. One
can approach estimation in a variety of ways, depending on further assumptions
placed on the various effects. In this paper, we will illustrate the utility
of (2.1) in the basic Gaussian setting, i.e., by assuming that the random
components have a joint Gaussian distribution, and by placing stationary
structures on the temporal and spatial covariance kernels. Specifically, for
the functions in (2.6), we assume that $A_{\eta}$ and $B_{j}$ are temporally
stationary covariance kernels, while $C_{j}$ is a spatially stationary
correlation kernel. That is, with $\delta_{ab}=\mathbf{1}(a=b)$ and
$\left\lVert\cdot\right\rVert_{2}$ denoting the Kronecker delta and Euclidean
norm, respectively,

$\begin{split}A_{\eta}(t,t^{\prime})&=k_{\eta}G(|t-t^{\prime}|;\tau_{\eta})+\sigma_{\eta}^{2}\delta_{tt^{\prime}},\\\
B_{j}(t,t^{\prime})&=k_{\gamma_{j}}H_{j}(|t-t^{\prime}|;\tau_{\gamma_{j}}),\\\
C_{j}(v,w)&=K_{j}(\left\lVert
v-w\right\rVert_{2};\phi_{\gamma_{j}}),\end{split}$ (2.8)

for valid stationary correlation kernels $G,H_{j},$ and $K_{j}$ defined on the
non-negative real line, and variances
$k_{\eta},k_{\gamma_{j}},\sigma^{2}_{\eta}>0,$ where the kernels are assumed
to be known up to the rate parameters
$\tau_{\eta},\tau_{\gamma_{j}},\phi_{\gamma_{j}}>0,$ $j=1,\ldots,J$. We model
the covariance structure of $\eta_{j}$ with a nugget effect
$\sigma^{2}_{\eta}$ to allow for short-scale variability that is under the
measurement unit of $t$. Numerically, this nugget effect also serves to
stabilize the estimation as the gradient of $\rho_{jj^{\prime}}$ is tied to
that of $A_{\eta}(t,t^{\prime})$, and small values of $k_{\eta}$ could result
in large variability in the estimated $\rho_{jj^{\prime}}$ if such a nugget
effect is not present.

Suppose now that $\\{(\eta_{1}(t),\ldots,\eta_{J}(t)):t\in[0,\mathcal{T}]\\}$
is a multivariate Gaussian process, $\gamma_{1},\ldots,\gamma_{J}$ are
independent Gaussian spatiotemporal random fields, and
$\epsilon_{ljm}\overset{\mathrm{iid}}{\sim}\mathbf{N}(0,\sigma^{2})$. Then
(2.8) implies that (2.7) is a Gaussian linear mixed effect model. However, as
Gaussian likelihood usually requires, with $n=M\sum_{j=1}^{J}L_{j}$ being the
total number of data points , $O(n^{3})$ computing time and $O(n^{2})$ memory,
the large size of the data set usually makes a straight shot at likelihood
evaluation prohibitive. For context, a typical human fMRI scan consists of a
voxel count on the order of 10,000 with $M$ on the order of 100, while the
corresponding parameters for a rate are around 3000 voxels and a few hundred
timepoints. Therefore, we will endeavor to break down the computation into
smaller, simpler pieces. In doing so, we will also draw a connection to the
usual estimator $\hat{\rho}^{\textrm{CA}}$ in the context of this model, and
propose two new alternative functional connectivity estimators.

As a full maximum likelihood approach is prohibitive, we will restrict our
attention to individual pairs of regions, motivated by the position of the
correlation $\rho_{jj^{\prime}}$ as the primary parameters of interest.
Specifically, we propose a two-stage estimation approach. In the first step,
one isolates data for each region in order to estimate the covariance
parameters associated with the intra-regional spatiotemporal noise
$\gamma_{j}$. In the second step, one isolates each pair of regions in order
to estimate the remaining parameters, which include the inter-regional
correlations.

####  2.3.1 Stage $1$: Estimating Region-Specific Parameters

In the first step, one considers all data for each region $\mathcal{R}_{j}$
individually. A benefit of this approach is that all regional parameters

$\bm{\theta}_{j}=\\{\phi_{\gamma_{j}},k_{\gamma_{j}},\tau_{\gamma_{j}}\\}$
(2.9)

can be estimated in parallel. At this stage, all signals within a same region
have a common component signal that is a single realization of $\eta_{j}$,
which will be treated as a fixed effect. In other words, we will condition the
observed data on the unobserved latent process $\eta_{j}$, and maximize the
conditional likelihood. The _intra-regional_ model of the observed BOLD
signals, conditional on $\bm{\eta}_{j}=(\eta_{j1},\ldots,\eta_{jM})$, is

$X_{jlm}\mid\bm{\eta}_{j}=\nu_{jm}+\gamma_{jlm}+\epsilon_{jlm},\quad\nu_{jm}=\mu_{j}+\eta_{jm}.$
(2.10)

Note that, since the spatiotemporal noise $\gamma_{j}$ and measurement errors
are assumed to be independent of $\eta_{j}$, their variance components are not
affected by conditioning. Thus, the sole effect of conditioning on $\eta_{j}$
is to treat it as a fixed effect in the first step of estimation. There are
two reasons for doing this, one practical and the other computational. As a
practical matter, by estimating the regional effects $\nu_{jm}$ as a single
fixed effect, one mimics the first step in the Correlation of Averages
estimator outlined in Section 2.1. Indeed, if one assumes independence in both
time and space, $\bar{X}_{jm}$ in (2.2) would be the maximum likelihood
estimator of $\nu_{jm}$. However, under dependence, the Correlation of
Averages approach will be less efficient, though computationally simpler, as
an estimator of $\nu_{jm},$ since it does not account for any spatiotemporal
dependencies. The conditional treatment of $\eta_{j}$ also has computational
advantages associated with the use of a Gaussian likelihood, as it renders a
Kronecker covariance structure that enables faster computation of inverse
covariance matrices that would not be possible to exploit if treating
$\eta_{j}$ as a random effect in _Stage $1$_.

Because the number of fixed effects now scales with the number time points
$M$, we reduce the dimension of the problem by using a B-spline basis
expansion. In terms of the neuroscience problem at hand, this implies that the
underlying latent signals that embody the collaborative operation of distinct
brain regions possesses some degree of smoothness, which is an intuitive
assumption. Specifically, suppose $\\{\psi_{k}\\}_{k=1}^{K}$ is a cubic
B-spline basis on a given set of interior knots $\\{s_{u}\\}_{u=1}^{K-4}$.
Then we approximate the fixed effects by
$\nu_{jm}\approx\sum_{k=1}^{K}v_{jk}\psi_{k}(t_{m})$. Let
$\bm{v}_{j}=[v_{j1},\dots,v_{jK}]^{\intercal}$. An initial estimate of the
fixed effect represented by $\bm{v}_{j}$ can be obtained by fitting a
regression spline on the sample pairs $\\{t_{m},X_{jlm}\\}$,
$l=1,\dots,L_{j}$, $m=1,\dots,M$, which yields $\hat{\bm{v}}_{j}.$ This
estimate can either be fixed throughout the minimization of the negative
likelihood function, or simply be used as a starting point that is updated
along with the variance components.

For $j=1,\dots,J$, denote the observed signals for region $\mathcal{R}_{j}$ by

$\bm{X}_{j}=[X_{j11},X_{j12},\dots,X_{j1M},X_{j21},\dots,X_{jL_{j}M}]^{\intercal},$
(2.11)

which is a column vector of stacked BOLD signals. Similarly, set

$\begin{split}\bm{\gamma}_{j}=&[\gamma_{j11},\gamma_{j12},\dots,\gamma_{j1M},\gamma_{j21},\dots,\gamma_{jL_{j}M}]^{\intercal},\\\
\bm{\epsilon}_{j}=&[\epsilon_{j11},\epsilon_{j12},\dots,\epsilon_{j1M},\epsilon_{j21},\dots,\epsilon_{jL_{j}M}]^{\intercal}.\end{split}$
(2.12)

Define the matrix

$\tilde{\bm{G}}=\begin{bmatrix}\psi_{1}(t_{1})&\psi_{2}(t_{1})&\dots&\psi_{K}(t_{1})\\\
\psi_{1}(t_{2})&\psi_{2}(t_{2})&\dots&\psi_{K}(t_{2})\\\
\vdots&\vdots&\ddots&\vdots\\\
\psi_{1}(t_{M})&\psi_{2}(t_{M})&\dots&\psi_{K}(t_{M})\\\ \end{bmatrix},$
(2.13)

and $\bm{G}=\bm{1}_{L_{j}}\otimes\tilde{\bm{G}}$, where $\otimes$ denotes the
Kronecker product and $\bm{1}_{L_{j}}$ is the column vector of $1$’s with
length $L_{j}$. Then the matrix form of (2.10) is

$\bm{X}_{j}\mid\bm{\eta}_{j}=\bm{G}_{j}\bm{v}_{j}+\bm{\gamma}_{j}+\bm{\epsilon}_{j}.$
(2.14)

Looking ahead, in order to estimate the regional parameters (2.9) in this
model, we construct and solve for restricted maximum likelihood (ReML)
problem, or equivalently, minimize the negative value of the restricted log
likelihood function of $\bm{X}_{j}|\bm{\eta}_{j}$. The reason that we prefer
ReML over the usual maximum likelihood (ML) is because ReML yields
asymptotically unbiased estimators for the covariance components by fitting
likelihood based on a set of contrasts that eliminates the fixed effects
(Harville, 1974; Jennrich and Schluchter, 1986; Lindstrom and Bates, 1988;
Pinheiro and Bates, 2006). We note that it is feasible and reasonable to
recover the fixed effects from ReML estimates. However, the caveat is that the
ReML estimated fixed effects are not invariant under change of basis of data,
which is a desired property possessed by the usual ML estimators (see Pinheiro
and Bates, 2006, page $76$). It is clear that the ReML function includes the
variance of measurement errors $\sigma^{2}$, which enters the covariance
matrix diagonally in an additive manner. If one factors out $\sigma^{2}$ from
the entire covariance structure of $\bm{X}_{j}|\bm{\eta}_{j}$,
$\hat{\sigma}^{2}$ can be analytically represented by the other parts of the
likelihood function, which makes the computation process faster and more
accurate; see (2.18) and (2.26).

Provided that one uses the analytical solution for $\hat{\sigma}^{2}$, the
model parameter estimates are invariant to the choice between estimating
$k_{\gamma_{j}}$ directly, or estimating the variance ratio
$\tilde{k}_{\gamma_{j}}=k_{\gamma_{j}}/\sigma^{2}$, as a whole. We adopt the
latter approach for improved computation efficiency. Indeed, one can observe
that the estimation of the primary parameters $\rho_{jj^{\prime}}$ only
requires estimating the variance ratios of the model components. Therefore,
letting
$\\{\tilde{\bm{B}}_{j}\\}_{mm^{\prime}}=\tilde{B}_{j}(t_{m},t_{m^{\prime}})=\tilde{k}_{\gamma_{j}}H_{j}(t_{m},t_{m^{\prime}};\tau_{\gamma_{j}})$,
$\\{\bm{C}_{j}\\}_{ll^{\prime}}=K_{j}(\left\lVert
v_{jl}-v_{jl^{\prime}}\right\rVert_{2};\phi_{\gamma_{j}}),$ and $\bm{I}_{A}$
be the $A\times A$ identity matrix, the marginal model of (2.10) is

$\bm{X}_{j}\sim\mathbf{N}(\bm{G}_{j}\bm{v}_{j},\sigma^{2}\bm{V}_{j}(\bm{\theta}_{j})),$
(2.15)

where

$\bm{V}_{j}(\bm{\theta}_{j})=\bm{C}_{j}\otimes\tilde{\bm{B}}_{j}+\bm{I}_{(ML_{j})}.$
(2.16)

Suppressing dependence of $V_{j}$ on $\bm{\theta}_{j},$ the corresponding
negative ReML function is

$\displaystyle\mathcal{L}_{j}^{\textrm{ReML}}(\bm{v}_{j},\bm{\theta}_{j}\mid\bm{X}_{j},\bm{\eta}_{j})=$
$\displaystyle\frac{1}{2}\log\left\lvert\bm{V}_{j}\right\rvert$ (2.17)
$\displaystyle+\frac{1}{2}\log\left\lvert\bm{G}_{j}^{\intercal}\bm{V}_{j}^{-1}\bm{G}_{j}\right\rvert$
$\displaystyle+\frac{1}{2}(ML_{j}-K)\log\left(\bm{X}_{j}-\bm{G}_{j}\bm{v}_{j}\right)^{\intercal}\bm{V}_{j}^{-1}\left(\bm{X}_{j}-\bm{G}_{j}\bm{v}_{j}\right),$

where $\sigma^{2}$ has been profiled using

$\sigma^{2}=\frac{1}{L_{j}M-K}\left(\bm{X}_{j}-\bm{G}_{j}\bm{v}_{j}\right)^{\intercal}\hat{\bm{V}}_{j}^{-1}\left(\bm{X}_{j}-\bm{G}_{j}\bm{v}_{j}\right).$
(2.18)

Estimating the regional covariance components $\bm{\theta}_{j}$ by minimizing
the negative ReML function in (2.17) is referred to as _Stage $1$_ estimation
in our two-stage approach. In particular, $\hat{\bm{\theta}}_{j}$ will be
plugged into the full ReML function in _Stage $2$_ (see Section 2.3.2, where
it can be further updated or held fixed to reduce computation cost. The
estimated parameters for fixed effects $\hat{\bm{v}}_{j}$ obtained by B-Spline
regression are used as the initial values of the optimization. In the case
where one updates these estimates, the analytical generalized least square
(GLS) estimator should be used rather than running the optimizer on
$\bm{v}_{j}$ directly, that is,
$\hat{\bm{v}}_{j}=(\bm{G}_{j}^{\intercal}\hat{\bm{V}}^{-1}_{j}\bm{G}_{j})^{-1}\bm{G}_{j}^{\intercal}\hat{\bm{V}}_{j}^{-1}\bm{X}_{j}.$
The Kronecker structure of the covariance matrix in (2.16) allows us to use
scalable inference approaches as shown in Saatchi (2011). See Section A.1 of
Appendix A for details on this scalable inference approach, including
Algorithm 1 that details the optimization steps.

As mentioned previously, a modified version of the correlation of averages
estimator might also be considered at the end of _Stage $1$_. Recall that one
obtains the fixed-effects estimates in _Stage $1$_ as

$\hat{\nu}_{jm}=\sum_{k=1}^{K}\hat{v}_{jk}\psi_{k}(t_{m}),\,\,j=1,\dots,L_{j},m=1,\dots,M.$
(2.19)

Then, letting $\check{\mu}_{j}$ and $\check{\zeta}^{2}_{j}$ be the empirical
mean and variance of $\hat{\nu}_{jm}$ across $m$ (time), this modified
estimator of $\rho_{jj^{\prime}}$ is

$\tilde{\rho}^{\textrm{CA}}_{jj^{\prime}}=\frac{M^{-1}\sum_{m=1}^{M}(\hat{\nu}_{jm}-\check{\mu}_{j})(\hat{\nu}_{j^{\prime}m}-\check{\mu}_{j^{\prime}})}{\check{\zeta}_{j}\check{\zeta}_{j^{\prime}}}.$
(2.20)

This estimator will be similar to the correlation of averages when the number
of knots is large (close to $M$) and the estimate of $\bm{v}_{j}$ is computed
only using the initial spline computation, without updating it in the
likelihood maximization step.

####  2.3.2 Stage $2$: Estimating Global and Inter-Regional Parameters

Having obtained estimates for the regional variance components in _Stage $1$_,
we advance to pairwise regional data. Without loss of generality, consider
regions $\mathcal{R}_{1}$, $\mathcal{R}_{2}$. Let
$\bm{X}=[\bm{X}^{\intercal}_{1},\bm{X}^{\intercal}_{2}]^{\intercal}$,
$\bm{\mu}=[\mu_{1},\mu_{2}]^{\intercal}$,
$\bm{\eta}=[\bm{\eta}^{\intercal}_{1},\bm{\eta}^{\intercal}_{2}]^{\intercal}$,
$\bm{\gamma}^{\intercal}=[\bm{\gamma}^{\intercal}_{1},\bm{\gamma}_{2}]^{\intercal}$,
$\bm{\alpha}=[\bm{\eta}^{\intercal},\bm{\gamma}^{\intercal}]^{\intercal}$, and
$\bm{\epsilon}=[\bm{\epsilon}^{\intercal}_{1},\bm{\epsilon}^{\intercal}_{2}]^{\intercal}$.
Furthermore, define the design matrices

$\bm{Z}=\begin{bmatrix}\bm{1}_{ML_{1}}&\bm{0}_{ML_{1}}\\\
\bm{0}_{ML_{2}}&\bm{1}_{ML_{2}}\end{bmatrix},$

and $\bm{U}=[\bm{Z},\,\bm{I}_{M(L_{1}+L_{2})}]$, where $\bm{0}_{L}$ is the
column vector of $0$’s of length $L$. Then the pairwise _inter-regional_ model
has the matrix form

$\bm{X}=\bm{Z}\bm{\mu}+\bm{U}\bm{\alpha}+\bm{\epsilon}.$ (2.21)

Following the discussion under (2.10), define
$\tilde{k}_{\eta}=k_{\eta}/\sigma^{2}$,
$\tilde{\sigma}^{2}_{\eta}=\sigma^{2}_{\eta}/\sigma^{2}$, and
$\\{\tilde{\bm{A}}\\}_{t_{m}t_{m^{\prime}}}=\tilde{A}_{\eta}(t_{m},t_{m^{\prime}})=\tilde{k}_{\eta}G(|t_{m}-t_{m^{\prime}}|;\tau_{\eta})+\tilde{\sigma}^{2}\delta_{mm^{\prime}}$.
Let $\bm{\theta}$ be the vector of parameters that characterizes the
covariance structure of $\bm{X}$, which is a superset of the regional
parameter vectors $\bm{\theta}_{1},\bm{\theta}_{2}$, given by

$\bm{\theta}=[\tau_{\eta},k_{\eta},\phi_{\gamma_{1}},\phi_{\gamma_{2}},\tilde{\tau}_{\gamma_{1}},\tilde{k}_{\gamma_{1}},\tilde{\tau}_{\gamma_{2}},\tilde{k}_{\gamma_{2}},\rho_{12},\tilde{\sigma}^{2}_{\eta}]^{\intercal}.$
(2.22)

Then, the corresponding marginal model is

$\bm{X}\sim\mathbf{N}(\bm{Z}\bm{\mu},\sigma^{2}\bm{V}(\bm{\theta})),$ (2.23)

where

$\bm{V}(\bm{\theta})=\begin{bmatrix}\bm{C}_{1}\otimes\tilde{\bm{B}}_{1}+\bm{J}_{L_{1},L_{1}}\otimes\tilde{\bm{A}}&\rho_{12}\bm{J}_{L_{1},L_{2}}\otimes\tilde{\bm{A}}\\\
\rho_{12}\bm{J}_{L_{2},L_{1}}\otimes\tilde{\bm{A}}&\bm{C}_{2}\otimes\tilde{\bm{B}}_{2}+\bm{J}_{L_{2},L_{2}}\otimes\tilde{\bm{A}}\end{bmatrix}+\bm{I}_{M(L_{1}+L_{2})}$
(2.24)

and $\bm{J}_{L,L^{\prime}}$ is the $L\times L^{\prime}$ matrix of $1$’s. To
simplify the notation, we will reduce $\bm{V}(\bm{\theta})$ to $\bm{V}$. One
can estimate the model parameters $\bm{\theta}$, $\bm{\mu}$ by minimizing the
negative value of the ReML function for inter-regional model, given by

$\displaystyle\mathcal{L}^{\textrm{ReML}}(\bm{\mu},\bm{\theta}\mid\bm{X})=$
$\displaystyle\frac{1}{2}\log\left\lvert\bm{V}\right\rvert$ (2.25)
$\displaystyle+\frac{1}{2}\log\left\lvert\bm{Z}^{\intercal}\bm{V}^{-1}\bm{Z}\right\rvert$
$\displaystyle+\frac{1}{2}(M(L_{1}+L_{2})-2)\log\left(\bm{X}-\bm{Z}\bm{\mu}\right)^{\intercal}\bm{V}^{-1}\left(\bm{X}-\bm{Z}\bm{\mu}\right).$

The estimates from _Stage $1$_ can either be fixed throughout the optimization
iterations, or be used as initial values for the full optimization procedure
and updated during the optimization. Similar to _Stage $1$_, $\hat{\bm{\mu}}$
is obtained by its generalized least squares (GLS) estimator. Again, profiling
has been applied via

$\sigma^{2}=\frac{1}{M(L_{1}+L_{2})-2}\left(\bm{X}-\bm{Z}\bm{\mu}\right)^{\intercal}\bm{V}^{-1}\left(\bm{X}-\bm{Z}\bm{\mu}\right).$
(2.26)

We refer the estimation of $\bm{\theta}$ as _Stage $2$_ of our two-stage
approach, and denote the estimated correlation parameter as
$\hat{\rho}^{\textrm{ReML}}_{12}$. While the covariance matrix in (2.24) no
longer possesses the Kronecker structure from _Stage $1$_, we exploit its
block structure by using the Schur complement (Petersen and Pedersen, 2008) to
efficiently compute inverses when evaluating the loss function in (2.25).
Details are included in Section A.2 of Appendix A, including the computational
details in Algorithm 2.

In our experiments in Sections 3 and 4, two optimizers were investigated to
carry out the two-stage estimation approach. The first one is the gradient
based limited-memory-Broyden–Fletcher–Goldfarb–Shanno (LBFGS) solver (Nocedal,
1980), which is a widely adapted optimizer that implements a quasi-Newton
method. By levering an approximate Hessian matrix, this method provides an
optimization algorithm that is efficient in both time and memory. The
potential pitfall is that the method can become unstable if elements of the
gradient become small. To mitigate this issue, we also used a trust region
based optimizer, the bound-optimization-by-quadratic-approximation (BOBYQA)
solver (Powell, 2009). Although slower than a gradient-based optimizer, this
derivative-free algorithm was found to be more stable in our experiments on
real data. As a result, we use the LBFGS optimizer in the simulation studies
for its efficiency and BOBYQA in real data application for its stability.

###  2.4 Asymptotic Inference for $\hat{\rho}^{\textrm{ReML}}$ estimators

Another advantage of the proposed model is that the influence of the
spatiotemporal dependence between the various signals can be incorporated into
the assessment of estimation uncertainty by using well-established asymptotic
properties of ReML estimators. In keeping with the nature of fMRI data, the
appropriate asymptotic regime is that of an expanding time window
$\mathcal{T}$, with the number of time measurements $M$ diverging accordingly.
While the spatial dimensions, or region sizes, remain fixed, we remark that
these still plays a crucial role since the spatial design is not a regular
lattice and the spatial correlation is only assumed to be stationary within a
region and not across the entire brain. General asymptotic properties for ReML
estimators under the spatial regression setting were investigated by Cressie
and Lahiri (1993) in which the main focus was on the parameters of the
residual covariance structure, matching the focus of our motivating problem.
Cressie and Lahiri (1996) further proposed practically verifiable sufficient
conditions for the asymptotic results of ReML estimators to hold under various
common settings.

For simplicity, as well as to be faithful to the computational approach
outlined above, we always consider the pairwise model (2.21), so that the
sample size is $n=M(L_{1}+L_{2})$ and the model parameters of interest is
$\bm{\theta}\in\mathbb{R}^{10}$, defined in (2.22), which contains the key
correlation parameter. Define

$\bm{\Pi}(\bm{\theta})=\bm{V}^{-1}(\bm{\theta})-\bm{V}^{-1}(\bm{\theta})\bm{Z}(\bm{Z}^{\intercal}\bm{V}^{-1}(\bm{\theta})\bm{Z})^{{}^{-1}}\bm{Z}^{\intercal}\bm{V}^{-1}(\bm{\theta}),$
(2.27)

which is the projection of $\bm{V}^{-1}(\bm{\theta})$ onto the orthogonal
complement column space of the design matrix $\bm{Z}$, arising from the ReML
technique of removing the fixed effects. Let
$\bm{V}_{i}(\bm{\theta})=\partial\bm{V}(\bm{\theta})/\partial\theta_{i}$,
where $\theta_{i}$ is the $i$th element of $\bm{\theta}$, $i=1,\dots,10.$
Similarly,
$\bm{V}_{ij}(\bm{\theta})=\partial^{2}\bm{V}(\bm{\theta})/\partial\theta_{i}\partial\theta_{j}$,
$i,j=1,\dots,10$. Furthermore, let $\bm{\mathcal{I}}(\bm{\theta})$ be the
$p\times p$ matrix of second-order partial derivatives of the negative ReML
function, that is,
$\\{\bm{\mathcal{I}}(\bm{\theta})\\}_{ij}=\partial^{2}\mathcal{L}^{\textrm{ReML}}(\bm{\theta})/\partial\theta_{i}\partial\theta_{j},1\leq
i,j\leq p$. It is straightforward to show that the $ij$th element of Fisher
information matrix is

$\mathbb{E}_{\bm{\theta}}(\bm{\mathcal{I}}(\bm{\theta}))_{ij}=\frac{1}{2}\operatorname*{Tr}\left\\{\bm{\Pi}(\bm{\theta})\bm{V}_{i}(\bm{\theta})\bm{\Pi}(\bm{\theta})\bm{V}_{j}(\bm{\theta})\right\\},\quad
1\leq i,j\leq p.$ (2.28)

Application of the results in Theorems 3.1 and 3.2 of Cressie and Lahiri
(1993) to our setting is then relatively straightforward. As time is the only
dimension that is involved in the asymptotic development, the assumptions used
in Cressie and Lahiri (1993) to derive these results amount to restrictions in
the temporal correlation decay represented by the matrices
$\tilde{\bm{B}}_{j}$ and $\tilde{\bm{A}}.$ When these are satisfied, one
immediately finds that

$\left\\{\mathbb{E}{\bm{\mathcal{I}}_{n}(\bm{\theta})}\right\\}^{1/2}(\hat{\bm{\theta}}^{\textrm{ReML}}-\bm{\theta})\rightsquigarrow\mathbf{N}(\bm{0},\bm{I}_{p}),$
(2.29)

with $\rightsquigarrow$ denoting convergence in distribution. Focusing only on
the primary parameter $\hat{\rho}^{\textrm{ReML}},$ a consequence of the above
is that

$\left\\{\frac{1}{2}\operatorname*{Tr}\left\\{\left(\bm{\Pi}(\bm{\theta})\begin{bmatrix}\bm{0}&\bm{J}_{L_{1},L_{2}}\\\
\bm{J}_{L_{2},L_{1}}&\bm{0}\end{bmatrix}\otimes\tilde{\bm{A}}\right)^{2}\right\\}\right\\}^{\frac{1}{2}}(\hat{\rho}^{\textrm{ReML}}_{12}-\rho_{12})\rightsquigarrow\mathbf{N}(\bm{0},1).$
(2.30)

The standardizing matrix in (2.30) can be readily estimated by plugging in the
relevant components of the computed estimate $\hat{\bm{\theta}}$ into
$\bm{\Pi}(\bm{\theta})$ and $\tilde{\bm{A}}.$

##  3 Simulation Study

To understand the performance of our model, we proceed with a simulation
study. We simultaneously simulate BOLD signals from 3 regions under different
settings by varying signal strength and intra-regional correlation. For the
temporal correlation kernels $G$ and $H_{j}$ in (2.8), we use the squared
exponential kernel (also called the Gaussian or radial basis function (RBF)
kernel)

$\kappa(u;\tau)=\exp\left(-\dfrac{\tau^{2}u^{2}}{2}\right),$ (3.1)

setting $G(u;\tau_{\eta})=\kappa(u;\tau_{\eta})$ and
$H_{j}(u;\tau_{\gamma_{j}})=\kappa(u;\tau_{\gamma_{j}})$. For the spatial
correlation kernel $K_{j}$ in (2.8), we used the Matérn-5/2 kernel

$K_{j}(d;\phi_{\gamma_{j}})=\left(1+\sqrt{5}\phi_{\gamma_{j}}d+\dfrac{5}{3}\phi_{\gamma_{j}}^{2}d^{2}\right)\exp\left(-\sqrt{5}\phi_{\gamma_{j}}d\right),$
(3.2)

that is a popular choice in spatial statistics (Stein, 1999).

One of the main goals is to study the performance of our estimator
$\hat{\rho}^{\textrm{ReML}}$ in comparison to $\hat{\rho}^{\textrm{CA}}$ under
different settings of signal strength relative to spatiotemporal noise, which
is mainly quantified by the values of $\tilde{k}_{\eta}$ relative to those of
$\tilde{k}_{\gamma}$. Therefore, we set $\tilde{k}_{\gamma}=2$,
$\tau_{\gamma}=0.5$ for all regions and fix $\sigma^{2}=1$,
$\tau_{\eta}=0.25$, and $\tilde{\sigma}^{2}_{\eta}=0.1$ for all simulations.
Additionally, we fix the mean signals to be $\mu_{1}=1,\mu_{2}=10,\mu_{3}=20$
and the pairwise inter-regional correlations are
$\rho_{12}=0.1,\rho_{13}=0.35,\rho_{23}=0.6$. For each setting, $Q=100$
independent simulations were run with $L_{j}\equiv 50$ voxels for each region
(randomly sampled from a 3D rectangular lattice with side length of 7) and
$M=60$ timepoints for all signals. In _Stage $1$_ estimation, $K=45$ cubic
B-spline basis functions were used to approximate the fixed effect using an
equally spaced set of 41 interior knots. The code for all simulations can be
found on Github (https://github.com/cbtran/qfuncMM-reproducible).

###  3.1 Comparison of estimators

(a) $\rho_{12}=0.1$

(b) $\rho_{13}=0.35$

(c) $\rho_{23}=0.6$

Figure 1: Empirical distribution of estimated inter-regional correlations: (a) Weak inter-regional correlation, (b) Medium inter-regional correlation, (c) Strong inter-regional correlation. In each subfigure, two intra-regional correlation rates are $\phi_{\gamma}=0.25$ (strong) and $\phi_{\gamma}=1$ (weak), and two levels of signal are $k_{\eta}=0.5$ (weak) and $k_{\eta}=1$ (strong). The three estimators are Pearson correlation of average (CA, $\hat{\rho}^{\textrm{CA}}$), Pearson correlation of fixed-effect estimate (FE, $\tilde{\rho}^{\mathrm{CA}}$), and the proposed ReML estimator (ReML, $\hat{\rho}^{\textrm{ReML}}$). Dashed line denotes the true value of inter-regional correlation. Table 1:  Estimates of correlation coefficients from simulated data. Two levels of intra-regional correlation, $\phi_{\gamma}=0.25$ (strong) and $\phi_{\gamma}=1$ (weak), and three levels of inter-regional correlation, $\rho=0.1,0.35,0.6$, and two levels of signal strength $k_{\eta}=0.5,1$ are used. The three estimators are Pearson correlation of average (CA, $\hat{\rho}^{\textrm{CA}}$), Pearson correlation of fixed-effect estimate (FE, $\tilde{\rho}^{\mathrm{CA}}$), and the proposed ReML estimator (ReML, $\hat{\rho}^{\textrm{ReML}}$).  Signal | Intra-correlation | $\rho_{\text{True}}$ | Type | RMSE | $\left\lvert\text{Bias}\right\rvert$ | SD
---|---|---|---|---|---|---
$k_{\eta}=0.5$ | $\phi_{\gamma}=1$ | 0.1 | ReML | 0.1584 | 0.0379 | 0.1546
FE | 0.3145 | 0.0101 | 0.3159
CA | 0.2793 | 0.0044 | 0.2806
0.35 | ReML | 0.1387 | 0.0057 | 0.1393
FE | 0.2925 | 0.0640 | 0.2869
CA | 0.2606 | 0.0576 | 0.2555
0.6 | ReML | 0.1231 | 0.0105 | 0.1232
FE | 0.2904 | 0.1587 | 0.2444
CA | 0.2645 | 0.1595 | 0.2121
$\phi_{\gamma}=0.25$ | 0.1 | ReML | 0.2255 | 0.0474 | 0.2216
FE | 0.2638 | 0.0205 | 0.2643
CA | 0.2271 | 0.0344 | 0.2256
0.35 | ReML | 0.1897 | 0.0136 | 0.1902
FE | 0.3537 | 0.2187 | 0.2793
CA | 0.3532 | 0.2321 | 0.2677
0.6 | ReML | 0.1550 | 0.0247 | 0.1538
FE | 0.4246 | 0.3547 | 0.2346
CA | 0.4701 | 0.4132 | 0.2254
$k_{\eta}=1$ | $\phi_{\gamma}=1$ | 0.1 | ReML | 0.1638 | 0.0133 | 0.1641
FE | 0.3104 | 0.0361 | 0.3098
CA | 0.2897 | 0.0263 | 0.2900
0.35 | ReML | 0.1593 | 0.0181 | 0.1591
FE | 0.2641 | 0.0032 | 0.2654
CA | 0.2445 | 0.0093 | 0.2455
0.6 | ReML | 0.1309 | 0.0147 | 0.1307
FE | 0.2765 | 0.1057 | 0.2568
CA | 0.2602 | 0.1108 | 0.2366
$\phi_{\gamma}=0.25$ | 0.1 | ReML | 0.2113 | 0.0034 | 0.2124
FE | 0.3012 | 0.0100 | 0.3026
CA | 0.2670 | 0.0306 | 0.2666
0.35 | ReML | 0.2083 | 0.0235 | 0.2080
FE | 0.2996 | 0.1393 | 0.2666
CA | 0.3153 | 0.1803 | 0.2600
0.6 | ReML | 0.1604 | 0.0178 | 0.1602
FE | 0.3917 | 0.2868 | 0.2681
CA | 0.4222 | 0.3345 | 0.2589

First, we compare estimation performance of the three estimators,
$\hat{\rho}^{\textrm{ReML}}$, $\hat{\rho}^{\textrm{CA}}$, and
$\tilde{\rho}^{\textrm{CA}}$, which are the ReML estimator of $\rho$, the
conventional Pearson correlation of averages (CA), and the Pearson correlation
of averages on B-spline fixed effect estimates from _Stage $1$_ (FE),
respectively. The calculation of the two correlation of averages type
estimators are detailed in (2.3) and (2.20). Four settings were considered,
corresponding to weak ($\phi_{\gamma}=1$) and strong ($\phi_{\gamma}=0.25$)
intra-regional correlation and weak ($k_{\eta}=0.5$) and strong ($k_{\eta}=1$)
signal.

The estimation performances of the different estimators under each setting are
shown in Figure 1 and Table 1. Overall, $\hat{\rho}^{\textrm{ReML}}$ has the
smallest standard deviation for all settings, which can be be observed
visually in Figure 1 and numerically in Table 1. In the case of weak intra-
regional correlation ($\phi_{\gamma}$ = 1), the true values of $\rho$ are
between the 1st and 3rd quartiles from all estimators. Notably,
$\hat{\rho}^{\textrm{CA}}$ and $\tilde{\rho}^{\textrm{CA}}$ are biased toward
0, their medians deviating strongly from the true value when $\rho=0.6$.
Unsurprisingly, the effect of bias seems to be less severe when signal gets
stronger ($k_{\eta}=1$), resulting in better concentrations around the true
value for both estimators. Similarly, $\hat{\rho}^{\textrm{CA}}$ and
$\tilde{\rho}^{\textrm{CA}}$ are heavily biased towards $0$ in the cases of
strong intra-regional correlation ($\phi_{\gamma}=0.25$). Especially when
$\rho=0.6$, their 3rd quartiles are smaller than the true value of inter-
regional correlation. On the other hand, $\hat{\rho}^{\textrm{ReML}}$ is still
robust with its median staying close to the true inter-regional correlations
even though its spread increases. Interestingly, $\tilde{\rho}^{\textrm{CA}}$
seems to performs better than $\hat{\rho}^{\textrm{CA}}$ for strong intra-
regional correlation and strong inter-regional correlation. This result
suggests that our fixed-effect estimator $\hat{\nu}$ from (2.19) might provide
a better estimate for $\eta_{j}$ signal than regional averages, particularly
when data are not independent and identically distributed. Overall, we can
observe that $\hat{\rho}^{\textrm{ReML}}$ is the most robust estimator across
all different simulation settings. See Appendix B for additional simulation
results with very strong signal ($k_{\eta}=1.5$).

###  3.2 Validation of large sample properties

Figure 2: Coverage of confidence intervals. From top to bottom: weak signal,
strong signal. From left to right: weak intra-regional correlation, strong
intra-regional correlation. Black solid line denotes the desired $1-\alpha$
coverage.

Next, we validate the asymptotic normality of $\hat{\rho}^{\textrm{ReML}}$
from Theorem 2.30 by constructing confidence intervals. Using Fisher
transformation, the $(1-\alpha)$ confidence interval then has the form

$\displaystyle
C^{ReML}_{1-\alpha}=\tanh\left(\operatorname*{arctanh}({\hat{\rho}^{\textrm{ReML}}})\pm\Phi^{-1}(1-\alpha/2)\times
SE(\operatorname*{arctanh}(\hat{\rho}^{\textrm{ReML}}))\right),$ (3.3)

where $\hat{\rho}^{\textrm{ReML}}$ is estimated correlation via ReML from
Section 3.1, $SE(\operatorname*{arctanh}(\hat{\rho}^{\textrm{ReML}}))$ is the
standard error of transformed parameter
$\operatorname*{arctanh}(\hat{\rho}^{\textrm{ReML}})$ computed via the
$\delta$-method, and $\Phi^{-1}(\cdot)$ is the quantile function of the
standard normal distribution.

The coverage of confidence intervals which is the percentage of intervals that
cover the true inter-regional correlation for all 3 pairs at different
significant levels $\alpha$ are shown Figure 2. Similar to estimation results,
confidence intervals get better coverage as signal strength increases and
intra-regional correlation decreases. When intra-regional correlation is
strong ($\phi_{\gamma}=0.25$), confidence intervals failed to achieve the
desired $1-\alpha$ coverage with all lines are under the reference $y=x$ line.
This result is unsurprising since in the presence of strong intra-regional
correlations among voxels, effective sample size decreases, thus more samples
are needed for the large sample properties to be satisfied.

##  4 Data Application

Having demonstrated the robust empirical performance of the proposed estimator
$\hat{\rho}^{\textrm{ReML}}$ in simulation studies, we proceed to the
motivating data application. The model was applied to one dead and one
anesthetized live rat. These data sets, among others, have been previously
described and studied in Becq et al. (2020). In particular, the live rat under
analysis in this paper is labeled as rat $20160616\\_145220$ whose BOLD
signals were collected after administration of the Etomidate anesthetic. The
dead rat analyzed is labeled as rat $0160524\\_153000$ and is used primarily
to illustrate performance of the model estimation in the null case of no
connectivity. For both rats, the BOLD signals are collected at the frequency
of $2$ Hz for the duration of $30$ minutes. Preprocessing steps include motion
removal and wavelet transformation of the original signals. As the wavelet
transformation is linear, model (2.1) applies to the raw BOLD signals as well
as the wavelet-transformed data, with the only difference being that $t_{m}$
indexes different wavelet coefficients rather than time, although the time
ordering is preserved. Wavelets provide a very adequate framework for
analysing BOLD signals because of their flexibility and ability to deal with
long-memory time series (Whitcher, Guttorp and Percival, 2000; Percival and
Walden, 2006).

Popular wavelet methods such as the discrete wavelet transform (DWT) and
maximal overlap discrete wavelet transform (MODWT) have been used for
quantifying functional connectivity (Achard et al., 2006; Deuker et al., 2009;
Vértes et al., 2012). For our analysis, we use DWT with Daubechies Least
Asymmetric filter of length 8 as this transformation reduces the number of
available wavelet coefficients. This makes DWT more computationally attractive
for the method proposed in this paper. The wavelet coefficients of level $4$
were chosen according to Becq et al. (2020) where systematic variables such as
heart rat or body temperature were found to perturb the functional
connectivity less relative to other choices. For more details on these
preprocessing steps, we refer readers to Becq et al. (2020).

We focus on $21$ regions of interest out of $51$ regions, defined by the
custom atlas used in Becq et al. (2020). Full regions names and corresponding
labels can be found in Table 2. These regions consist of 10 bi-lateral region
pairs (one each in the left and right hemispheres), and one central region
(ACC). For each pair of regions $(j,j^{\prime}),1\leq j<j^{\prime}\leq 21$, we
first obtain $\hat{\rho}^{\textrm{ReML}}_{jj^{\prime}}$, then use the
asymptotic distribution derived in Theorem 2.30 to calculate the $z$-scores of
these estimators. In particular, we parameterized the inter and intra-regional
covariance matrices with radial basis function (RBF), also known as Gaussian,
kernels and the spatial correlation matrix with the Matérn-5/2 kernel. A cubic
B-spline basis of size $K=30$ was used in _Stage $1$_ of the estimation. As
previously mentioned, we use the BOBYQA for real data application for its
superior stability. Detailed steps of the optimization are outlined in
Algorithms 1 and 2 in Appendix A, and were implemented in an R package that
can be found on Github (https://github.com/cbtran/qfuncMM). We further adopt
the Benjamini–Yekutieli (BY) procedure (Benjamini and Yekutieli, 2001) to
perform hypothesis testing for significant correlations between pairs of
regions. The BY procedure is used to control the false discovery rate (FDR)
under arbitrary dependence assumptions (Logan and Rowe, 2004), which is
suitable for this application. Furthermore, since discoveries are expected for
the live rat, controlling FDR means that false discoveries are not penalized
as harshly as the controlling the family-wise error rate, which is equivalent
to FDR only under the global null in which all inter-regional correlations are
zero, making FDR less conservative and more reasonable for this application.

Table 2: Names of the 21 regions of the rat brain chosen for study, and their corresponding labels. Region name | Region label
---|---
anterior cingulate cortex | ACC
amygdalopiriform transition area right | APir r
auditory cortex right | AU r
entorhinal area right | Ent r
parietal association cortex right | Par r
primary motor cortex right | M1 r
supplementary motor cortex right | M2 r
retrosplenial cortex right | RSC r
somatosensory 1 right, barrel field | S1BF r
amygdalopiriform transition area left | APir l
auditory cortex left | AU l
entorhinal area left | Ent l
parietal association cortex left | Par l
primary motor cortex left | M1 l
supplementary motor cortex left | M2 l
retrosplenial cortex left | RSC l
somatosensory 1 left, barrel field | S1BF l
hypothalamic region right | H r
hippocampus right | HIP r
hypothalamic region left | H l
hippocampus left | HIP l

As an initial check on the model, data from the dead rat were analyzed. The
smallest unadjusted $p$-value was 0.157, so that no significant correlations
were discovered even without the BY adjustment, using any conventional level
of significance. In comparison with Becq et al. (2020), this results lines up
with their group analysis of several dead rats, although their methodology did
lead to a few spurious edges when applied to individual rats as done here.

Furthermore, in the case of a live rat, we compare the result with the network
of the same rat in Becq et al. (2020). Here, significant edges were detected
by applying the BY procedure with level $q=0.01.$ Figure 3(a) shows the
estimated $\hat{\rho}^{\textrm{ReML}}$ for all $21$ studied regions. Without
quantifying their uncertainty, these raw estimates should not be used directly
to represent a brain network since large values might also correspond to large
variations. Therefore, the network is constructed by performing the BY
procedure and is represented in Figure 3(b) by showing the estimates
$\hat{\rho}^{\textrm{ReML}}$ for significant edges as the upper triangular
portion of the weighted adjacency matrix.

To further summarize the role of each region in the network, we also
constructed node degree (ND) and functional connectivity strength (FCS) plots
in Figure 4. Here, ND is the total number of edges for each region, while FCS
is the average correlation estimate across detected edges for the region. At
rest, networks should be roughly symmetric, so Figure 4 stacks these
quantities vertically for bilateral regions, which only excludes ACC. The
minor degree of asymmetry in both ND and FCS metrics between left and right
regions is similar to that observed in Becq et al. (2020). We found that our
result is mostly consistent with the previous study. For example, under the
same Etomidate anesthetic, regions such as ACC, M1, M2, AU were found there to
be highly connected, and the proposed model also estimated a high connectivity
level for these regions. On the other hand, for Ent and RSC, which are
commonly found to be disconnected for several live rats under this anaesthetic
(Becq et al., 2020), a low level of connectivity is also observed in our
result. Overall, the network constructed using the proposed model was more
sparse, having 30 fewer edges than the network constructed using
$\hat{\rho}^{\textrm{CA}}$ after applying the same BY procedure. The ReML
network discovered 22 edges that do not appear in the Correlation of Averages
network, 17 of which belong to one of the regions ACC, M1, M2, and AU.

(a) Raw $\hat{\rho}^{\textrm{ReML}}$ values.

(b) Thresholded network by BY procedure.

Figure 3: Brain network of anesthetized rat $20160616\\_145220$, Eto-L.

Figure 4: Node degree (a) and FCS (b) of anesthetized rat $20160616\\_145220$,
Eto-L.

##  5 Discussion

A novel spatiotemporal model for voxel-level BOLD signals has been proposed
that explicitly models the regional signals that drive both inter-regional and
intra-regional correlation, local spatiotemporal fields that allow for
idiosyncratic, heterogeneous intra-regional correlation, and measurement
error. By fully addressing the spatiotemporal dependency structure, a coherent
quantification of inter-regional brain connectivity emerges, together with a
well-established and justified estimation approach. Large sample properties
were established with mild conditions on the smoothness of covariance
structure and the decaying rate of dependency. Simulation studies showed that
the key inter-regional correlation parameters can be accurately estimated,
with stark advantages compared conventional Pearson correlation of averages in
the presence of strong spatiotemporal noise. Furthermore, the application on
rat data showed that our findings are mostly consistent with previously
established results, which indicates that our model can produce reasonable
results while being superior in certain noisy regimes. In particular, where
the network constructed by the proposed model differed from the established
correlation of averages network is precisely those cases in which our
simulation studies indicate the superior reliability of
$\hat{\rho}^{\textrm{ReML}}$ compared to $\hat{\rho}^{\textrm{CA}}.$ Further
development of the method to make it applicable to a wider array of data sets
includes computational tools for a more flexible collection of dependence
structures. While this paper has focused on the Gaussian case, alternative
estimators in non-Gaussian linear mixed models can readily be applied in the
proposed two-stage approach.

Although designed for the purpose of estimating functional connectivity, the
proposed model could be useful for quantifying dependence amongst spatially
grouped time series in other contexts. In climate studies, for example, the
data sets consist of time series measuring the temperature for a given area.
Similar models, as proposed here, were used to predict the temperature on a
large area (Hengl et al., 2012; Gräler, Pebesma and Heuvelink, 2016). Another
example can be found in disease propagation, where the data sets consist of
time series measuring the number of positive cases of COVID in a specific area
(Bartolucci and Farcomeni, 2022). The model used is a Markov Chain where for
each area the sequence of latent variables across time follows a Markov chain
with initial and transition probabilities that also depend on latent variables
in neighboring areas.

Lastly, the data illustrations have utilized a pre-specified set of regions
for functional connectivity analysis (Biswal et al., 1995). This approach
relies on a predefined brain parcellation, commonly in form of an atlas, to
determine regions of interest. As a result, different subjects are also
assumed to have the same functional network nodes. Recently, data-driven
methods for functional connectivity analysis have been proposed (Van Den
Heuvel and Pol, 2010; Lv et al., 2018), thus allowing nodes in functional
networks to be subject-specific (Cui et al., 2020; Michon et al., 2022). It
will be interesting to investigate if our proposed methods for quantifying
functional connectivity can be combined with subject-specific region discovery
in order to simultaneously study variability in the spatial distribution of
functional connectivity nodes as well as connections between them.

[Acknowledgments]

The authors acknowledge the following computational facilities for providing
computational resources and technical support that have contributed to the
results reported within this publication:

  1. 1.

The Office of Research Computing at Brigham Young University. URL:
https://rc.byu.edu.

  2. 2.

The Center for Scientific Computing (CSC), which administrated funds from the
National Science Foundation (CNS-1725797) for the purchase of computational
facilities. The CSC is supported by the California NanoSystems Institute and
the Materials Research Science and Engineering Center (MRSEC; NSF DMR 1720256)
at UC Santa Barbara. URL: https://csc.cnsi.ucsb.edu.

The authors also gratefully acknowledge the support of this research provided
by grants from the US National Science Foundation’s Collaborative Research in
Computational Neuroscience program (Award IIS-2135859) and the French National
Research Agency.

##  Appendix A Scalable inference approaches

Common computation techniques, such as direction inversion or Cholesky
decomposition of the covariance matrices, result in $O(N^{3})$ time
complexity, where $N=L_{j}\times M$ for a single, and
$N=(L_{j}+L_{j^{\prime}})\times M$ for a region pair. Thus, we exploit the
structures of the covariance matrices to imrpove scalability.

###  A.1 _Stage $1$_: Region-specific Model

Recall that in region-specific model, the covariance matrix has the form
$\bm{V}_{j}(\bm{\theta}_{j})=\bm{C}_{j}\otimes\bm{B}_{j}+\sigma^{2}\bm{I}$. As
shown in Saatchi (2011), this Kronecker structure leads to an efficient
algorithm for evaluating Gaussian (restricted) log-likelihood. In particular,
let $\bm{C}_{j}=\bm{Q}_{C}{\bm{\Lambda}_{C}}{\bm{Q}_{C}^{\top}}$ and
$\bm{B}_{j}=\bm{Q}_{B}{\bm{\Lambda}_{C}}{\bm{Q}_{B}^{\top}}$ be
eigendecompositions of $\bm{C}_{j}$ and $\bm{B}_{j}$, respectively. Then, for
a vector $\bm{z}\in\mathbb{R}^{N},N=L_{j}\times M$, we have

$\displaystyle\bm{V}_{j}^{-1}\bm{z}$
$\displaystyle=\left(\bm{C}_{j}\otimes\bm{B}_{j}+\sigma^{2}\bm{I}\right)^{-1}\bm{z}$
$\displaystyle=(\bm{Q}_{C}\otimes\bm{Q}_{B})\left({\bm{\Lambda}_{C}}\otimes{\bm{\Lambda}_{B}}+\sigma^{2}\bm{I}\right)^{-1}({\bm{Q}_{C}^{\top}}\otimes{\bm{Q}_{B}^{\top}})\bm{z},$
(A.1)

where
$\left({\bm{\Lambda}_{C}}\otimes{\bm{\Lambda}_{B}}+\sigma^{2}\bm{I}\right)$ is
a diagonal matrix, and $({\bm{Q}_{C}^{\top}}\otimes{\bm{Q}_{B}^{\top}})\bm{z}$
can be efficiently computed using

$({\bm{Q}_{C}^{\top}}\otimes{\bm{Q}_{B}^{\top}})\bm{z}=\text{vectorize}({\bm{Q}_{B}^{\top}}{\bm{Z}}\bm{Q}_{C}),$

where $\bm{z}=\text{vectorize}(\bm{Z})$ for $\bm{Z}\in\mathbb{R}^{M\times
L_{j}}$. In addition,

$\log\det({{\bm{V}_{j}}})=\sum_{l=1}^{L_{j}}\sum_{m=1}^{M}\left[({\bm{\Lambda}_{C}})_{l}({\bm{\Lambda}_{B}})_{m}+\sigma^{2}\right].$
(A.2)

Since the most computationally intensive step is computing eigendecompositions
of $\bm{C}_{j}$ and $\bm{B}_{j}$, time complexity of this inference approach
is $O(L_{j}^{3}+M^{3})$ for ML and $O(L_{j}^{3}+M^{3}+K^{3})$ for ReML where
$K$ is the B-spline degrees of freedom instead of $O(L_{j}^{3}M^{3})$ when
using direct inversion or Cholesky decomposition. We then estimate region-
specific parameters using Algorithm 1.

Algorithm 1 Region-specific model

Input: Data matrix $\bm{X}_{j}$, design matrix $\bm{G}_{j}$, initial values
$\bm{\theta}_{j}^{(0)}$ and $\bm{v}_{j}^{(0)}$, number of iterations $T$

$\bm{r}_{j}^{(0)}\leftarrow\bm{X}_{j}-\bm{G}_{j}\bm{v}_{j}^{(0)}$

for $t\in 1:T$ do

Compute $\bm{C}_{j},\bm{B}_{j}$

Compute eigendecompositions
$\bm{C}_{j}=\bm{Q}_{C}{\bm{\Lambda}_{C}}{\bm{Q}_{C}^{\top}}$ and
$\bm{B}_{j}=\bm{Q}_{B}{\bm{\Lambda}_{C}}{\bm{Q}_{B}^{\top}}$

Evaluate (restricted) log-likelihood

if Convergence then

Break

else

Update $\bm{\theta}_{j}^{(t)}$ and $\bm{v}_{j}^{(t)}$

end if

end for

Output: $\hat{\bm{\theta}}_{j}$ and $\hat{\bm{v}}_{j}$

###  A.2 _Stage $2$_: Inter-regional model

For inter-regional model, we use the Schur complement (Petersen and Pedersen,
2008) to efficiently evaluate the loss function. For simplicity, let

$\displaystyle\bm{V}(\bm{\theta})$
$\displaystyle=\begin{bmatrix}\bm{V}_{11}&\bm{V}_{12}\\\
\bm{V}_{12}^{\top}&\bm{V}_{22}\end{bmatrix},$ (A.3) $\displaystyle\bm{W}$
$\displaystyle=\bm{V}_{11}-\bm{V}_{12}\bm{V}_{22}^{-1}\bm{V}_{12}^{\top}.$
(A.4)

Then we have

$\displaystyle\log\det(\bm{V})$
$\displaystyle=\log\det(\bm{W})+\log\det(\bm{V}_{22}).$ (A.5)

For a vector $\bm{z}\in\mathbb{R}^{N},N=(L_{j}+L_{k})\times M$, we want to
compute $\bm{V}^{-1}\bm{z}$, which is equivalent to solving linear systems

$\begin{bmatrix}\bm{V}_{11}&\bm{V}_{12}\\\
\bm{V}_{12}^{\top}&\bm{V}_{22}\end{bmatrix}\begin{bmatrix}\bm{z}_{1}\\\
\bm{z}_{2}\end{bmatrix}=\begin{bmatrix}\bm{b}_{1}\\\ \bm{b}_{2}\end{bmatrix},$
(A.6)

where $\bm{z}=\begin{bmatrix}\bm{z}_{1}\\\ \bm{z}_{2}\end{bmatrix}$. We first
solve for $\bm{b}_{1}$ using

$\bm{b}_{1}=\bm{W}^{-1}\left(\bm{z}_{1}-\bm{V}_{12}\bm{V}_{22}^{-1}\bm{z}_{2}\right).$
(A.7)

Then we use $\bm{b}_{1}$ to solve for $\bm{b}_{2}$ in (A.6).
Note that these computations of log-determinant and matrix inversion can be
further improved by Cholesky decomposition of block components. For example,
let $\bm{W}^{Chol}$ be the Choleskly decomposition of $\bm{W}$ such as

$\bm{W}=\bm{W}^{Chol}(\bm{W}^{Chol})^{\top}.$ (A.8)

Then

$\displaystyle\log\det(\bm{W})$ $\displaystyle=2\times\sum_{i=1}^{L_{j}\times
M}\log(\bm{W}^{Chol}_{ii}),$ (A.9)

and for a vector $\bm{z}\in\mathbb{R}^{L_{j}\times M}$, $\bm{W}^{-1}\bm{z}$
can be computed using forward-backward substitution.
Finally, inter-regional parameters are estimated using Algorithm 2.

Algorithm 2 Inter-regional model

Input: Data matrix $\bm{X}=[\bm{X}_{j}^{\top}\bm{X}_{k}^{\top}]^{\top}$,
design matrix $\bm{Z}$, region-specific estimates
$\bm{C}_{j},\bm{B}_{j},\bm{C}_{k},\bm{B}_{k}$, initial values
$\bm{\theta}^{(0)}$ and $\bm{\mu}^{(0)}$, number of iteration $T$

$\bm{r}^{(0)}\leftarrow\bm{X}-\bm{Z}\bm{\mu}^{(0)}$

for $t\in 1:T$ do

Compute $\bm{A}$

Compute $\bm{V}_{22}$ and $\bm{W}$

Evaluate (restricted) log-likelihood

if Convergence then

Break

else

Update $\bm{\theta}^{(t)}$ and $\bm{\mu}^{(t)}$

end if

end for

Output: $\hat{\bm{\theta}}$ and $\hat{\bm{\mu}}$

###  A.3 Hardware acceleration

Despite our two-stage approach effectively reducing the dimension of parameter
space, likelihood evaluation remains expensive due to the large dimensionality
of our covariance matrix. To speed up the computation, we utilize parallel
computing to run estimations simultaneously. Since it is well-known that
Graphics Processing Units (GPU) may provide a significant speed-up when
compared with Central Processing Unit (CPU) for large-scale computations
(Krizhevsky, Sutskever and Hinton, 2012; Gardner et al., 2018), we delegate
level 3 Basic Linear Algebra Subprograms (BLAS) matrix operations to GPU using
NVIDIA’s CUDA framework to accelerate matrix calculations. These operations
are sufficiently complex (e.g., matrix-matrix product, solving triangular
matrix equations) to ensure the speed gain of using GPU outshines the overhead
of transferring data between CPU and GPU. To demonstrate this point, we
benchmarked matrix inversion and log determinant calculations, which are
commonly encountered in log likelihood evaluation, using CPU and GPU for
$n\times n$ matrices, $n=1000,2000,\dots,8000$. Figure 5 shows a decisive
advantage of GPU calculation when matrix size grows over a certain level
($4000\times 4000$ in our experiment). In our data application, the covariance
matrix could easily go beyond $10000$ rows (columns) so using GPU acceleration
could greatly improve computation efficiency.

Figure 5: Comparison between time complexity of matrix operations, inversion
and log determinant, using CPU and GPU.

##  Appendix B Additional Simulation Results

This section reports additional results under the simulation design discussed
in Section 3. While the results reported in the text correspond to weak
($k_{\eta}=0.5$) and strong ($k_{\eta}=1)$ signal strengths, those reported in
Table 3 and Figure 6 below present results for the case of a very strong
signal ($k_{\eta}=1.5$).

Table 3:  Estimate of correlation coefficients from simulated data. Two levels of intra-regional correlation, $\phi_{\gamma}=0.25$ (strong) and $\phi_{\gamma}=1$ (weak), and three levels of inter-regional correlation, $\rho=0.1,0.35,0.6$ are used along with a very strong signal $k_{\eta}=1.5$. The three estimators are Pearson correlation of average (CA), Pearson correlation of fixed-effect estimate (FE), and the proposed ReML estimator (ReML).  Signal | Intra-correlation | $\rho_{\text{True}}$ | Type | RMSE | $\left\lvert\text{Bias}\right\rvert$ | SD
---|---|---|---|---|---|---
$k_{\eta}=1.5$ | $\phi_{\gamma}=1$ | 0.1 | ReML | 0.1550 | 0.0189 | 0.1547
FE | 0.2949 | 0.0570 | 0.2908
CA | 0.2837 | 0.0484 | 0.2809
0.35 | ReML | 0.1521 | 0.0014 | 0.1529
FE | 0.2724 | 0.0444 | 0.2701
CA | 0.2613 | 0.0454 | 0.2586
0.6 | ReML | 0.1104 | 0.0196 | 0.1092
FE | 0.2404 | 0.0726 | 0.2303
CA | 0.2315 | 0.0788 | 0.2188
$\phi_{\gamma}=0.25$ | 0.1 | ReML | 0.1961 | 0.0256 | 0.1954
FE | 0.2961 | 0.0242 | 0.2966
CA | 0.2657 | 0.0068 | 0.2669
0.35 | ReML | 0.1843 | 0.0230 | 0.1838
FE | 0.2965 | 0.1301 | 0.2678
CA | 0.3243 | 0.1720 | 0.2763
0.6 | ReML | 0.1456 | 0.0027 | 0.1463
FE | 0.3144 | 0.1959 | 0.2472
CA | 0.3591 | 0.2594 | 0.2496

(a) $\rho_{12}=0.1$

(b) $\rho_{13}=0.35$

(c) $\rho_{23}=0.6$

Figure 6: Empirical distribution of estimated inter-regional correlations with
very strong signal $k_{\eta}=1.5$: (a) Weak inter-regional correlation, (b)
Medium inter-regional correlation, (c) Strong inter-regional correlation. In
each subfigure: two intra-regional correlation rates are $\phi_{\gamma}=0.25$
(strong) and $\phi_{\gamma}=1$ (weak); three estimators are Pearson
correlation of average (CA), Pearson correlation of fixed-effect estimate
(FE), ReML estimator (ReML). Dashed line denotes the true value of inter-
regional correlation.

## References

  * Achard and Gannaz (2019) [author] Achard, SophieS. and Gannaz, IrèneI. (2019). Wavelet-based and Fourier-based multivariate Whittle estimation: multiwave. Journal of Statistical Software 89 1–31.
  * Achard et al. (2006) [author] Achard, SophieS., Salvador, RaymondR., Whitcher, BrandonB., Suckling, JohnJ. and Bullmore, EDE. (2006). A resilient, low-frequency, small-world human brain functional network with highly connected association cortical hubs. Journal of Neuroscience 26 63–72.
  * Achard et al. (2011) [author] Achard, SophieS., Coeurjolly, Jean-FrançoisJ.-F., Marcillaud, RomainR. and Richiardi, JonasJ. (2011). fMRI functional connectivity estimators robust to region size bias. In 2011 IEEE Statistical Signal Processing Workshop (SSP) 813–816. IEEE.
  * Afyouni, Smith and Nichols (2019) [author] Afyouni, SorooshS., Smith, Stephen MS. M. and Nichols, Thomas ET. E. (2019). Effective degrees of freedom of the Pearson’s correlation coefficient under autocorrelation. NeuroImage 199 609–625. https://doi.org/10.1016/j.neuroimage.2019.05.011
  * Azevedo et al. (2022) [author] Azevedo, TiagoT., Campbell, AlexanderA., Romero-Garcia, RafaelR., Passamonti, LucaL., Bethlehem, Richard AIR. A., Liò, PietroP. and Toschi, NicolaN. (2022). A deep graph neural network architecture for modelling spatio-temporal dynamics in resting-state functional MRI data. Medical Image Analysis 79 102471\. https://doi.org/10.1016/j.media.2022.102471
  * Bartolucci and Farcomeni (2022) [author] Bartolucci, FrancescoF. and Farcomeni, AlessioA. (2022). A spatio-temporal model based on discrete latent variables for the analysis of COVID-19 incidence. Spatial Statistics 49 100504\.
  * Becq et al. (2020) [author] Becq, GJCG., Habet, TarikT., Collomb, NoraN., Faucher, MargauxM., Delon-Martin, ChantalC., Coizet, VéroniqueV., Achard, SophieS. and Barbier, Emmanuel LE. L. (2020). Functional connectivity is preserved but reorganized across several anesthetic regimes. NeuroImage 219 116945\.
  * Benjamini and Yekutieli (2001) [author] Benjamini, YoavY. and Yekutieli, DanielD. (2001). The Control of the False Discovery Rate in Multiple Testing under Dependency. The Annals of Statistics 29 1165–1188.
  * Biswal et al. (1995) [author] Biswal, BharatB., Zerrin Yetkin, F.F., Haughton, Victor M.V. M. and Hyde, James S.J. S. (1995). Functional connectivity in the motor cortex of resting human brain using echo-planar MRI. Magnetic Resonance in Medicine 34 537-541. https://doi.org/10.1002/mrm.1910340409
  * Bowman (2007) [author] Bowman, F DuboisF. D. (2007). Spatiotemporal models for region of interest analyses of functional neuroimaging data. Journal of the American Statistical Association 102 442–453.
  * Chaimow et al. (2018) [author] Chaimow, DenisD., Yacoub, EssaE., Uğurbil, KâmilK. and Shmuel, AmirA. (2018). Spatial specificity of the functional MRI blood oxygenation response relative to neuronal activity. NeuroImage 164 32–47.
  * Cressie and Lahiri (1993) [author] Cressie, NoelN. and Lahiri, Soumendra NathS. N. (1993). The asymptotic distribution of REML estimators. Journal of Multivariate Analysis 45 217–233.
  * Cressie and Lahiri (1996) [author] Cressie, NoelN. and Lahiri, Soumendra NathS. N. (1996). Asymptotics for REML estimation of spatial covariance parameters. Journal of Statistical Planning and Inference 50 327–341.
  * Cui et al. (2020) [author] Cui, ZaixuZ., Li, HongmingH., Xia, Cedric H.C. H., Larsen, BartB., Adebimpe, AzeezA., Baum, Graham L.G. L., Cieslak, MattM., Gur, Raquel E.R. E., Gur, Ruben C.R. C., Moore, Tyler M.T. M., Oathes, Desmond J.D. J., Alexander-Bloch, Aaron F.A. F., Raznahan, ArminA., Roalf, David R.D. R., Shinohara, Russell T.R. T., Wolf, Daniel H.D. H., Davatzikos, ChristosC., Bassett, Danielle S.D. S., Fair, Damien A.D. A., Fan, YongY. and Satterthwaite, Theodore D.T. D. (2020). Individual Variation in Functional Topography of Association Networks in Youth. Neuron 106 340-353.e8. https://doi.org/10.1016/j.neuron.2020.01.029
  * Delicado et al. (2010) [author] Delicado, PedroP., Giraldo, RamónR., Comas, CarlosC. and Mateu, JorgeJ. (2010). Statistics for spatial functional data: some recent contributions. Environmetrics: The Official Journal of The International Environmetrics Society 21 224–239.
  * Deuker et al. (2009) [author] Deuker, LorenaL., Bullmore, Edward TE. T., Smith, MarieM., Christensen, SorenS., Nathan, Pradeep JP. J., Rockstroh, BrigitteB. and Bassett, Danielle SD. S. (2009). Reproducibility of graph metrics of human brain functional networks. Neuroimage 47 1460–1468.
  * Eickhoff, Yeo and Genon (2018) [author] Eickhoff, Simon BS. B., Yeo, BTB. and Genon, SarahS. (2018). Imaging-based parcellations of the human brain. Nature Reviews Neuroscience 19 672–686.
  * French and Kokoszka (2021) [author] French, Joshua PJ. P. and Kokoszka, Piotr SP. S. (2021). A sandwich smoother for spatio-temporal functional data. Spatial Statistics 42 100413\.
  * Gardner et al. (2018) [author] Gardner, JacobJ., Pleiss, GeoffG., Weinberger, Kilian QK. Q., Bindel, DavidD. and Wilson, Andrew GA. G. (2018). GPyTorch: Blackbox Matrix-Matrix Gaussian Process Inference with GPU Acceleration. In Advances in Neural Information Processing Systems (S.S. Bengio, H.H. Wallach, H.H. Larochelle, K.K. Grauman, N.N. Cesa-Bianchi and R.R. Garnett, eds.) 31\. Curran Associates, Inc.
  * Gräler, Pebesma and Heuvelink (2016) [author] Gräler, BenediktB., Pebesma, Edzer JE. J. and Heuvelink, Gerard BMG. B. (2016). Spatio-temporal interpolation using gstat. The R Journal 8 204\.
  * Gromenko et al. (2012) [author] Gromenko, OleksandrO., Kokoszka, PiotrP., Zhu, LieL. and Sojka, JanJ. (2012). Estimation and testing for spatially indexed curves with application to ionospheric and magnetic field trends. The Annals of Applied Statistics 6 669–696.
  * Harville (1974) [author] Harville, David AD. A. (1974). Bayesian inference for variance components using only error contrasts. Biometrika 61 383–385.
  * Hengl et al. (2012) [author] Hengl, TomislavT., Heuvelink, GerardG., Perčec Tadić, MelitaM. and Pebesma, Edzer JE. J. (2012). Spatio-temporal prediction of daily temperatures using time-series of MODIS LST images. Theoretical and Applied Climatology 107 265–277.
  * Hörmann and Kokoszka (2011) [author] Hörmann, SiegfriedS. and Kokoszka, PiotrP. (2011). Consistency of the mean and the principal components of spatially distributed functional data. In Recent Advances in Functional Data Analysis and Related Topics 169–175. Springer.
  * Jennrich and Schluchter (1986) [author] Jennrich, Robert IR. I. and Schluchter, Mark DM. D. (1986). Unbalanced repeated-measures models with structured covariance matrices. Biometrics 42 805–820. https://doi.org/10.2307/2530695
  * Kaiser (2011) [author] Kaiser, MarcusM. (2011). A tutorial in connectome analysis: topological and spatial features of brain networks. NeuroImage 57 892–907. https://doi.org/10.1016/j.neuroimage.2011.05.025
  * Krizhevsky, Sutskever and Hinton (2012) [author] Krizhevsky, AlexA., Sutskever, IlyaI. and Hinton, Geoffrey EG. E. (2012). ImageNet Classification with Deep Convolutional Neural Networks. In Advances in Neural Information Processing Systems (F.F. Pereira, C. J.C. J. Burges, L.L. Bottou and K. Q.K. Q. Weinberger, eds.) 25\. Curran Associates, Inc.
  * Lindstrom and Bates (1988) [author] Lindstrom, Mary JM. J. and Bates, Douglas MD. M. (1988). Newton-Raphson and EM algorithms for linear mixed-effects models for repeated-measures data. Journal of the American Statistical Association 83 1014–1022.
  * Logan and Rowe (2004) [author] Logan, Brent R.B. R. and Rowe, Daniel B.D. B. (2004). An evaluation of thresholding techniques in fMRI analysis. NeuroImage 22 95-108. https://doi.org/10.1016/j.neuroimage.2003.12.047
  * Lohmann et al. (2018) [author] Lohmann, GabrieleG., Stelzer, JohannesJ., Lacosse, EricE., Kumar, Vinod JV. J., Mueller, KarstenK., Kuehn, EstherE., Grodd, WolfgangW. and Scheffler, KlausK. (2018). LISA improves statistical analysis for fMRI. Nature Communications 9 1–9.
  * Lv et al. (2018) [author] Lv, H.H., Wang, Z.Z., Tong, E.E., Williams, L. M.L. M., Zaharchuk, G.G., Zeineh, M.M., Goldstein-Piekarski, A. N.A. N., Ball, T. M.T. M., Liao, C.C. and Wintermark, M.M. (2018). Resting-state functional MRI: everything that nonexperts have always wanted to know. American Journal of Neuroradiology 39 1390–1399. 10.3174/ajnr.A5527
  * McIntosh, Chau and Protzner (2004) [author] McIntosh, Anthony RA. R., Chau, WKW. and Protzner, Andrea BA. B. (2004). Spatiotemporal analysis of event-related fMRI data using partial least squares. NeuroImage 23 764–775. https://doi.org/10.1016/j.neuroimage.2004.05.018
  * Meskaldji et al. (2011) [author] Meskaldji, Djalel EddineD. E., Ottet, Marie-ChristineM.-C., Cammoun, LeilaL., Hagmann, PatricP., Meuli, RetoR., Eliez, StephanS., Thiran, Jean PhilippeJ. P. and Morgenthaler, StephanS. (2011). Adaptive strategy for the statistical analysis of connectomes. PloS One 6 e23009. https://doi.org/10.1371/journal.pone.0023009
  * Michon et al. (2022) [author] Michon, Katherine J.K. J., Khammash, DaliaD., Simmonite, MollyM., Hamlin, Abbey M.A. M. and Polk, Thad A.T. A. (2022). Person-specific and precision neuroimaging: Current methods and future directions. NeuroImage 263 119589\. https://doi.org/10.1016/j.neuroimage.2022.119589
  * Moghimi et al. (2022) [author] Moghimi, PanteaP., Dang, Anh TheA. T., Do, QuanQ., Netoff, Theoden I.T. I., Lim, Kelvin O.K. O. and Atluri, GowthamG. (2022). Evaluation of functional MRI-based human brain parcellation: a review. Journal of Neurophysiology 128 197-217. PMID: 35675446. 10.1152/jn.00411.2021
  * Nocedal (1980) [author] Nocedal, JorgeJ. (1980). Updating quasi-Newton matrices with limited storage. Mathematics of Computation 35 773–782. https://doi.org/10.1090/S0025-5718-1980-0572855-7
  * Park and Fiecas (2022) [author] Park, Jun YoungJ. Y. and Fiecas, MarkM. (2022). CLEAN: Leveraging spatial autocorrelation in neuroimaging data in clusterwise inference. NeuroImage 255 119192\.
  * Percival and Walden (2006) [author] Percival, Donald BD. B. and Walden, Andrew TA. T. (2006). Wavelet Methods for Time Series Analysis 4\. Cambridge University Press.
  * Petersen and Pedersen (2008) [author] Petersen, Kaare BrandtK. B. and Pedersen, Michael SyskindM. S. (2008). The Matrix Cookbook. Technical University of Denmark 7 510\.
  * Petersen et al. (2016) [author] Petersen, AlexanderA., Zhao, JianyangJ., Carmichael, OwenO. and Müller, Hans-GeorgH.-G. (2016). Quantifying individual brain connectivity with functional principal component analysis for networks. Brain Connectivity 6 540–547.
  * Pinheiro and Bates (2006) [author] Pinheiro, JoséJ. and Bates, DouglasD. (2006). Mixed-effects models in S and S-PLUS. Springer Science & Business Media.
  * Powell (2009) [author] Powell, Michael JDM. J. (2009). The BOBYQA algorithm for bound constrained optimization without derivatives. Cambridge NA Report NA2009/06, University of Cambridge, Cambridge 26\.
  * Richiardi et al. (2013) [author] Richiardi, JonasJ., Achard, SophieS., Bunke, HorstH. and Van De Ville, DimitriD. (2013). Machine learning with brain graphs: predictive modeling approaches for functional imaging in systems neuroscience. IEEE Signal Processing Magazine 30 58–70.
  * Saatchi (2011) [author] Saatchi, YunusY. (2011). Scalable inference for structured Gaussian process models, PhD thesis, University of Cambridge.
  * Stein (1999) [author] Stein, Michael LM. L. (1999). Interpolation of Spatial Data: Some Theory for Kriging. Springer Science & Business Media.
  * Taylor and Worsley (2007) [author] Taylor, Jonathan EJ. E. and Worsley, Keith JK. J. (2007). Detecting sparse signals in random fields, with an application to brain mapping. Journal of the American Statistical Association 102 913–928.
  * Termenon et al. (2016) [author] Termenon, MaïtéM., Jaillard, AssiaA., Delon-Martin, ChantalC. and Achard, SophieS. (2016). Reliability of graph analysis of resting state fMRI using test-retest dataset from the Human Connectome Project. NeuroImage 142 172–187.
  * Van Den Heuvel and Pol (2010) [author] Van Den Heuvel, Martijn PM. P. and Pol, Hilleke E HulshoffH. E. H. (2010). Exploring the brain network: a review on resting-state fMRI functional connectivity. European Neuropsychopharmacology 20 519–534.
  * Vértes et al. (2012) [author] Vértes, Petra EP. E., Alexander-Bloch, Aaron FA. F., Gogtay, NitinN., Giedd, Jay NJ. N., Rapoport, Judith LJ. L. and Bullmore, Edward TE. T. (2012). Simple models of human brain functional networks. Proceedings of the National Academy of Sciences 109 5868–5873.
  * Whitcher, Guttorp and Percival (2000) [author] Whitcher, B.B., Guttorp, P.P. and Percival, D. B.D. B. (2000). Wavelet analysis of covariance with application to atmospheric time series. Journal of Geophysical Research 105(D11) 941-962.
  * Woolrich et al. (2004) [author] Woolrich, Mark WilliamM. W., Jenkinson, MarkM., Brady, J MichaelJ. M. and Smith, Stephen MS. M. (2004). Fully Bayesian spatio-temporal modeling of fMRI data. IEEE Transactions on Medical Imaging 23 213–231.
  * Yarger, Stoev and Hsing (2022) [author] Yarger, DrewD., Stoev, StilianS. and Hsing, TailenT. (2022). A functional-data approach to the Argo data. The Annals of Applied Statistics 16 216–246.
  * Zalesky, Fornito and Bullmore (2012) [author] Zalesky, AndrewA., Fornito, AlexA. and Bullmore, EdE. (2012). On the use of correlation as a measure of network connectivity. NeuroImage 60 2096–2106.
  * Zhang and Li (2022) [author] Zhang, HaozheH. and Li, YehuaY. (2022). Unified principal component analysis for sparse and dense functional data under spatial dependency. Journal of Business & Economic Statistics 40 1523-1537. 10.1080/07350015.2021.1938085
  * Zhang et al. (2014) [author] Zhang, LinlinL., Guindani, MicheleM., Versace, FrancescoF. and Vannucci, MarinaM. (2014). A spatio-temporal nonparametric Bayesian variable selection model of fMRI data for clustering correlated time courses. NeuroImage 95 162–175.