# An efficient approach to characterize spatio-temporal dependence in cortical
surface fMRI data

Huy Dang
Dept. of Statistics, Pennsylvania State University, USA
and
Marzia A. Cremona
Dept. of Operations and Decision Systems, Université Laval, Canada
CHU de Québec – Université Laval Research Center, Canada
and
Francesca Chiaromonte
Dept. of Statistics, Pennsylvania State University, USA
Inst. of Economics and L’EMbeDS, Sant’Anna School of Advanced Studies, Italy
and
Nicole Lazar
Dept. of Statistics, Pennsylvania State University, USA

###### Abstract

Functional magnetic resonance imaging (fMRI) is a neuroimaging technique known
for its ability to capture brain activity non-invasively and at fine spatial
resolution (1-3mm). Cortical surface fMRI (cs-fMRI) is a recent development of
fMRI that restricts attention to signals exclusively from tissue types that
have neuronal activities, as opposed to the whole brain. cs-fMRI data is
plagued with non-stationary spatial correlations and long temporal dependence
which, if inadequately accounted for, can hinder various types of downstream
statistical analyses. We propose a fully integrated approach that captures
both spatial non-stationarity and varying ranges of temporal dependence across
regions of interest. More specifically, we impose non-stationary spatial
priors on the latent activation fields and model temporal dependence via
fractional Gaussian errors of varying Hurst parameters, which can be studied
through a wavelet transformation and its coefficients’ variances at different
scales. We demonstrate the performance of our proposed approach via
simulations and an application to a visual working memory task cs-fMRI
dataset.

Keywords: cs-fMRI, spatio-temporal dependence, SPDE, wavelets, fractional
Gaussian noise.

##  1 Introduction and motivation

Functional magnetic resonance imaging (fMRI) is a neuroimaging technique known
for its ability to non-invasively measure Blood Oxygen Level Dependent (BOLD)
signals at fine spatial resolution (typically 2-3mm). The BOLD signal is
defined as the changes in the ratio of oxygenated to deoxygenated blood,
either in response to a task/stimulus in an experiment or as a consequence of
spontaneous neural metabolism. In fMRI literature, the BOLD signal is widely
taken to be a proxy of brain activity, and used to capture associations
between local brain areas with different functions. Traditionally, fMRI data
collected during a task are represented by 4D arrays: signals are measured at
discrete locations, each representing a small cube in a partition of the 3D
brain (these are known as voxels), along a time course – which produces the
4th dimension. This type of “volumetric” fMRI data contains both signals from
tissues that have neuronal activities (such as gray matter) and noise from
tissues that do not (such as white matter and cerebral spinal fluid). Because
of spatial contiguity between active and inactive tissues, distance-based
analyses are often affected by spurious noise from the latter. Cortical
surface fMRI data (cs-fMRI) is the result of a recently developed technique
that uses vertices on a 2D surface mesh, instead of traditional volumetric
voxels, to represent the folded, sheet-like geometry of the cerebral cortex.
In multi-subject studies surface meshes from different individuals can be
aligned to a common template based on cortical folding patterns and areal
features (Glasser et al., 2013, 2016). By virtue of being confined exclusively
to cortical gray matter, cs-fMRI data is not affected by signal contamination
from inactive tissues. Increased homogeneity due to the fact that data are
collected only from cortical gray matter thus improves distance-based
analyses. Figure 1 illustrates the surface mesh used to achieve the 2D
representation of the cerebral cortex at different degrees of “inflation” –
ranging from “white matter surface” (least inflated and preserving folding
patterns) to “inflated”, “very inflated”, and “spherical” (most inflated, no
folding patterns). By inflating the surface mesh, distances between mesh
vertices lose their original 3D Euclidean interpretation and become more
geodesic-like. Readers interested in a more comprehensive review of volumetric
and cs-fMRI data are referred to Mejia et al. (2020). From here on, “voxel”
refers to the smallest spatial unit in volumetric fMRI data, and “vertex” to
its cs-fMRI counterpart.

|
---|---
|
Fig. 1:  The 2D cortical surface mesh used in cs-fMRI data, visualized with
increasing degrees of inflation: white matter surface, inflated, very
inflated, spherical. The mesh is colored according to a parcellation of the
brain into 100 regions (50 regions per hemisphere) introduced in Schaefer et
al. (2018). Subcortical regions are plotted in white and are not analyzed.

###  1.1 Temporal dependence

It is well known that fMRI data exhibit temporal correlations. Existing
literature provides evidence of both short- and long-range temporal dependence
that can be attributed to spontaneous, non task-related fluctuations in
neuronal activity (see e.g. Bullmore et al., 2003; Park et al., 2010; Li &
Huang, 2014). Moreover, Park et al. (2010); Li & Huang (2014) showed that the
ranges of dependence can be affected by brain structures and functions.

Fig. 2: Autocorrelation function (left), partial autocorrelation function
(middle) and spectrum-frequency plot (right) for time series with short-range
(top) and long-range dependence (bottom).

As a real data illustration, Figure 2 shows the autocorrelation function
(ACF), the partial autocorrelation function (PACF) and the spectrum against
frequency for two resting state BOLD time series. These are taken at different
brain locations from an individual in the Human Connectome Project (HCP)
healthy young adults data set (Van Essen et al., 2013). The upper panels
illustrate a time series with short-range temporal dependence, as evidenced by
the presence of only a few significant ACF and PACF coefficients at small
lags, while the bottom panels illustrate a time series with long-range
temporal dependence, with a slowly decaying autocorrelation. The latter also
shows very large estimated spectral density at small frequencies, in contrast
to the more evenly distributed spectral density across frequencies of the
former. Since the spectral density is the Fourier transform of the covariance
function, large spectral density estimates corresponds to large contributions
to the covariance function’s second moment, implying significant dependence at
frequencies near $0$, or equivalently, at large lags.

The vast majority of fMRI literature employs autoregressive models of order
$p$ – that we denote as $AR(p)$ – to prewhiten the time series separately at
each voxel, or vertex, prior to further analyses. Some recent examples include
the use of $AR(1)$ in Derado et al. (2010), $AR(2)$ in Castruccio et al.
(2018), $AR(3)$ in Risk et al. (2016), and $AR(6)$ in Mejia et al. (2020).
However, ordinary $AR(p)$ models are designed to capture autocorrelations up
to a finite order $p$, and therefore are not suitable for modeling long-range
dependence. In addition, we show with a simple simulation that prewhitening
time series with an AR process may cause loss of information and lead to
biased estimates. Let
$\boldsymbol{x}=\\{x_{t}\\}=\int_{-\infty}^{\infty}h(u)s(t-u)du$ be a stimulus
$h(\cdot)$ convolved with the canonical hemodynamic response function
$s(\cdot)$ observed along times, say, $t=1,\ldots,256$ (more detail on this
convolutional representation, which is typical in fMRI data, will be provided
in Section 3). We simulate the responses $\boldsymbol{y}$ and $\boldsymbol{z}$
by adding to $\boldsymbol{x}$ a long-range dependent noise
$\boldsymbol{\epsilon}$ and an auto-regressive noise $\boldsymbol{\eta}$,
respectively; that is

$\displaystyle\boldsymbol{y}=\beta\boldsymbol{x}+\boldsymbol{\epsilon}\;,\;\boldsymbol{z}=\beta\boldsymbol{x}+\boldsymbol{\eta}\
.$

To simulate the long-range dependent $\boldsymbol{\epsilon}$, we use
fractional Gaussian noise with Hurst parameter $H=0.8$. For the auto-
regressive $\boldsymbol{\eta}$, we fix $p=6$ and use AR coefficients estimated
from actual BOLD signals (from the same HCP individual used above, at randomly
selected vertices/locations). In both cases, we set $\beta=2$ and noise
variance equal to $1$. We then carry out the estimation of $\beta$ following
the procedure used in Mejia et al. (2020): we fit $AR(6)$ models to both
$\boldsymbol{\epsilon}$ and $\boldsymbol{\eta}$, obtain the AR coefficients in
each of the two cases and use them to pre-whiten $\boldsymbol{y}$ and
$\boldsymbol{z}$, and finally regress the pre-whitened time series
$\boldsymbol{y}_{pw}$ and $\boldsymbol{z}_{pw}$ on $\boldsymbol{x}$ to
estimate $\beta$.

Dependence | Method | Min | Q1 | Median | Mean | Q3 | Max | SD
---|---|---|---|---|---|---|---|---
| Ordinary | 1.189 | 1.831 | 2.015 | 2.008 | 2.186 | 3.010 | 0.263
Long | $AR(6)$ prewhitening | 0.688 | 1.111 | 1.242 | 1.239 | 1.353 | 1.8285 | 0.182
| Ordinary | 0.482 | 1.770 | 2.013 | 2.009 | 2.244 | 2.997 | 0.343
$AR(6)$ | $AR(6)$ prewhitening | 0.476 | 1.106 | 1.278 | 1.283 | 1.455 | 2.025 | 0.246
Table 1: Summaries for $\hat{\beta}$ estimates in long- and short-range
dependence scenarios, obtained with $AR(6)$ pre-whitening and with an ordinary
regression without pre-whitening (results based on 1000 simulation runs).

Table 1 shows estimation results in the long- and short-range dependence
scenarios with the $AR(6)$ prewhitening scheme, as well as with an ordinary
regression without prewhitening. The results, based on 1000 simulation runs,
indicate that ordinary regression (without prewhitening) produces unbiased
estimates of $\beta$. On the contrary, $AR(6)$ prewhitening of the time series
produces a marked underestimation of $\beta$ in both scenarios, including when
the noise is in fact generated with an $AR(6)$ process.

This simple but informative simulation exercise illustrates the problems that
can be induced by prewhitening time series with autoregressive models prior to
further analyses. In order to account for temporal correlations, we strive for
a more flexible modeling approach that can accommodate both short- and long-
range dependence across brain regions. Such approach should not cause loss of
information or biases, and should be implementable through a computationally
efficient algorithm. In particular, to contain computational cost, it should
produce sparse temporal correlations.

###  1.2 Spatial dependence

The sheer size of fMRI data presents a challenge when modeling correlations
between voxels or vertices. In addition to sparsity in representing temporal
dependences, sparsity assumptions on the spatial dependence structure are
absolutely critical; without them, most standard analyses (e.g., fitting
linear models) will face an immense computational burden due to the need to
invert a dense $VT\times VT$ covariance matrix, where $V$ is the number of
voxels or vertices, and $T$ is the number of time points. In traditional
volumetric fMRI data $V\approx 120,000$, whereas in cs-fMRI $V\approx 30,000$
for each hemisphere, and $T$ typically ranges from $200$ to $1,200$. In
existing fMRI literature, the most common (and unrealistic) assumption is that
of no spatial dependence, which means fitting a model for each time series at
each voxel/vertex independently – an approach that is often referred to as
massively univariate  (Lazar, 2008). However, untreated spatial correlations
often lead to underestimated uncertainty and inflated type I error rates in
hypothesis tests (Castruccio et al., 2018).

While inverting the aforementioned covariance matrix has been and remains
computationally unrealistic, recent years have seen a rise in efforts to
acknowledge and treat spatial dependence. This is typically achieved via a
combination of more sophisticated statistical models, approximation
techniques, and/or down-sampling of the data. An early example of spatial
dependence treatment is the work of Bowman (2005), employing a two-stage
hierarchical Bayesian approach. First, standard regressions were fitted
independently at each voxel to estimate activation coefficients; then,
coefficients of voxels belonging to the same functional region were used to
model within region correlations via a spatial autoregressive model.
Computation was made feasible by subsampling the 3D volumetric data to 2D
slice data, and by assuming stationarity of the covariance. A different
approach was taken by Kang et al. (2013), which applied a Fourier transform to
the fMRI time series and fitted a spatio-spectral mixed-effect model, using
random effects to capture spatial correlation within and between regions of
interest. Although some degrees of non-stationarity was allowed in this
approach, activation was restricted to regional level instead of individual
voxels. More recently, Castruccio et al. (2018) used $L$ Gaussian processes
with different Matern covariances to model local dependence within different
regions, and a region-specific random effect to account for between-region
dependence. While the model incorporated non-stationary spatial dependence on
full 3D volumetric data, estimation at each spatial scale (local and regional)
was carried out sequentially and only a portion of the data was used at each
estimation step. For example, local dependence was estimated region by region,
and regional dependence was estimated using only regional averages. On a
different front, Mejia et al. (2020) assumed a latent task activation model
and represented spatial dependence via spatial priors on the activation
fields. The proposed method was computationally very efficient, as estimation
was carried out within the Gaussian Markov Random Fields (GMRFs) framework,
and inference was made possible in continuous space thanks to an established
connection between GMRFs and Gaussian Fields (GFs) with Matern covariance
functions (Lindgren et al., 2011). Their approach, however, did not allow for
non-stationarity, which is an important spatial feature of fMRI data.

###  1.3 On things to come

To improve on the existing treatment of fMRI data, we propose a fully
integrated approach that captures both spatial non-stationarity and varying
ranges of temporal dependence across regions of interest, focusing in
particular on the new class of cs-fMRI data. More specifically, we impose
spatial priors on the latent activation fields as in Mejia et al. (2020);
however, our approach allows for non-stationarity by letting the prior
hyperparameters be driven by local smoothness in the data. We model temporal
dependence using fractional Gaussian errors of varying Hurst parameters; the
Hurst parameters can then be studied through wavelet transformation and the
wavelet coefficients’ variances at different scales. Bayesian inference is
carried out approximating the marginal posteriors with the Integrated Nested
Laplace Approximation (INLA) approach (Rue & Martino, 2009).

The remainder of this article is organized as follows. Section 2 provides the
theoretical background. Section 3 details our modeling approach. Section 4
demonstrates the performance of our proposal through simulations and an
application to real cs-fMRI data concerning a visual working memory task.
Section 5 contains final remarks.

##  2 Theoretical background

###  2.1 Fractional Gaussian noise and wavelet transformation

Long- and short-range dependence can be flexibly and sparsely modeled with
fractional Gaussian noise (fGn). fGn is defined as the difference between
consecutive values of fractional Brownian motion (Fadili & Bullmore, 2002).
Given a lag $\ell>0$, the autocovariance of a fractional Gaussian process $G$
is given by

$\displaystyle
C_{G}(\ell;H)=\frac{\sigma^{2}}{2}[(\ell+1)^{2H}-2\ell^{2H}+(\ell-1)^{2H}]\underset{\ell\rightarrow\infty}{\approx}\sigma^{2}H(2H-1)\ell^{2H-2}$
(1)

where $H$ is the Hurst parameter and $\sigma^{2}$ is the variance. It is easy
to see that when $H=\frac{1}{2}$, the process corresponds to white noise. When
$0<H<\frac{1}{2}$, or equivalently, $-2<2H-2<-1$, the autocovariance decays
exponentially fast with lag $\ell$ and the process has short-range dependence.
On the other hand, when $\frac{1}{2}<H<1$, the autocovariance decays at
hyperbolic rate $-1<2H-2<0$, causing
$\sum_{\ell=0}^{\infty}C_{G}(\ell;H)=\infty$, and the process has long-range
dependence. When $\sum_{\ell=0}^{\infty}C_{G}(\ell;H)=\infty$, individual
correlations at large lags may be small, but cumulatively, their effect is
significant. A consequence of untreated correlations are biases in variance
estimates that do not vanish as the sample size increases, which in turn
affects confidence intervals and hypothesis tests.

To treat time series with varying ranges of dependence, one efficient approach
is wavelet transformation. It has the desirable property of representing an
autocorrelated process with coefficients that are approximately uncorrelated
(Flandrin, 1992; Fadili & Bullmore, 2002). In addition, there is no loss of
information as wavelet bases are orthonormal, and the transformation is
calculated with $\mathcal{O}(n\log_{2}n)$ operations – comparable to fast
Fourier transform (Mallat, 2009). As detailed below, the Hurst parameter can
be estimated from the decaying rate of the wavelet coefficients’ variances
across scales.

Let $(\epsilon_{1},\ldots,\epsilon_{n})$ be a process with covariance function
given by Equation 1. The discrete wavelet transform decomposes the process
into a set of detail coefficients $\\{d_{jk}:j=1,\dots,J,\
k=1,\ldots,n/2^{j}\\}$, and approximation coefficients
$\\{a_{Jk}:k=1,\ldots,n/2^{J}\\}$, where $J$ is the coarsest level of the
transformation. Let $\gamma=2H-1$, and $c_{\gamma}=(2\pi)^{-2H}\sin{(\pi
H)}\Gamma(2H+1)$. According to Flandrin (1992); Fadili & Bullmore (2002), the
covariance matrix of the wavelet transformed process is approximately
diagonal, and the detail coefficients $\\{d_{jk}:k=1,\ldots,n/2^{j}\\}$ and
the approximation coefficients $\\{a_{Jk}:k=1,\ldots,n/2^{J}\\}$ have
variances

$\displaystyle
S_{d_{j}}\approx\frac{\sigma^{2}c_{\gamma}2^{j\gamma}(2-2^{\gamma})}{(2\pi)^{\gamma}(1-\gamma)}\;,\;S_{a_{J}}\approx\frac{\sigma^{2}c_{\gamma}2^{(J+1)\gamma}}{(2\pi)^{\gamma}(1-\gamma)}\;.$
(2)

Thus, $\gamma$, or equivalently the Hurst parameter $H$, can be estimated
through the linear relationship $\log_{2}(S_{d_{j}})=\gamma j+\mathcal{O}(1)$.

Fig. 3: Autocorrelation function (left), partial autocorrelation function
(middle) and spectrum-frequency plot (right) for a time series with long-range
dependence (top) and its discrete wavelet transform (bottom).

Figure 3 demonstrates the decorrelating effect of the discrete wavelet
transform on a time series with long range dependence. Compared to the
original time series (top row), its wavelet transformation (bottom row) has
insignificant autocorrelation coefficients, and much more evenly distributed
spectral density across frequencies.

###  2.2 The link between GMRFs and GFs

In most analyses of interest, estimation requires calculations with precision
matrices ($\boldsymbol{Q}=\boldsymbol{\Sigma}^{-1}$), not covariance matrices
($\boldsymbol{\Sigma}$). For example, if
$\boldsymbol{\beta}=(\beta_{1},\cdots,\beta_{V})\sim
N(\boldsymbol{\mu},\boldsymbol{Q}^{-1})$, then the likelihood is

$L(\theta|\boldsymbol{\beta})=\sqrt{(2\pi)^{-V}|\boldsymbol{Q}|}\;\exp{\bigg{(}-\frac{(\boldsymbol{\beta}-\boldsymbol{\mu})^{\top}\boldsymbol{Q}(\boldsymbol{\beta}-\boldsymbol{\mu})}{2}\bigg{)}}.$

Thus, GMRFs’ sparse, banded precision matrix structure can significantly speed
up model estimation. By definition of GMRFs, such structure is a result of the
conditional independence assumption. Note that while it is entirely
unrealistic to assume any two data points to be independent, it is more
reasonable to assume any two data points outside of their respective
neighborhoods to be independent conditional on all others. In this case, a
sparse, banded precision matrix still has a dense inverse, i.e. a dense
covariance matrix. Even though estimation can be performed speedily, there
remains a serious challenge. Most GMRFs’ precision matrices do not have an
established correspondence with a closed-form covariance function. Real life
data, including fMRI, are observed discretely from continuous processes; and
extrapolating beyond the discrete grid of observations requires knowing the
functional form of the covariance. For example, for an unobserved vector
$\boldsymbol{\beta}^{*}$, the conditional expectation
$E(\boldsymbol{\beta}^{*}|\boldsymbol{\beta})=\boldsymbol{\mu}_{\beta^{*}}-\boldsymbol{\Sigma}_{\beta^{*}\beta}\boldsymbol{Q}_{\beta\beta}(\boldsymbol{\beta}-\boldsymbol{\mu}_{\beta})$
is not available without knowing how to relate $\boldsymbol{Q}_{\beta\beta}$
to $\boldsymbol{\Sigma}_{\beta^{*}\beta}$.

Lindgren et al. (2011) solved this problem for a class of covariance
functions, namely the stationary Matérn covariances. Their solution relies on
a result by Whittle (1954, 1963) that links such class to the solutions of a
particular stochastic partial differential equation (SPDE). Specifically, the
stationary solutions $\beta(\boldsymbol{s})$ to the linear fractional SPDE

$\displaystyle(\kappa^{2}-\Delta)^{\frac{\nu+d/2}{2}}\beta(\boldsymbol{s})=\mathcal{W}(\boldsymbol{s}),\quad\boldsymbol{s}\in\mathbb{R}^{d}$
(3)

have Matérn covariance functions

$C(\beta(\boldsymbol{0}),\beta(\boldsymbol{s}))=\frac{\sigma^{2}}{2^{\nu-1}\Gamma(\nu)}(\kappa\left\lVert
s\right\rVert)^{\nu}K_{\nu}(\kappa\left\lVert s\right\rVert)$

where $K_{\nu}(\cdot)$ is the modified Bessel function of the second kind and
order $\nu$, $\Delta=\sum_{i=1}^{d}\partial^{2}/\partial s_{i}^{2}$ is the
Laplacian operator, $\mathcal{W}$ is Gaussian noise with unit variance,
$\kappa$ is the spatial scale parameter, $\alpha=\nu+d/2$ controls the
smoothness, and

$\displaystyle\sigma^{2}=\frac{\Gamma(\nu)}{\Gamma(\nu+d/2)(4\pi)^{d/2}\kappa^{2\nu}}.$

To complete their proposal, Lindgren et al. (2011) formulated GMRF
representations of the solutions to the SPDE in Equation 3. The precision
matrices $\boldsymbol{Q}$ of such GMRF representations are explicitly
parameterized by the SPDE parameters, or equivalently by the stationary Matérn
covariance parameters. We start by giving a finite dimensional representation
of the solutions to the SPDE; in symbols

$\displaystyle\beta(\boldsymbol{s})=\sum_{l=1}^{L}\psi_{l}(\boldsymbol{s})w_{l}$
(4)

for some deterministic basis functions $\\{\psi_{l}\\}$ and weights $w_{l}$.
In our applications, the bases are piecewise linear, constructed by
partitioning the domain into a set of non-intersecting triangles (i.e. a
triangular mesh). Specifically, given a triangle with vertices $i,j$ and $k$
located at $\boldsymbol{s}_{i},\boldsymbol{s}_{j}$ and $\boldsymbol{s}_{k}$ in
$\mathbb{R}^{2}$, for any inner location $\boldsymbol{s}$,
$\psi_{l}(\boldsymbol{s})$ is equal to the area of the triangle formed by
$\boldsymbol{s}$, $\boldsymbol{s}_{j},\boldsymbol{s}_{k}$ divided by the area
of the original triangle. It follows that $\psi_{i}$ takes a value of $1$ at
the $i^{th}$ vertex and a value of $0$ at all other vertices. Let
$\boldsymbol{C},\boldsymbol{G}$, $\boldsymbol{K}$ and
$\boldsymbol{Q_{\alpha}}$ be $L\times L$ matrices such that
$C_{ij}=\langle\psi_{i}\psi_{j}\rangle$,
$G_{ij}=\langle\nabla\psi_{i}\nabla\psi_{j}\rangle$,
$K_{ij}=\kappa^{2}C_{ij}+G_{ij}$ and

$\displaystyle\boldsymbol{Q}_{\alpha}=\begin{cases}\boldsymbol{K},\;\text{if
}\alpha=1\\\ \boldsymbol{KC}^{-1}\boldsymbol{K},\;\text{if }\alpha=2\\\
\boldsymbol{KC}^{-1}\boldsymbol{Q}_{\alpha-2}\boldsymbol{C}^{-1}\boldsymbol{K},\;\text{if
}\alpha=3,4,...\end{cases}$ (5)

Then, if $\boldsymbol{w}=\\{w_{l}\\}\sim N(0,\boldsymbol{Q}^{-1}_{\alpha})$,
the finite dimensional representations of the solutions to the SPDE in
Equation 3 are GFs with precisions $\boldsymbol{Q}_{\alpha}$. If
$\boldsymbol{C}$ is replaced by the diagonal matrix $\tilde{\boldsymbol{C}}$
where $\tilde{\boldsymbol{C}}_{ii}=\langle\psi_{i},1\rangle$, we obtain GMRF
representations instead of GFs. Note that $\nu$ (and hence $\alpha$) is
usually fixed since it is typically not identifiable in applications. Thus, if
$\boldsymbol{s}\in\mathbb{R}^{2}$, fixing $\nu=1,2,...$ gives $\alpha=2,3,...$

###  2.3 Integrated Nested Laplace Approximation

For Bayesian hierarchical models that involve latent Gaussian fields, closed-
form posterior distributions are, in general, unavailable. Sampling-based
Markov-chain Monte Carlo (MCMC) methods are used instead of seeking analytical
solutions. While such methods are asymptotically exact, in reality,
computational cost and time may prevent one from achieving asymptotic
guarantees; MCMC samples thus remain just repeated approximations of posterior
distributions. The issue is particularly marked when MCMC methods are applied
to latent Gaussian models; despite their versatility, they can be extremely
slow and the samples may not converge (Rue & Martino, 2009). To address this,
Rue & Martino (2009) introduced the Integrated Nested Laplace Approximation
(INLA), a method that performs direct numerical approximation of posterior
distributions instead of sampling from them. Given observed data
$\boldsymbol{y}$, the latent Gaussian field
$\boldsymbol{\beta}=\\{\beta_{v},v=1,\ldots,V\\}$ and a hyperparameter vector
$\boldsymbol{\theta}$, the joint posterior distribution of
$\boldsymbol{\beta}$ and $\boldsymbol{\theta}$ is

$\displaystyle\pi(\boldsymbol{\beta},\boldsymbol{\theta}|\boldsymbol{y})$
$\displaystyle\propto\pi(\boldsymbol{\theta})\pi(\boldsymbol{\beta}|\boldsymbol{\theta})\pi(\boldsymbol{y}|\boldsymbol{\beta},\boldsymbol{\theta})$
$\displaystyle\propto\pi(\boldsymbol{\theta})|\boldsymbol{Q}(\boldsymbol{\theta})|^{1/2}\exp{\left[-\frac{1}{2}\boldsymbol{\beta}^{\top}\boldsymbol{Q}(\boldsymbol{\theta})\boldsymbol{\beta}+\log{(\pi(\boldsymbol{y}|\boldsymbol{\beta},\boldsymbol{\theta}))}\right]}\
.$

Here, the goal is to approximate the posterior marginals

$\displaystyle\pi(\beta_{v}|\boldsymbol{y})$
$\displaystyle=\int\pi(\beta_{v},\boldsymbol{\theta}|\boldsymbol{y})d\boldsymbol{\theta}=\int\pi(\beta_{v}|\boldsymbol{\theta},\boldsymbol{y})\pi(\boldsymbol{\theta}|\boldsymbol{y})d\boldsymbol{\theta}\propto\int\frac{\pi(\boldsymbol{\beta},\boldsymbol{\theta}|\boldsymbol{y})}{\pi(\boldsymbol{\beta}_{-v}|\beta_{v},\boldsymbol{\theta},\boldsymbol{y})}\frac{\pi(\boldsymbol{\beta},\boldsymbol{\theta}|\boldsymbol{y})}{\pi(\boldsymbol{\beta}|\boldsymbol{\theta},\boldsymbol{y})}d\boldsymbol{\theta}$
$\displaystyle\pi(\theta_{l}|\boldsymbol{y})$
$\displaystyle=\int\pi(\boldsymbol{\theta}|\boldsymbol{y})d\boldsymbol{\theta}_{-l}\propto\int\frac{\pi(\boldsymbol{\beta},\boldsymbol{\theta}|\boldsymbol{y})}{\pi(\boldsymbol{\beta}|\boldsymbol{\theta},\boldsymbol{y})}d\boldsymbol{\theta}_{-l}$

where $\boldsymbol{\beta}_{-v}=\boldsymbol{\beta}\setminus\\{\beta_{v}\\}$ and
$\boldsymbol{\theta}_{-l}=\boldsymbol{\theta}\setminus\\{\theta_{l}\\}$. The
computational cost of numerical integration increases significantly with the
dimension of the hyperparameter vector $\boldsymbol{\theta}$; INLA can
accommodate up to about $10$ hyperparameters, but becomes unfeasible in higher
dimension. For a given $\boldsymbol{\theta}$, the denominator densities
$\pi(\boldsymbol{\beta}_{-v}|\beta_{v},\boldsymbol{\theta},\boldsymbol{y})$
and $\pi(\boldsymbol{\beta}|\boldsymbol{\theta},\boldsymbol{y})$ are
approximated by
$\tilde{\pi}(\boldsymbol{\beta}_{-v}|\beta_{v},\boldsymbol{\theta},\boldsymbol{y})$
and $\tilde{\pi}(\boldsymbol{\beta}|\boldsymbol{\theta},\boldsymbol{y})$,
respectively, using Laplace approximation. The marginals
$\pi(\beta_{v}|\boldsymbol{y})$ and $\pi(\theta_{l}|\boldsymbol{y})$ are then
approximated by numerical integration over $\boldsymbol{\theta}$; that is

$\displaystyle\pi(\beta_{v}|\boldsymbol{y})$
$\displaystyle\approx\int\frac{\pi(\boldsymbol{\beta},\boldsymbol{\theta}|\boldsymbol{y})}{\tilde{\pi}(\boldsymbol{\beta}_{-v}|\beta_{v},\boldsymbol{\theta},\boldsymbol{y})}\frac{\pi(\boldsymbol{\beta},\boldsymbol{\theta}|\boldsymbol{y})}{\tilde{\pi}(\boldsymbol{\beta}|\boldsymbol{\theta},\boldsymbol{y})}d\boldsymbol{\theta}$
$\displaystyle\pi(\theta_{l}|\boldsymbol{y})$
$\displaystyle\approx\int\frac{\pi(\boldsymbol{\beta},\boldsymbol{\theta}|\boldsymbol{y})}{\tilde{\pi}(\boldsymbol{\beta}|\boldsymbol{\theta},\boldsymbol{y})}d\boldsymbol{\theta}_{-l}\
.$

The computational speed-up provided by INLA with respect to MCMC algorithms is
on the scale of seconds/minutes compared to hours/days, with comparable
approximation errors. For more details, see Rue & Martino (2009).

##  3 Modeling approach

A standard model for an fMRI data set is

$\displaystyle y_{vt}=\sum_{k=0}^{K}x_{kt}\beta_{kv}+\epsilon_{vt}$ (6)

where $v\in\\{1,\cdots,V\\}$ indexes location, i.e. vertices in cs-fMRI,
$t\in\\{1,\cdots,T\\}$ indexes time, and $k\in\\{1,\cdots,K\\}$ indexes a task
or stimulus. $x_{kt}$ is the convolution of a so-called stimulus time course
$s_{k}(\cdot)$ for task $k$ with the canonical hemodynamic response function
$h(\cdot)$; namely

$\displaystyle x_{kt}=\int_{-\infty}^{\infty}h(u)s_{k}(t-u)du\ .$ (7)

$s_{k}(\cdot)$ takes value $1$ when the task is active, and $0$ otherwise. The
canonical hemodynamic response function $h(\cdot)$, visualized in Figure 5,
characterizes the temporal change in oxygenated blood flow for regions of the
brain that are affected by the task. This typically consists of a 2 seconds
delay at the onset of the stimulus, then a dip below baseline followed by a
gradual increase that peaks after 4 seconds, a slow decay to below baseline
level, and a return to baseline level. The duration of each phase may vary
depending on the task, and the total amount of response time is approximately
15 to 20 seconds (Lazar, 2008). We can rewrite the above model for all $V$
vertices and $T$ time points in vector form as

$\displaystyle\boldsymbol{y}=\sum_{k=0}^{K}\boldsymbol{X}_{k}\boldsymbol{\beta}_{k}+\boldsymbol{\epsilon},\quad\boldsymbol{\epsilon}\sim
N(0,\boldsymbol{\Sigma})$ (8)

where $\boldsymbol{y}$ is now a $VT\times 1$ vector created by stacking $V$
time series, each of length $T$, $\boldsymbol{X}_{k}$ is a $VT\times V$ design
matrix containing activation information for task $k$, and
$\boldsymbol{\beta}_{k}$ is a $V\times 1$ vector of activation amplitudes. As
we shall see, the spatial dependence is modeled via non-stationary Matern
priors on the task activation fields $\boldsymbol{\beta}_{k}$. The $VT\times
VT$ error covariance matrix $\boldsymbol{\Sigma}$ is block diagonal, with $V$
$T\times T$ blocks $\\{\boldsymbol{\Sigma}_{v}:v=1,\dots,V\\}$, each capturing
the temporal dependence at a particular vertex.

###  3.1 Modeling varying-range temporal dependence

Temporal dependence can be modeled separately at each vertex, as the
covariance matrix $\boldsymbol{\Sigma}$ in Equation 8 is block diagonal. To
this end, the fMRI time series at a generic vertex $v$ is modeled as a
fractional Gaussian noise process, whose covariance $\boldsymbol{\Sigma}_{v}$
is parameterized by the Hurst parameter $H_{v}$ (see Section 2.1). After
discrete wavelet transformation of the time series at vertex $v$, Equation 6
becomes

$\boldsymbol{y}_{v}^{(w)}=\sum_{k=0}^{K}\boldsymbol{x}^{(w)}_{k}\beta_{kv}+\boldsymbol{\epsilon}_{v}^{(w)},\quad\boldsymbol{\epsilon}_{v}^{(w)}\sim
N(0,\boldsymbol{\Sigma}_{v}^{(w)})$

where $\boldsymbol{x}_{k}^{(w)}$, $\boldsymbol{y}_{v}^{(w)}$ and
$\boldsymbol{\epsilon}_{v}^{(w)}$ are the discrete wavelet transforms of
$x_{kt},y_{vt}$ and $\epsilon_{vt}$, $t=1,\ldots,T$. Because of the
approximately decorrelating property of the discrete wavelet transformation,
the covariance matrix $\boldsymbol{\Sigma}_{v}^{(w)}$ is, to a good
approximation, diagonal. As described in Section 2.1, the diagonal entries
depend on the scale of wavelet transform, and on the range of temporal
dependence at vertex $v$. With the goal of accommodating different ranges of
dependence across brain regions of interest, while keeping the model
parsimonious, we devise a data-driven scheme that groups the ranges into just
$n_{H}$ ($<<V$) distinct levels, thus leaving us with only $n_{H}$ Hurst
parameters to estimate. We proceed as follows:

  1. Step 3.1.1

At each vertex $v$, we obtain preliminary estimates of coefficients
$\\{\hat{\beta}_{kv}:k=1,\ldots,K\\}$ and residuals
$\\{\hat{\epsilon}_{vt}:t=1,\ldots,T\\}$ using linear regression;

  2. Step 3.1.2

At each vertex $v$, we take the discrete wavelet transform of the residuals
$\\{\hat{\epsilon}_{vt}:t=1,\ldots,T\\}$, producing detail coefficients
$\\{d_{jk}:j=1,\dots,J,k=1,\ldots,n/2^{j}\\}$ and approximation coefficients
$\\{a_{Jk}:k=1,\ldots,n/2^{J}\\}$, where $J$ is the coarsest scale of
transformation;

  3. Step 3.1.3

At each vertex $v$, we estimate the variances $\\{S_{d_{j}}:j=1,\dots,J\\}$ in
Equation 2 computing the variances of $\\{d_{jk}:j=1,\ldots,n/2^{j}\\}$;

  4. Step 3.1.4

At each vertex $v$, we estimate $\gamma$ via the linear relationship
$\log_{2}(S_{d_{j}})=\gamma j+\mathcal{O}(1)$ and produce a preliminary
estimate of the Hurst parameter as $\hat{H}_{v}=(\hat{\gamma}+1)/2$;

  5. Step 3.1.5

For each region of interest $r=1,\ldots,R$, we obtain the median of such such
estimates across vertices in the region,
$M_{r}^{(H)}=\text{med}\\{\hat{H}_{v}:v\in r\\}$;

  6. Step 3.1.6

We group the $R$ regions into $n_{H}$ clusters based on their median estimates
$\\{M_{r}^{(H)}:r=1,\ldots,R\\}$, producing cluster memberships
$C_{r}\in\\{1,\ldots,n_{H}\\}$ for each region $r$. Since no two regions share
a vertex, cluster memberships $C_{v}\in\\{1,\ldots,n_{H}\\}$ can also be
attributed to each individual vertex $v\in\\{1,\ldots,V\\}$ based on the
membership of the regions it belongs to;

  7. Step 3.1.7

At each vertex $v$, we assume the time series to be distributed as a
fractional Gaussian process parameterized by the cluster-specific Hurst
parameter corresponding to the vertex membership.

At the end of this procedure, any two vertices $v$ and $v^{\prime}$ such that
$C_{v}=C_{v^{\prime}}$ will have the same Hurst parameter, and thus their
$\boldsymbol{\Sigma}_{v}^{(w)}$ and $\boldsymbol{\Sigma}_{v^{\prime}}^{(w)}$
will have the same diagonal entries expressed in Equation 2. Note that we are
not using the procedure to estimate the parameters, we are just employing
rough, preliminary estimates $\hat{H}_{v}$ to group regions and vertices, and
thus reduce the number of parameters to be estimated. Note also that, in Step
3.1.2, the number of available coefficients to estimate $S_{d_{j}}$ is
$n/2^{j}$. As $j$ increases, this number decreases by a factor of 2, causing
estimation of $S_{d_{j}}$ to be volatile at large scales. Hence, we recommend
using only scales with at least 16 detail coefficients to generate the
preliminary $\hat{H}_{v}$’s in Step 3.1.3. Finally, we implement the
clustering in Step 3.1.6 using a simple K-mean algorithm with the number of
clusters $n_{H}$ fixed beforehand (more clusters correspond to more distinct
dependence ranges, and thus to a less parsimonious model, with more parameters
to be estimated). While more sophisticated approaches could be employed, in
our experience a K-means with $n_{H}\leq 5$ works sufficiently well (larger
$n_{H}$ values do not provide any added improvements in model fit). The choice
to perform clustering on regions instead of directly on vertices is based on
the assumption that vertices in the same region have similar temporal
dependence. This is supported by Park et al. (2010), in which the authors
estimated Hurst parameters separately for each brain voxel, after removing
spatial covariance from the data; even with spatial covariance removed, Hurst
parameter estimates still mapped nicely into obvious brain structures.
Clustering regions has the added advantage of reducing spurious temporal
noise, which may arise naturally in the data, and/or from previous estimation
steps in the algorithm. We also note that, while beyond the scope of the
current article, size and interpretation of the regions being clustered can be
flexibly altered by changing the choice of parcellation.

Fig. 4: Boxplots of $\hat{H}_{v}$ values grouped by $50$ left hemisphere
cortical regions from the Schaefer parcellation (Schaefer et al., 2018). The
numbering of the regions does not correspond to any particular arrangement or
order. The histogram on the secondary y-axis (right) shows the overall
distribution. The horizontal red line marks $H=0.5$, which corresponds to a
white noise process.

Figure 4 provides an illustration of the preliminary estimates $\hat{H}_{v}$
which are represented as boxplots for each region of interests in the Schaefer
100 parcellation (see Figure 1, each hemisphere’s cortical surface is
partitioned into 50 regions; Schaefer et al., 2018). Specifically, the
illustration uses resting state time series ($1,024$ time points) for $6,000$
vertices in the left hemisphere of an individual from the healthy young adult
data set of the Human Connectome Project (HCP; Van Essen et al., 2013). Of
course, in task-related (as opposed to resting state) fMRI data, temporal
dependence, and thus Hurst parameter estimates, may be influenced by the task.
The horizontal red line in Figure 4 marks $H=0.5$, which corresponds to a
white noise process. Although many regions have similar $\hat{H}_{v}$ values
(e.g., regions 32–39), we do see evidence of varying ranges of dependence
(values between $H=0.5$ and $H=1$) even in this resting state data. The fact
that for some regions almost all coarse estimates are above $H=0.5$ (e.g.,
regions 3, 8, 12, 20, 43 and 47), together with the skew of the overall
distribution (histogram on the right of Figure 4), strongly suggest the
existence of long-range dependence.

###  3.2 Modeling non-stationary spatial dependence

In this section, the task indices are suppressed for convenience, so that
$\boldsymbol{\beta}_{k}=\boldsymbol{\beta}$ and $\beta_{kv}=\beta_{v}$. The
SPDE framework can be extended to accommodate non-stationary processes by
introducing location-dependent SPDE parameters. Following Lindgren et al.
(2011), a parameter $\tau$ is used to scale $\beta$ in the SPDE in (3), while
keeping the variance of the Gaussian noise constant; that is

$\displaystyle(\kappa^{2}(\boldsymbol{s})-\Delta)^{\frac{\nu+d/2}{2}}\tau(\boldsymbol{s})\beta(\boldsymbol{s})=\mathcal{W}(\boldsymbol{s}),\quad\boldsymbol{s}\in\mathbb{R}^{d}$

where $\boldsymbol{s}$ indicates the location in a space of generic dimension
$d$. Note the distinction between this continuous coordinates notation, which
is necessary for distance calculations, and the vertex indices
$v\in\\{1,\ldots,V\\}$. The marginal variance becomes

$\displaystyle\sigma^{2}(\boldsymbol{s})=\frac{\Gamma(\nu)}{\Gamma(\nu+d/2)(4\pi)^{d/2}\kappa^{2\nu}(\boldsymbol{s})\tau(\boldsymbol{s})}\
.$ (9)

As noted in Lindgren & Rue (2015), it is often more intuitive to
reparameterize the SPDE using the standard deviation $\sigma$ and the range
$\rho=(8\nu)^{0.5}/\kappa$, where $\rho$ is the distance at which 2
observations are approximately independent (correlation approximately $0.13$)
for all $\nu>0.5$. We model $\sigma(\boldsymbol{s})$ and
$\rho(\boldsymbol{s})$ as functions of the location and translate them back to
$\tau(\boldsymbol{s})$ and $\kappa(\boldsymbol{s})$. Specifically, we set

$\displaystyle\begin{split}&\log(\sigma(\boldsymbol{s}))=\log(\sigma_{0})+\theta_{1}\delta(\boldsymbol{s})\\\
&\log(\rho(\boldsymbol{s}))=\log\rho=\log(\rho_{0})+\theta_{2}\end{split}$
(10)

where $\sigma_{0}$ and $\rho_{0}$ are baseline standard deviation and range,
and $\delta(\boldsymbol{s})$ is some chosen local variability score. Using the
relationship $\rho=(8\nu)^{0.5}/\kappa$ and Equation 9 we thus get

$\displaystyle\log\kappa(\boldsymbol{s})$
$\displaystyle=\log\kappa=\frac{\log(8\nu)}{2}-\log(\rho_{0})-\theta_{2}$
$\displaystyle=\log\kappa_{0}-\theta_{2}$
$\displaystyle\log\tau(\boldsymbol{s})$
$\displaystyle=\frac{1}{2}\log\left(\frac{\Gamma(\nu)}{\Gamma(\nu+\frac{d}{2})(4\pi)^{\frac{d}{2}}}\right)-\log(\sigma_{0})-\nu\left(\frac{\log(8\nu)}{2}-\log(\rho_{0})\right)-\theta_{1}\log(\delta(\boldsymbol{s}))+\theta_{2}\nu$
$\displaystyle=\log\tau_{0}-\theta_{1}\delta(\boldsymbol{s})+\theta_{2}\nu\ .$

Let $\boldsymbol{T}=\text{diag}(\tau(\boldsymbol{s_{v}}):v=1,\ldots,V)$, then
the precision matrix in Equation 5 becomes
$\boldsymbol{T}\boldsymbol{Q}_{\alpha}\boldsymbol{T}$. Lindgren et al. (2011);
Lindgren & Rue (2015) remarked that the link between SPDE and Matern
parameters is no longer valid in the non-stationary adaptation. However, when
$\kappa(\boldsymbol{s})$ and $\tau(\boldsymbol{s})$ vary slowly over the
domain, Equations LABEL:model_sigma_rho provide a valid approximation to the
local variances and correlation ranges. This requirement is achieved if the
local variability score $\delta(\boldsymbol{s})$ is also slowly varying with
$\boldsymbol{s}$. For our applications, we calculate the local variability
score as follows:

  1. Step 3.2.1

At each vertex $v$, we obtain initial estimates $\hat{\beta}_{v}$ using linear
regression;

  2. Step 3.2.2

At each vertex $v$, we calculate $\delta_{v}=\delta(\boldsymbol{s}_{v})$ as
the standard deviation of the values $\hat{\beta}_{i}$, $i\in N_{v}$, where
$N_{v}$ is a neighborhood of $v$.

The size of the neighborhood $N_{v}$ in Step 3.2.2 can be chosen such that
$\delta(\boldsymbol{s})$ is smooth in $\boldsymbol{s}$. For sufficiently dense
data such as ours, nearest neighbors suffice.

###  3.3 Tying it all together: joint posterior

Given the wavelet transformed data $\boldsymbol{y}^{(w)}$, we can now write
the joint posterior of marginal temporal variance $\sigma$, temporal Hurst
parameters $H_{1},\ldots,H_{n_{H}}$ (recall these are reduced in number with
respect to the number of vertices and are not task dependent) and, for each
task $k=1,\dots,K$, activation fields
$\boldsymbol{\beta}_{k}=\\{\beta_{k1},\ldots,\beta_{kV}\\}$ and spatial
hyperparameters $\boldsymbol{\theta}_{k}=\\{\theta_{k1},\theta_{k2}\\}$. Such
joint posterior is

$\displaystyle\pi(\boldsymbol{\beta}_{1},\ldots,\boldsymbol{\beta}_{K},\boldsymbol{\theta}_{1},\ldots,\boldsymbol{\theta}_{K},H_{1},\ldots,H_{n_{H}},\sigma|\boldsymbol{y}^{(w)})$
$\displaystyle\propto\pi(\sigma)\left[\prod_{k=1}^{K}\pi(\boldsymbol{\theta}_{k})\pi(\boldsymbol{\beta}_{k}|\boldsymbol{\theta}_{k})\right]\left[\prod_{i=1}^{n_{H}}\pi(H_{i})\right]\left[\prod_{v=1}^{V}\pi(\boldsymbol{y}_{v}^{(w)}|\boldsymbol{\beta}_{1},\ldots,\boldsymbol{\beta}_{K},H_{1},\ldots,H_{n_{H}},\sigma)\right]$
$\displaystyle\propto\pi(\sigma)\left[\prod_{k=1}^{K}\pi(\boldsymbol{\theta}_{k})|\boldsymbol{Q}(\boldsymbol{\theta}_{k})|^{1/2}\exp\left(-\frac{1}{2}\boldsymbol{\beta}_{k}^{\top}\boldsymbol{Q}(\boldsymbol{\theta}_{k})\boldsymbol{\beta}_{k}\right)\right]\left[\prod_{i=1}^{n_{H}}\pi(H_{i})\right]\times$
$\displaystyle\hskip
34.14322pt\left[\prod_{v=1}^{V}|\boldsymbol{\Sigma}_{v}^{(w)}|^{-1/2}\exp\left(\frac{1}{2}\bigg{(}\boldsymbol{y}_{v}^{(w)}-\sum_{k=0}^{K}\boldsymbol{x}^{(w)}_{k}\beta_{kv}\bigg{)}^{\top}{\boldsymbol{\Sigma}_{v}^{(w)}}^{-1}\bigg{(}\boldsymbol{y}_{v}^{(w)}-\sum_{k=0}^{K}\boldsymbol{x}^{(w)}_{k}\beta_{kv}\bigg{)}\right)\right]$

where $\pi(\sigma)$ is Gamma(1,1); $\\{\pi(H_{i}):i=1,\ldots,n_{H}\\}$ are
Uniform(0,1); $\\{\pi(\theta_{kj}):k=1,\ldots,K;\ j=1,2\\}$ are Gaussian with
mean $0$ and precision $0.3$;
$\\{\pi(\boldsymbol{\beta}_{k}|\boldsymbol{\theta}_{k}):k=1,\ldots,K\\}$ are
GMRFs with mean $0$ and precision
$\boldsymbol{Q}(\boldsymbol{\theta})=\boldsymbol{T}\boldsymbol{Q}_{\alpha}\boldsymbol{T}$,
as described in Section 3.2. Based on this joint posterior distribution,
marginal posteriors for different parameters are approximated using INLA, as
described in Section 2.3. We shall call our approach _non-stationary, varying-
range Bayesian General Linear Model_ , hereinafter _NSVR-BayesGLM_.

##  4 Results

###  4.1 Simulations

Since the main target of our proposed approach are cs-fMRI data, which provide
a rendition of the cerebral cortex via a 2D surface mesh, it makes sense to
consider 2D simulations. Here, to facilitate comparisons with prior
literature, we build simulations on a traditional 2D brain slice.
Specifically, we follow the set-up of Mejia et al. (2020) and compare results
between our _NSVR-BayesGLM_ and their _BayesGLM_ approach. A $46\times 55$
image is constructed from a brain mask provided by the SPM8 software.
Activation is placed at four different sites, with varying signal strengths,
spreads and degrees of smoothness. Signal strength is defined as the ratio of
activation magnitude at the center of the site, $M=|\beta_{center}|$, to noise
standard deviation, $\sigma$. Spread is defined as the radius $r$ of the site.
Smoothness is defined by how fast the signal decays from the center to the
edge of the site, via the exponential function $e^{-\lambda d}$, where $d$ is
the distance to the center and $\lambda$ the parameter that controls
smoothness. Thus, for an activation site with signal strength $M/\sigma$,
radius $r$ and smoothness parameter $\lambda$, the activation magnitude at a
location $d$ distance away from the center is $Me^{-\lambda d}$ if $d<r$ and 0
otherwise. Then, at different sites, we add error time series with different
dependence ranges. Specifically, we place fractional Gaussian noise processes
with $H=0.8$ in two active regions to represent long temporal dependence,
$H=0.4$ in the remain two active regions to represent short temporal
dependence, and $H=0.5$ outside of active regions to represent white noise. We
assume no overlapping between activation sites. For convenience, we choose
$\sigma=1$ and, since the discrete wavelet transform requires $T$ to be a
power of two, $T=2^{9}=512$. In reality, for any $T$, one can perform some
“padding” to meet this requirement. A summary of the simulation design can be
found in Table 2 – the two tasks, which are depicted in Figure 5, alternate in
a single-block design scheme.

| Signal strength | Activation center | Smoothness | H
---|---|---|---|---
|  | top | 0.2 | 0.8
|  | right | 0.2 | 0.4
|  | bottom | 0.05 | 0.4
Task 1 | 2 | left | 0.05 | 0.8
|  | top | 0.2 | 0.8
|  | right | 0.2 | 0.4
|  | bottom | 0.05 | 0.4
Task 2 | 3 | left | 0.05 | 0.8
Table 2: Details of simulation design. Fig. 5: Stimulus function (left),
canonical hemodynamic response function (middle), and the design matrix that
results from their convolution (right).

Both _NSVR-BayesGLM_ and _BayesGLM_ use SPDE priors to model spatial
dependence of the activation fields, and both are implemented with the R
package R-INLA using the parallel PARDISO solver (Rue et al., 2023). The
varying temporal dependence component of _NSVR-BayesGLM_ is written and
implemented in C as a custom add-on that operates within R-INLA. While
_BayesGLM_ prewhitens all time-series with an $AR(6)$ model and assumes
stationary spatial priors, _NSVR-BayesGLM_ allows varying ranges of temporal
dependence in wavelet space and non-stationary spatial priors in one, unifying
framework.

Figure 6 shows the true activation magnitudes (left) and estimates by _NSVR-
BayesGLM_ (middle) and _BayesGLM_ (right), for all simulation settings.
Different gradients of colors at activation sites portray different degrees of
spatial smoothness. Whereas the non-stationary feature of _NSVR-BayesGLM_
allows it to capture this varying spatial smoothness, _BayesGLM_ ’s stationary
priors assume the same smoothness level across all locations. This is evident
in the latter’s tendency to oversmooth activation magnitudes in the bottom and
the left activation sites. _BayesGLM_ also underestimates activation
magnitudes, which is most evident at the centers of activation sites. This is
likely a combined effect of prewhitening with $AR(6)$ models, as discussed in
Section 1.1, and of assuming stationarity. In addition to the activation
magnitudes captured in Figure 6, each activation site has a different temporal
range of dependence. Here we use $n_{H}=3$; _NSVR-BayesGLM_ automatically
groups regions with similar temporal dependence ranges – such as the right and
bottom activation sites – and produces Hurst parameter estimates
$\hat{H}=0.373,0.825,0.486$ for the true $H=0.4,0.8,0.5$, respectively. Note
that, if instead we use $n_{H}=5$, _NSVR-BayesGLM_ produces estimates
$\hat{H}=0.390,0.396,0.853,0.871,0.487$ for the true $H=0.4,0.4,0.8,0.8,0.5$,
respectively. _BayesGLM_ does not accommodate varying ranges of dependence and
therefore does not produce estimates of such ranges.

Fig. 6: Simulated data (left) and activation magnitude estimates by NSVR-
BayesGLM (middle) and BayesGLM (right). Top and bottom plots concern Tasks 1
and 2, respectively. The ventricles, which house cerebral spinal fluid and are
located near the center of the brain, are plotted in white and are not
analyzed. Fig. 7: Simulated data (left) and activation region estimates by
_NSVR-BayesGLM_ (middle) and _BayesGLM_ (right). Top and bottom plots concern
Tasks 1 and 2, respectively. False positives are plotted in yellow and false
negatives in red (as in Figure 6, ventricles are plotted in white and are not
analyzed).

Figure 7 shows the true activation sites (left), and estimates by _NSVR-
BayesGLM_ (middle) and _BayesGLM_ (right), for all simulation settings. In
both approaches, an estimated activation site is determined as the largest set
$D$ of contiguous locations such that the joint posterior probability
$P(\beta_{v}>0\text{ for all }v\in D\,|\,\boldsymbol{y})>1-\alpha$, where
$\alpha$ is a posterior probability threshold which we fixed at $0.05$. This
way of estimating activation sites is called the excursions set approach; it
was introduced by Bolin & Lindgren (2015) and is implemented in the R package
excursions (Bolin & Lindgren, 2018). As detailed in Mejia et al. (2020), since
the joint posterior probability is over all vertices, there is no need to
adjust for multiple comparisons – and the approach appears to be fairly
insensitive to the choice of $\alpha$. Based on the top plots of Figure 7, in
Task 1, both _NSVR-BayesGLM_ and _BayesGLM_ produce false negatives around the
edge of the top and right activation sites, likely an effect of low signal
strength in those areas. Notably, _BayesGLM_ is also particularly prone to
false positives, even in the case of strong signals (Task 2, bottom plots).
The stationarity assumption of _BayesGLM_ acts like a smoother of activation
magnitudes, dampening those that are larger and boosting those that are
smaller. This has the effect of making false positives more likely around the
edges of the bottom and left activation sites, while making false negatives
less likely around the edges of the top and right activation sites.

###  4.2 An application to visual working memory

In this application, we use cs-fMRI data to characterize activation and
temporal dependence in a visual working memory task. The data is part of a
working memory study included in the healthy young adult data set of the Human
Connectome Project (HCP; Van Essen et al., 2013). Preprocessing is conducted
with the HCP minimal preprocessing pipeline (Glasser et al., 2013), and
individual subjects’ cortical areas are registered to a common template for
the purpose of group analysis and intersubject comparisons (Robinson et al.,
2014). In the study, subjects are presented with blocks of trials consisting
of different images. Subjects indicate if an image is a 2-back repeat (i.e.,
occurred 2 trials before) in a ‘2-back’ condition, or if an image matches a
cue stimulus in a ‘0-back’ condition (as a control for working memory). In
summary, four ‘2-back’ blocks and four ‘0-back’ blocks alternate randomly, for
a total of eight blocks per run, each block consisting of 10 trials of 2.5
seconds each. For all subjects, a working memory score is also measured, using
the NIH Toolbox List Sorting Working Memory Test (Tulsky et al., 2014). For
the purposes of our analysis, we select 10 subjects from the 100 unrelated
subjects (no family relations) available from the HCP; namely, the 5 highest-
and the 5 lowest-ranking in the List Sorting Working Memory Test.

We analyze left and right hemispheres separately, using the approach described
in Section 3. To reduce computational burden, we use the Connectome Workbench
to subsample vertices in each hemisphere from $32,492$ to $\approx 3,000$.
Thus, for our modeling exercise $V\approx 3,000$ and $T=401$. The BOLD
response and design matrix are centered and scaled at each vertex for
numerical stability. In addition, spurious variability caused by subject
movements and scanner drift is removed from the BOLD signal, using $6$ motion
parameters estimates (translation and rotation parameters, each in 3
directions) and their temporal derivatives as regressors. We carry out the
modelling on the ‘very-inflated’ surface, and visualize results on the
‘inflated’ surface (see Figure 1). As in Section 4.1, _NSVR-BayesGLM_ is
implemented with the R package R-INLA, using the parallel PARDISO solver (Rue
et al., 2023), and with the varying temporal dependence component written in C
as a custom add-on within R-INLA. Manipulation of cortical surface fMRI data
in R is made possible through the ciftiTools package (Pham et al., 2022).
Using 8 parallel threads, computation takes on average 30 minutes per
hemisphere per subject, on a Macbook Pro with M2 Max chip and 96GB memory.

Fig. 8: Activation field estimates for the ‘0-back’ (left) and the ‘2-back’
(right) tasks, averaged across sessions and subjects within the high (top) and
low (bottom) working memory scoring groups.

Figure 8 shows the activation field estimates for the ‘0-back’ task (left) and
‘2-back’ task (right), taken as averages across sessions and subjects
belonging to the same group (high working memory score, top, and low working
memory score, bottom). Because _NSVR-BayesGLM_ allows non-stationarity, we can
capture activation patterns with noticeably different degrees of smoothness
across regions. In both tasks, activation is, as expected, most prominent in
the striate and extrastriate visual cortex (BA V1, BA V2), followed by
subregions of several functional divisions such as precentral ventral
attention (BA 6, BA 44), salience/ventral attention (BA 44, BA 6),
sommatosensory motor (BA 1, BA 2, BA 4p), prefrontal cortex (BA 45),
intraparietal sulcus (IPS), and superior parietal lobule (SPL). This is
consistent with the conclusion in a recent meta analysis of visual working
memory task (Li et al., 2022). Notably, the IPS and SPL have been documented
to play a crucial role in visuo-spatial attention and the analysis of spatial
elements (see, e.g., Papadopoulos et al., 2018; MacSweeney et al., 2008).
Functional regions that contain a larger, and/or more significantly activated
areas during the ‘2-back’ task include the precentral dorsal attention,
salience/ventral attention, prefrontal cortex, inferior parietal lobule and
precuneus. While the prefrontal cortex is key in regulating a large number of
higher-order executive functions, including working memory, the precuneus is
well-known for its role in visuo-spatial imagery and episodic memory retrieval
(see, e.g., Cavanna & Trimble, 2006). When contrasting high- and low-scoring
groups, the same regions seem to be activated in each task, but interestingly
activation magnitudes are higher across the board in the low scoring group.
This may suggest an increased effort in performing the tasks from subjects who
score lower in the working memory test.

Fig. 9: Hurst parameter estimates, averaged across subjects within the low
(left) and high (right) working memory scoring groups. No regions exhibit
short-range dependence in either group, and high-scoring subjects exhibit
longer range dependences than low-scoring subjects.

Figure 9 shows Hurst parameter estimates for different brain regions, taken
again as averages across subjects belonging to the same group. We note that
using $n_{H}=5$, all regions report Hurst estimates larger than $0.5$, which
strengthen the case for an approach that, like _NSVR-BayesGLM_ , can
accommodate long range dependences. Interestingly, subjects who score high in
the working memory test have noticeably higher Hurst estimates across almost
all brain regions. This suggest an intriguing association between temporal
dependence ranges and working memory performance – and one that our region-
specific Hurst estimates, unlike the single score value of the test, allow us
to articulate spatially. While this is only a preliminary result, it may point
towards a promising, dependence range-based approach to the study of memory-
related diseases. It also further strengthens our recommendation against
prewhitening, which erases potentially valuable information contained in model
residuals.

##  5 Discussion

In this article we propose _NSVR-BayesGLM_ , a first-of-its-kind approach that
fully integrates the modelling of non-stationary spatial dependence and
varying ranges of temporal dependence in cs-fMRI data. The former is
implemented through an extension of the spatial SPDE prior approach in Mejia
et al. (2020), which allows for spatial non-stationarity driven by local
spatial features of the data and results in a sparsely parameterized model.
Different ranges of temporal dependence are also modeled in a data-driven
fashion, combining the use of several fractional Gaussion processes of varying
persistence and wavelet transformations. We demonstrate via simulations that
_NSVR-BayesGLM_ can better accommodate local activation features compared to a
stationary counterpart, while retaining comparable power and improving false
positive control. Notably, in addition to assuming stationarity, the competing
method employs autoregressive models of order $p=6$ to prewhiten the data,
while _NSVR-BayesGLM_ accurately estimates and incorporates different ranges
of temporal dependence at different locations. Finally, we apply _NSVR-
BayesGLM_ to a visual working memory cs-fMRI data set. The activation patterns
are found to be consistent with existing literature. Moreover, _NSVR-BayesGLM_
allows us to unveil different degrees of smoothness in the signals
characterizing different regions, and an intriguing association between
estimated temporal dependence profiles and working memory performance –
something that may point towards further applications in the study of memory-
related diseases.

A limitation of the current proposal is that, while advancing the modeling of
spatial and temporal dependence with respect to existing approaches, it still
treats them as separate. Spatial dependence is embedded within the activation
coefficients’ prior hyper-parameters, while temporal dependence is imposed
upon the error terms and is reflected directly in the likelihood. In reality,
it is likely that the spatial correlations characterizing activation fields
change with time, e.g., as task stimuli switch on and off. Ideally, one should
explicitly incorporate interactions between spatial and temporal dependence –
but pursuing such an extension would still be prohibitive in terms of
computational burden.

Indeed, computational cost remains a major challenge that influences modeling
choices. As reported in Section 4.2, the current model estimation has a
runtime of approximately 30 minutes per hemisphere, compared to _BayesGLM_ ’s
45 seconds. The main burden comes from inefficient memory allocation for large
spatio-temporal models; R-INLA takes up nearly 85GB of RAM during the
estimation stage for a data set of $\approx 3000$ locations/vertices and
$\approx 500$ time points. Note that, if we use prewhitening and do not
explicitly model temporal dependence (as in the case of _BayesGLM_), the non-
stationary spatial model is just as fast and memory-efficient. For the purpose
of circumventing memory issue and reducing computational burden, we could
consider an approach that combines INLA with the MCMC algorithm. This works by
iterating between drawing MCMC samples from the posterior marginals of
temporal dependence hyperparameters, and fitting only the non-stationary
spatial component with INLA, conditioning on the MCMC samples. This “fusion”
approach, also known as the Metropolis-Hastings with INLA, was proposed by
Gomez-Rubio & Rue (2018) and shown to have a performance comparable to INLA
via numerical examples and heuristic arguments, for some selected models. For
our model, potential complications with parameter support boundaries and
approximation errors of the conditional marginal likelihood might exist, and
will need to be studied in future work. Complete code implementing the current
version of _NSVR-BayesGLM_ as well as data for reproducing our simulations are
available at https://github.com/hqd1/SPfmri.

## Acknowledgments

We utilized data from the Human Connectome Project, WU-Minn Consortium
(Principal Investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657)
funded by the 16 NIH Institutes and Centers that support the NIH Blueprint for
Neuroscience Research, and by the McDonnell Center for Systems Neuroscience at
Washington University. We utilized substantial portions of the R code
developed in Mejia et al. (2020), we thank the authors for making such code
available and for providing advice. The work of M.A. Cremona was partially
supported by the Natural Sciences and Engineering Research Council of Canada
(NSERC), by the Fonds de recherche du Québec Health (FRQS), and by the Faculty
of Business Administration, Université Laval. The work of F. Chiaromonte was
partially supported by the Huck Institute of the Life Sciences of the
Pennsylvania State University.
M.A. Cremona is also affiliated with CHU de Québec – Université Laval Research
Center, Canada. F. Chiaromonte is also affiliated with the Sant’Anna School of
Advanced Studies, Italy.

## References

  * (1)
  * Bolin & Lindgren (2015) Bolin, D. & Lindgren, F. (2015), ‘Excursion and Contour Uncertainty Regions for Latent Gaussian Models’, Journal of the Royal Statistical Society, Series B 77, 85–106.
  * Bolin & Lindgren (2018) Bolin, D. & Lindgren, F. (2018), ‘Calculating Probabilistic Excursion Sets and Related Quantities Using Excursions’, Journal of Statistical Software 86, 1–20.
  * Bowman (2005) Bowman, F. (2005), ‘Spatio-temporal Modeling of Brain Activity’, Biostatistics 6, 558–575.
  * Bullmore et al. (2003) Bullmore, E., Fadili, J., Breakspear, M., Salvadore, R., Suckling, J. & Brammer, M. (2003), ‘Wavelets and Statistical Analysis of Functional Magnetic Resonance Images of the Human Brain’, Statistical Methods in Medical Research 12, 375–399.
  * Castruccio et al. (2018) Castruccio, S., Ombao, H. & Genton, M. G. (2018), ‘A Scalable Multi-Resolution Spatio-Temporal Model for Brain Activation and Connectivity in fMRI Data’, Biometrics 74, 823–833.
  * Cavanna & Trimble (2006) Cavanna, A. E. & Trimble, M. R. (2006), ‘The Precuneus: A Review of Its Functional Anatomy and Behavioural Correlates’, Brain 129(3), 564–583.
  * Derado et al. (2010) Derado, G., Bowman, F. D. & Kilts, C. D. (2010), ‘Modeling the Spatial and Temporal Dependence in fMRI Data’, Biometrics 66, 949–957.
  * Fadili & Bullmore (2002) Fadili, M. & Bullmore, E. (2002), ‘Wavelet-generalized least squares: A new BLU estimator of linear regression models with 1/f errors’, NeuroImage 15, 217–232.
  * Flandrin (1992) Flandrin, P. (1992), ‘Wavelet Analysis and Synthesis of Fractional Brownian Motion’, IEEE Transactions on Information Theory 38, 910–917.
  * Glasser et al. (2016) Glasser, M. F., Smith, S. M., Marcus, D. S., Andersson, J. L. R., Auerbach, E. J., Behrens, T. E. J., Coalson, T. S., Harms, M. P., Jenkinson, M., Moeller, S., Robinson, E. C., Sotiropoulos, S. N., Xu, J., Yacoub, E., Ugurbil, K. & Van Essen, D. C. (2016), ‘The human connectome project’s neuroimaging approach’, Nature Neuroscience 19(9), 1175–1187.
  * Glasser et al. (2013) Glasser, M. F., Sotiropoulos, S. N., Wilson, J. A., Coalson, T. S., Fischl, B., Andersson, J. L., J., X., Jbabdi, S., Webster, M., Polimeni, J. R., Van Essen, D. C. & Jenkinson, M. (2013), ‘The Minimal Preprocessing Pipelines for the Human Connectome Project’, NeuroImage 80, 105–124.
  * Gomez-Rubio & Rue (2018) Gomez-Rubio, V. & Rue, H. (2018), ‘Markov chain monte carlo with the integrated nested laplace approximation’, Statistics and Computing 28, 1033–1051.
  * Kang et al. (2013) Kang, H., Yang, X., Bryan, F. W., Tripp, C. M. & Landman, B. A. (2013), ‘Whole brain functional connectivity using multi-scale spatio-spectral random effects model’, Multimodal Brain Image Analysis. Lecture Notes in Computer Science 8159.
  * Lazar (2008) Lazar, N. (2008), The statistical analysis of functional MRI data, Springer.
  * Li et al. (2022) Li, X., O’Sullivan, M. J. & Mattingley, J. B. (2022), ‘Delay activity during visual working memory: A meta analysis of 30 fmri experiments’, NeuroImage 255, 119204.
  * Li & Huang (2014) Li, Y. C. & Huang, Y. A. (2014), ‘Fractal analysis of spontaneous fluctuations of the bold signal in the human brain networks’, Journal of Magnetic Resonance Imaging 39, 1118–1125.
  * Lindgren & Rue (2015) Lindgren, F. & Rue, H. (2015), ‘Bayesian spatial modelling with r-inla’, Journal of Statistical Software 63(19).
  * Lindgren et al. (2011) Lindgren, F., Rue, H. & Lindström, J. (2011), ‘An explicit link between gaussian fields and gaussian markov random fields: The stochastic partial differential equation approach’, Journal of the Royal Statistical Society, Series B 73, 423–498.
  * MacSweeney et al. (2008) MacSweeney, M., Capek, C. M., Campbell, R. & Woll, B. (2008), ‘The signing brain: the neurobiology of sign language’, Trends in Cognitive Sciences 12(11), 432–440.
https://www.sciencedirect.com/science/article/pii/S1364661308002192

  * Mallat (2009) Mallat, S. (2009), A wavelet tour of signal processing: the sparse way, Elsevier /Academic Press.
  * Mejia et al. (2020) Mejia, A. F., Yue, Y., Bolin, D., Lindgren, F. & Lindquist, M. (2020), ‘A bayesian general linear modeling approach to cortical surface fmri data analysis’, Journal of the American Statistical Association 115(530), 501–520.
  * Papadopoulos et al. (2018) Papadopoulos, A., Sforazzini, F., Egan, G. & Jamadar, S. (2018), ‘Functional subdivisions within the human intraparietal sulcus are involved in visuospatial transformation in a non‐context‐dependent manner’, Human Brain Mapping 39(1), 354–368.
  * Park et al. (2010) Park, C., Lazar, N. A., Ahn, J. & Sornborger, A. (2010), ‘A multiscale analysis of the temporal characteristics of resting-state fmri data’, Journal of Neoroscience Methods 193, 334–342.
  * Pham et al. (2022) Pham, D. D., Muschelli, J. & Mejia, A. F. (2022), ‘ciftitools: A package for reading, writing, visualizing, and manipulating cifti files in r’, NeuroImage 250, 118877.
  * Risk et al. (2016) Risk, B. B., Matteson, D. S., Spreng, R. N. & Ruppert, D. (2016), ‘Spatiotemporal mixed modeling of multi-subject task fmri via method of moments’, NeuroImage 142, 280–292.
  * Robinson et al. (2014) Robinson, E., Jbabdi, S., Glasser, M., Andersson, J., Burgess, G., Harms, M., Smith, S., Van Essen, D. & Jenkinson, M. (2014), ‘Msm: a new flexible framework for multimodal surface matching’, NeuroImage 100, 414–426.
  * Rue & Martino (2009) Rue, H. & Martino, S. (2009), ‘Approximate bayesian inference for latent gaussian models by using integrated nested laplace approximations’, Journal of the Royal Statistical Society, Series B 71, 319–392.
  * Rue et al. (2023) Rue, H., Martino, S., Lindgren, F., Simpson, D. & Riebler, A. (2023), ‘Code repository for r-inla: Approximate bayesian inference using integrated nested laplace approximations’.
https://www.r-inla.org/download-install

  * Schaefer et al. (2018) Schaefer, A., Kong, R., Gordon, E. M., Laumann, T. O., Zuo, X.-N., Holmes, A. J., Eickhoff, S. B. & Yeo, B. T. T. (2018), ‘Local-global parcellation of the human cerebral cortex from intrinsic functional connectivity mri’, Cerebral Cortex 28, 3095–3114.
  * Tulsky et al. (2014) Tulsky, D. S., Carlozzi, N., Chiaravalloti, N. D., Beaumont, J. L., Kisala, P. A., Mungas, D., Conway, K. & Gershon, R. (2014), ‘Nih toolbox cognition battery (nihtb-cb): The list sorting test to measure working memory’, Journal of the International Neuropsychological Society 20(6), 599–610.
  * Van Essen et al. (2013) Van Essen, D. C., Smith, S. M., Barch, D. M., Behrens, T. E. J., Yacoub, E. & Ugurbil, K. (2013), ‘The wu-minn human connectome project: An overview’, NeuroImage 80, 62–79.
  * Whittle (1954) Whittle, P. (1954), ‘On stationary processes in the plane’, Biometrika 41, 434–449.
  * Whittle (1963) Whittle, P. (1963), ‘Stochastic processes in several dimensions’, Bulletin of the International Statistical Institute 40, 974–994.