A Bayesian General Linear Modeling Approach to Cortical Surface fMRI Data Analysis

TL;DR

This paper introduces a Bayesian general linear model for cortical surface fMRI data analysis, offering improved spatial modeling, activation detection, and multi-subject analysis over classical methods.

Contribution

It presents the first Bayesian spatial model for cs-fMRI, utilizing INLA for efficient computation and an excursions set method for activation detection without multiple comparisons correction.

Findings

Smoother activation estimates with the Bayesian model.

More accurate false positive control.

Higher power to detect true activations.

Abstract

Cortical surface fMRI (cs-fMRI) has recently grown in popularity versus traditional volumetric fMRI, as it allows for more meaningful spatial smoothing and is more compatible with the common assumptions of isotropy and stationarity in Bayesian spatial models. However, as no Bayesian spatial model has been proposed for cs-fMRI data, most analyses continue to employ the classical, voxel-wise general linear model (GLM) (Worsley and Friston 1995). Here, we propose a Bayesian GLM for cs-fMRI, which employs a class of sophisticated spatial processes to flexibly model latent activation fields. We use integrated nested Laplacian approximation (INLA), a highly accurate and efficient Bayesian computation technique (Rue et al. 2009). To identify regions of activation, we propose an excursions set method based on the joint posterior distribution of the latent fields, which eliminates the need for multiple comparisons correction. Finally, we address a gap in the existing literature by proposing a novel Bayesian approach for multi-subject analysis. The methods are validated and compared to the classical GLM through simulation studies and a motor task fMRI study from the Human Connectome Project. The proposed Bayesian approach results in smoother activation estimates, more accurate false positive control, and increased power to detect truly active regions.

Figures

Equations

$$\mathbf{y}=\sum _{k=0}^K{\mathbf{X}}_k{\mathbf{\beta }}_k+\sum _{j=1}^J{\mathbf{Z}}_j{\mathbf{b}}_j+\mathbf{\epsilon },\mathbf{\epsilon }\sim N\left(0,\mathbf{V}\right).$$

$${\theta }_1=\log (\xi ),{\theta }_k=\log \left(\frac{1+{\varphi }_{k-1}}{1-{\varphi }_{k-1}}\right),$$

$$\text{mbox}Cov(\mathbf{u},\mathbf{v})=\frac{{\sigma }^2}{2^{\nu -1}\Gamma (\nu )}(\kappa \Vert \mathbf{u}-\mathbf{v}\Vert {)}^{\nu }{K}_{\nu }(\kappa \Vert \mathbf{u}-\mathbf{v}\Vert ),$$

$${\left({\kappa }^2-\Delta \right)}^{\alpha /2}\left(\tau \beta (\mathbf{u})\right)=\mathcal{W}(\mathbf{u}),\mathbf{u}\in \text{mbox}{\mathrm{I}\mathrm{R}}^d$$

$$\beta (\mathbf{u})\approx \sum _{i=1}^n{\psi }_i(\mathbf{u})w_i.$$

$${\mathbf{Q}}_{\kappa ,\tau }={\tau }^2\left({\kappa }^4\mathbf{C}+2{\kappa }^2\mathbf{G}+\mathbf{G}{\mathbf{C}}^{-1}\mathbf{G}\right).$$

$$\mathbf{\beta }=\mathbf{\Psi }\mathbf{w},\mathbf{w}\mid \kappa ,\tau \sim N\left(0,{\mathbf{Q}}_{\kappa ,\tau }^{-1}\right),$$

$$\mathbf{y}\mid {\mathbf{\beta }}_k,{\mathbf{b}}_j,\mathbf{\theta }\sim N\left({\mathbf{\mu }}_y,\mathbf{V}\right),{\mathbf{\mu }}_y=\sum _{k=0}^K{\mathbf{X}}_k{\mathbf{\beta }}_k+\sum _{j=1}^J{\mathbf{Z}}_j{\mathbf{b}}_j,$$

$${\mathbf{\beta }}_k={\mathbf{\Psi }}_k{\mathbf{w}}_k,{\mathbf{w}}_k\mid \mathbf{\theta }\sim N\left(0,{\mathbf{Q}}_{{\kappa }_k,{\tau }_k}^{-1}\right),$$

$${\mathbf{b}}_j\sim N(0,\delta \mathbf{I}),\mathbf{\theta }\sim \pi (\mathbf{\theta }),$$

$$\tilde{\pi }(\mathbf{\theta }|\mathbf{y})\propto {\frac{\pi (\mathbf{f},\mathbf{\theta },\mathbf{y})}{{\pi }_G(\mathbf{f}|\mathbf{\theta },\mathbf{y})}|}_{\mathbf{f}={\mathbf{f}}^{\star }(\mathbf{\theta })},$$

$$\tilde{\pi }({\theta }_j|\mathbf{y})$$

$$\tilde{\pi }(f_i|\mathbf{y})$$

$$\tilde{\pi }(f_i|\mathbf{y})\approx \sum _{\ell }{\lambda }_{\ell }\tilde{\pi }(f_i|{\mathbf{\theta }}_{\ell },\mathbf{y}),$$

$$E_{\gamma ,\alpha }^{+}(f)=\underset{D}{\mathrm{argmax}}\{|D|:P(D\subseteq A_{\gamma }^{+}(f))\ge 1-\alpha \},$$

$$\pi (\mathbf{f}\mid \mathbf{y})=\int \pi (\mathbf{f}\mid \mathbf{y},\theta )\pi (\theta \mid \mathbf{y})d\theta .$$

$${\mathbf{y}}_m$$

$${\tilde{\mathbf{y}}}_m={\tilde{\mathbf{X}}}_m{\mathbf{\beta }}_m+{\mathbf{\epsilon }}_m,{\mathbf{\epsilon }}_m\sim N(0,\mathbf{V}),$$

$$\tilde{\mathbf{y}}=\tilde{\mathbf{X}}\mathbf{\beta }+\mathbf{\epsilon },\mathbf{\epsilon }\sim N(0,\mathbf{\Sigma }),$$

$$\mathbf{\beta }={\mathbf{X}}_G{\mathbf{\beta }}_G+{\mathbf{\epsilon }}_G,{\mathbf{\epsilon }}_G\sim N\left(0,{\sigma }_G^2\mathbf{I}\right)$$

$$\tilde{\mathbf{y}}=\tilde{\mathbf{X}}{\mathbf{X}}_G{\mathbf{\beta }}_G+\tilde{\mathbf{X}}{\mathbf{\epsilon }}_G+\mathbf{\epsilon }.$$

$${\hat{\mathbf{\beta }}}_G=\frac{1}{L}\sum _{\ell =1}^L{\hat{\mathbf{\beta }}}_G^{(\ell )}\text{mbox}and{F}_{\gamma }=\frac{1}{L}\sum _{\ell =1}^L{F}_{\gamma }^{(\ell )},$$

$${\mathbf{\mu }}_m={\mathbf{Q}}_m^{-1}{\mathbf{\Psi }}_m^{\prime }{\tilde{\mathbf{X}}}_m^{\prime }{\mathbf{V}}^{-1}{\tilde{\mathbf{y}}}_m,{\mathbf{Q}}_m={\mathbf{Q}}_{\theta }+{\mathbf{\Psi }}_m^{\prime }{\tilde{\mathbf{X}}}_m^{\prime }{\mathbf{V}}^{-1}{\tilde{\mathbf{X}}}_m{\mathbf{\Psi }}_m.$$

$${\mathbf{\beta }}_G\mid \mathbf{y},\mathbf{\theta }\sim N\left(\mathbf{A}\mathbf{\Psi }\mathbf{\mu },\mathbf{A}\mathbf{\Psi }{\mathbf{Q}}^{-1}{\mathbf{\Psi }}^{\prime }{\mathbf{A}}^{\prime }\right)$$

$$\pi ({\mathbf{\beta }}_G\mid \mathbf{y})=\int \pi ({\mathbf{\beta }}_G\mid \mathbf{y},\mathbf{\theta })\pi (\mathbf{\theta }\mid \mathbf{y})d\mathbf{\theta }.$$

$$\pi (\mathbf{y}\mid \mathbf{\theta })\propto \pi (\mathbf{\theta }{)}^{-M}\prod _{m=1}^M\pi (\mathbf{\theta }\mid {\mathbf{y}}_m),$$

$$\tilde{\pi }({\mathbf{\beta }}_G\mid \mathbf{y})=\int \pi ({\mathbf{\beta }}_G\mid \mathbf{y},\mathbf{\theta })\pi (\mathbf{\theta }{)}^{1-M}q(\mathbf{\theta }\mid \mathbf{y})d\mathbf{\theta }.$$

$$F_{\gamma }=\frac{1}{{\sum }_{\ell =1}^L{\lambda }_{\ell }}\sum _{\ell =1}^L{\lambda }_{\ell }{F}_{\gamma }^{(\ell )},$$

$$y_v(t)={\beta }_{0v}+a_{1v}{z}_1(t){\beta }_{1v}+a_{2v}{z}_2(t){\beta }_{2v}+{\epsilon }_v(t),$$

$${\mathbf{y}}_m=\sum _{k=0}^6{\mathbf{X}}_{mk}{\mathbf{\beta }}_{mk}+\sum _{j=1}^{14}{\mathbf{Z}}_{mj}{\mathbf{b}}_{mj}+{\mathbf{\epsilon }}_m,{\mathbf{\epsilon }}_m\sim N(0,{\sigma }_m^2\mathbf{I}),$$

Full text

Keywords: spatial statistics; smoothing; integrated nested Laplace approximation; stochastic partial differential equation; brain imaging

1   INTRODUCTION

Functional magnetic resonance imaging (fMRI) is a popular noninvasive neuroimaging technique commonly used to localize regions of the brain activated by a task or stimulus (Lindquist 2008, Poldrack et al. 2011). Functional MRI indirectly measures neuronal activity through the BOLD (blood oxygenation level dependent) response, an indicator of hemodynamic changes that occur following neuronal activation (Lindquist et al. 2009). Traditional or volumetric fMRI data consists of a time series of three-dimensional brain volumes, each composed of hundreds of thousands of equally sized volumetric elements (voxels). While neuronal activity is known to occur in gray matter, volumetric fMRI data consists largely of a number of other tissue classes, including white matter and cerebral spinal fluid.

Recently, an alternative representation of fMRI data in which the cortical gray matter is represented as a 2-dimensional manifold surface has experienced a rise in popularity (Fischl 2012, Glasser et al. 2013). The process of transforming volumetric to cortical surface fMRI (cs-fMRI) is illustrated in Figure LABEL:fig:vol_to_surf. First, a high-dimension structural image is used to identify the cortical gray matter ribbon (Dale et al. 1999). Second, a mesh is applied to the white matter surface, the internal boundary of the cortical gray matter, to form a 2-dimensional manifold within each hemisphere, which is then geometrically smoothed. Third, the surface is inflated to a sphere while minimizing distance distortions, in which format subjects are registered to a standard template space by aligning anatomical folding patterns (Fischl et al. 1999). Finally, the same volume-to-surface transformation is applied to each fMRI volume to obtain a cs-fMRI time series. The resulting data is a triangular mesh consisting of approximately 30,000 vertices in each hemisphere.

Cs-fMRI offers several advantages over volumetric fMRI, including better whole-brain visualization, dimension reduction, removal of extraneous tissue types, and improved alignment of cortical areas across subjects. However, perhaps the most important benefit of cs-fMRI is the greater neurobiological significance of distances in the cortical surface space. Nearby locations in cs-fMRI are close in terms of distance along the cortex, and therefore tend to exhibit similar patterns of neuronal activity, while in volumetric fMRI locations that are close in terms of Euclidean distance may be neurobiologically quite dissimilar, coming from different areas of the cortex or even different tissue classes as illustrated in Figure LABEL:fig:volumetric. Thus, cs-fMRI is more appropriate for use with methods that pool information across neighbors, such as spatial smoothing or Bayesian spatial models.

The traditional task analysis method for both volumetric and cs-fMRI data is the classical general linear model (GLM), also known as the “massive univariate” approach, in which a linear regression model relating the observed fMRI data to the expected BOLD response to each task is fit separately at each location (e.g. voxel or vertex) in the brain (Friston et al. 1994, Hagler et al. 2006). Prior to model fitting, the fMRI data is typically smoothed using a fixed-width Gaussian kernel in order to increase the signal-to-noise ratio (SNR) of the data. The coefficients of the model are estimates of the task-related activation at each location, the significance of which is tested using a $t$ or $F$ statistic. The corresponding $p$-values are then plotted at each location to form a statistical parametric map (SPM). To identify the areas of true activation, the null hypothesis of no activation at each location is tested by thresholding the SPM at significance level $\alpha$ (Worsley and Friston 1995), chosen to control the family-wise error rate (FWER) or false discovery rate (FDR) at some predetermined level (Genovese et al. 2002, Lindquist and Mejia 2015). To account for correlations between tests performed at neighboring locations, popular solutions include parametric methods such as random field theory (RFT) (Adler 1981) and nonparametric methods such as permutation tests (Nichols and Holmes 2002). To avoid identifying very small regions of spurious activation, cluster-based methods first threshold the SPM at a fixed level determined by the researcher (e.g. $p=0.01$ or $p=0.001$) then use permutation tests, random field theory, or ad-hoc methods to determine significant clusters (Poline and Mazoyer 1993).

While the effort to properly correct for multiple comparisons in the classical GLM is a noble and necessary one—one which is unfortunately still not universally practiced (Carp 2012)—most of the traditional correction methods have been shown to suffer from various pitfalls. In high-dimensional settings, methods that control the FWER or FDR have been shown to suffer from a lack of power to detect true effects (Ishwaran and Rao 2003, Marchini and Presanis 2004). In the context of fMRI, parametric methods such as RFT have been found to be inaccurate due to departures of the data from the parametric assumptions (Nichols and Hayasaka 2003, Wager et al. 2009, Eklund et al. 2012; 2016). Further, as voxels (or vertices) lack biological meaning as a unit of measure, controlling the voxel-wise (or vertex-wise) FWER is somewhat arbitrary and is sensitive to voxel size, which is gradually shrinking due to technological advances. While cluster-based methods avoid this limitation, they have been found to be sensitive to the choice of initial threshold (Woo et al. 2014, Eklund et al. 2016). Finally, since the corrected significance threshold becomes more conservative as the number of tests increases, inference is sensitive to the size of the search volume.

These issues are symptomatic of some of the fundamental limitations of the classical GLM, described previously by Friston et al. (2002) and Friston and Penny (2003), among others. First, while it is well-known that the activation amplitude of one voxel depends on its neighbors, the classical GLM does not account for such dependence, since the model is estimated separately at every location. Second, while spatial smoothing of the fMRI data prior to model fitting can increase SNR and help satisfy the assumptions of RFT, when applied to volumetric fMRI data it may also combine signal from different tissue types and across discontiguous regions of the cortex (see Figure LABEL:fig:volumetric), contaminating the signal of interest and leading to inaccurate identification of truly active regions. Smoothing of cs-fMRI data is less problematic but also tends to blur boundaries between active and non-active areas. Furthermore, smoothed data may still fail to exhibit the Gaussian spatial autocorrelations required for RFT (Wald and Polimeni 2016). Smoothing also increases dependence between tests, complicating the problem of correcting for multiple comparisons.

To remedy these problems, several Bayesian alternatives to the classical GLM have been proposed (Friston and Penny 2003, Zhang et al. 2015). In a Bayesian GLM, specific prior distributions are assumed for the latent task activation fields and other unknown parameters in the model, and they, together with the likelihood, form a Bayesian hierarchical model. A variety of Bayesian computation techniques can be used to fit the model and obtain posterior estimates. For each location, the posterior probability that the corresponding amplitude is greater than some biologically meaningful activation threshold (often a percentage of global mean signal) is calculated, and active locations are identified by thresholding the resulting posterior probability map (PPM) at a certain level (e.g., 0.95). The main challenges for the successful implementation of the Bayesian GLM are selecting an appropriate prior on latent task activation fields, performing the Bayesian computation efficiently, using the joint posterior distribution of each latent field to identify regions of activation, and performing multi-subject analysis. We now address each of these in turn.

To account for the spatial dependence in activation levels, the model coefficients are often assumed to follow a spatial process prior, which for computational purposes should have a sparse inverse covariance structure. Several such priors have been advocated for volumetric fMRI, including a first-order Gaussian Markov random field (GMRF) prior (Gössl et al. 2001, Quirós et al. 2010), Laplacian GMRF prior (Penny et al. 2005), sparse spatial basis function prior (Flandin and Penny 2007), diffusion-based spatial prior (Harrison et al. 2008), Ising MRF prior (Smith and Fahrmeir 2007), and nonstationary GMRF prior (Yue and Speckman 2010, Yue et al. 2010; 2012). Many priors that have been applied to volumetric fMRI are designed for data on a regular lattice and may not be applicable to cs-fMRI data, which is in the form of a triangular mesh. In this paper, we employ a novel class of GMRF priors introduced by Lindgren et al. (2011). The priors are obtained by solving certain stochastic partial differential equations (SPDEs) using a finite element method. Compared to the previously proposed spatial priors, the SPDE priors possess several advantages for analyzing fMRI data. First, unlike regular GMRF priors, the SPDE priors are explicit mappings of Matérn Gaussian fields, which have been extensively used in statistical modeling of spatial data (e.g., Guttorp and Gneiting 2006). Therefore, they combine the computational advantages of a GMRF with the flexibility of the Matérn covariance structure in representing the dependence between locations. The Matérn parametrization also provides an intuitive interpretation of the spatial properties of the random fields. Second, the SPDE priors are constructed on a flexible triangular mesh, the structure of cs-fMRI data, and can apply appropriate smoothing along boundaries. Finally, SPDE priors are consistent under re-triangulations of the surface, which is not true in general of GMRF models that are not defined as discretisations of continuous models (Lindgren et al. 2011).

The second challenge comes from the feasibility of Bayesian computation. Functional MRI data are often massive, consisting of 100,000-200,000 voxels, or 60,000 vertices in cs-fMRI, observed at hundreds of time points for each subject. In our motor task fMRI study, we have $284$ time points, so that each subject’s data consists of over $17$ million observations. Therefore, standard Markov chain Monte Carlo (MCMC) methods are typically too time-consuming to be practically useful for data of this size (e.g., Woolrich et al. 2004a). Much attention has instead been directed towards variational Bayesian (VB) techniques (Woolrich et al. 2004b, Penny et al. 2005), which achieve computational advances by assuming independence between parameters, locations, or both. However, it has been well-established that VB techniques tend to severely underestimate posterior variance in latent Gaussian models (Wang and Titterington 2005, Bishop 2006, Rue et al. 2009, Sidén et al. 2017). We instead employ a more recently developed Bayesian inference tool based on integrated nested Laplace approximation (INLA) (Rue et al. 2009). The INLA method can directly compute accurate approximations to the posterior distributions and is able to handle large data sets by taking advantage of the sparsity of GMRFs. It is much faster than MCMC (Rue et al. 2016) and can be easily implemented using the R-INLA package (http://www.r-inla.org).

The third challenge lies in activation identification. The usual PPM approach introduced by Friston and Penny (2003) has several limitations. First, as the marginal rather than the joint posterior distribution is used to determine active locations, spatial dependencies are not fully reflected in the posterior activation probabilities. Second, using the marginal posterior distribution introduces a multiple comparisons problem, as each location is considered separately. Finally, while each of the locations in an identified active region has marginal activation probability greater than $1-\alpha$, where $\alpha$ is a predetermined significance level, the probability that the entire region is active (i.e., all locations in the region are active) may be less than $1-\alpha$. Using the joint posterior distribution of each latent activation field would eliminate these issues, but doing so is typically considered computationally infeasible. Here, we propose a computationally efficient joint PPM approach based on the excursions set method introduced by Bolin and Lindgren (2015) to identify the active region, defined as largest region such that with probability $1-\alpha$ all locations in the region are active.

Finally, while group-level inference is often a primary goal in task fMRI studies, existing spatial Bayesian models for volumetric fMRI data are typically designed for single-subject analysis due to computational limitations (Sidén et al. 2017). To address this gap in the literature, we propose a novel approache for combining data from multiple subjects to estimate the group-level posterior distribution of each latent activation field. Group-level regions of activation can then also be obtained using the proposed joint PPM approach.

The remainder of this paper is organized as follows. The Bayesian GLM method is introduced in Section 2, where the SPDE priors, INLA algorithm and joint PPM approach are presented. The method is then extended to multi-subject analysis in Section 3. We assess the accuracy of the proposed Bayesian methods and compare their performance to the classical GLM in a simulation study described in Section 4, followed by an application to a motor task fMRI study in Section 5. We conclude with a discussion in Section 6.

2   SINGLE SUBJECT BAYESIAN GLM

Let $T$ be the number of time points in the fMRI timeseries and let $N$ be the number of surface vertices in each hemisphere of the brain. For a subject and hemisphere, we have the following model:

[TABLE]

Here $\bm{y}$ is an $TN\times 1$ vector containing the fMRI time series of all vertices, and the $\bm{X}_{k}$ and $\bm{Z}_{j}$ are $TN\times N$ design matrices for the activation amplitudes $\bm{\beta}_{k}$ (including baseline $\bm{\beta}_{0}$) and nuisance signals $\bm{b}_{j}$, respectively. The matrix $\bm{V}=\bm{I}_{N}\otimes\bm{\Sigma}(\xi,\bm{\phi})$, where $\bm{\Sigma}(\xi,\bm{\phi})$ is a $T\times T$ covariance matrix for an AR($p$) process with marginal precision $\xi$ and partial autocorrelation functions $\bm{\phi}=(\phi_{1},,\dots,\phi_{p})^{\prime}$ assumed for each time series, and $\otimes$ denotes the Kronecker product. For fully Bayesian inference, prior distributions are assumed on the unknown parameters in model (1). For nuisance parameters in $\bm{b}_{j}$, we take independent and diffuse Gaussian priors, that is $\bm{b}_{j}\sim N({\bm{0}},\delta\bm{I})$ where $\delta$ is a fixed large number. For $\xi$ and $\bm{\phi}$, we first reparameterize them as

[TABLE]

for $k=2,\ldots,{p+1}$. We then let $\theta_{1}$ follow a log gamma prior and $(\theta_{2},\ldots,\theta_{p+1})$ a multivariate normal prior. The priors for $\bm{\beta}_{k}$ are described below.

2.1   SPDE Spatial Priors

To account for spatial homogeneity, we need to take a spatial prior on each $\bm{\beta}_{k}$ for $k=0,\ldots,K$. A good candidate is the class of Matérn Gaussian fields that has been extensively used in spatial statistics due to its flexible covariance function between locations. We say $\beta(\bm{u})$ is a Matérn Gaussian process if the covariance between $\bm{u}$ and $\bm{v}$ ($\bm{u},\bm{v}\in\mbox{\rm I\kern-1.99997ptR}^{d}$) is given by

[TABLE]

where $K_{\nu}(\cdot)$ is the modified Bessel function of the second kind with order $\nu>0$, $\Gamma(\cdot)$ is the gamma function, $\kappa>0$ is the spatial scale, and $\sigma^{2}>0$ is the variance. However, a Matérn spatial process is not computationally feasible for large data sets because its covariance matrix is completely dense and therefore difficult to invert. Lindgren et al. (2011) addressed this issue by deriving an explicit GMRF representation for Matérn Gaussian fields through solving the following stochastic partial differential equation (SPDE)

[TABLE]

where $\Delta=\sum_{i=1}^{d}{\partial^{2}}/{\partial u_{i}^{2}}$ is the Laplacian operator, $\alpha$ is the parameter that affects the smoothness, and $\tau$ relates to the variance of $\beta$. On the right hand side of the equation, $\mathcal{W}(\bm{u})$ is the Gaussian white noise process. The stationary solution $\beta$ to this SPDE is a Matérn Gaussian field, and the link to the smoothness $\nu$ and variance $\sigma^{2}$ is $\nu=\alpha-d/2$ and $\sigma^{2}=\Gamma(\nu)\left[\Gamma(\alpha)(4\pi)^{d/2}\kappa^{2\nu}\tau^{2}\right]^{-1}$. Spectral theory shows that an integer $\alpha$ must be chosen to obtain a Markov field. We thus let $\alpha=2$, resulting in $\nu=1$ for the smoothness of a two-dimensional field.

To obtain a Markov structure, we approximate $\beta(\bm{u})$ using the following basis expansion

[TABLE]

Here $\psi_{i}$ is the piecewise linear function defined on a triangular mesh, taking a value of 1 at the $i^{\text{th}}$ vertex and 0 at all other vertices; $w_{i}$s are the random weights that need to be estimated, and $n$ is the number of vertices in the mesh. A typical mesh is chosen to maximize the minimum interior triangle angle to ensure smooth transitions between small and large triangles. The vertices are often chosen to be the data locations, and extra vertices are added heuristically to minimize the total number of triangles needed to fulfill the size and shape constraints of the function domain. The interpretation of (3) is that the weights $\bm{w}=(w_{1},w_{2},\ldots,w_{n})^{\prime}$ determine the value of the field at each vertex, and the values in the interior of the triangles are determined by linear interpolation.

The joint distribution of $\bm{w}$ is chosen so that the distribution of the functions $\beta(\bm{u})$ approximates the distribution of solutions to the SPDE (2). The result is that $\bm{w}$ is Gaussian with zero means and a sparse precision matrix given by

[TABLE]

where $\bm{G}$ is a sparse symmetric $n\times n$ matrix with non-zero entries in cells corresponding to neighboring locations, and $\bm{C}$ is a diagonal matrix (Bolin and Lindgren 2013). Consider $N$ data locations $\bm{u}_{i}~{}(i=1,\ldots,N)$ and let vector $\bm{\beta}$ contain a realization of the random field at those locations. Then, based on (3) the SPDE prior on $\bm{\beta}$ is given by

[TABLE]

where $\bm{\Psi}$ is the $N\times n$ sparse matrix of the basis functions. Note that $\bm{\Psi}$ is the identity matrix if the data locations are the vertices in the mesh. Regarding the hyperparameters $\kappa$ and $\tau$, we take independent log-normal priors with zero mean and unit variance, that is $\log(\kappa)\sim N(0,1)$ and $\log(\tau)\sim N(0,1)$.

As shown in Lindgren et al. (2011), this is the best piecewise linear approximation to the continuous solution to the SPDE (2) given a triangular mesh. Since it is a GMRF representation of Matérn Gaussian fields, the SPDE prior allows us to capture both the spatial correlation and spatial smoothness that exist in a spatial process. The SPDE prior is particularly well-suited for cs-fMRI data, since the data are structured as a triangular mesh, in which connected vertices represent nearest neighbors along the cortical surface that would be expected to display similar pattens of neuronal activation. In SPDE priors, dissimilarity between locations is related to the distance along the cortical surface, rather than Euclidean distance as in regular Matérn models. The Laplacian prior used in Penny et al. (2005) is a special case of the SPDE prior on a regular lattice.

2.2   Approximate Inference by INLA

Based on model (1) and the priors specified in the previous sections, we may construct a Bayesian hierarchical model for our fMRI analysis as follows

[TABLE]

where $\bm{\theta}=(\xi,\phi_{1},\ldots,\phi_{p},\kappa_{1},\ldots,\kappa_{k},\tau_{1},\ldots,\tau_{k})$ contains all the hyperparameters, and $\pi(\bm{\theta})$ denotes the joint density of their hyperpriors.

It is possible to derive the full conditional distribution of each unknown parameter, and then use MCMC-based algorithms to obtain quite a few samples from their posterior distributions and make Bayesian inferences using those samples. However, the MCMC may have mixing problems and be slow to converge, considering the size of our data and the complexity of our model. As an alternative to MCMC, Rue et al. (2009) introduced a novel Bayesian computation tool based on integrated nested Laplace approximations (INLA), which is implemented in the R-INLA package (Martins et al. 2013). The INLA method can handle general latent Gaussian hierarchical models, including the model proposed in this paper. It accurately approximates marginal posterior densities and computes all necessary estimates faster than MCMC techniques.

A typical latent Gaussian hierarchical model has a set of hyperparameters $\bm{\theta}$ with prior $\pi(\bm{\theta})$, a set of latent Gaussian variables $\bm{f}$ with prior $\pi(\bm{f}|\bm{\theta})$ and a set of response variables $\bm{y}$ with likelihood $\pi(\bm{y}|\bm{f},\bm{\theta})$. The joint posterior distribution is then given by $\pi(\bm{f},\bm{\theta}|\bm{y})\propto\pi(\bm{y}|\bm{f},\bm{\theta})\pi(\bm{f}|\bm{\theta})\pi(\bm{\theta})$. The INLA method first approximates the marginal posterior of $\bm{\theta}$ as follows:

[TABLE]

where $\pi_{G}(\bm{f}|\bm{\theta},\bm{y})$ is the Gaussian approximation to the full conditional of $\bm{f}$, and $\bm{f}^{\star}(\bm{\theta})$ is the mode of that distribution. Then the approximated marginals are given by

[TABLE]

where $\bm{\theta}_{-j}$ denotes a subvector of $\bm{\theta}$ without element $\theta_{j}$. The approximated marginal of $\theta_{j}$ can be obtained by summing out the remaining variables $\bm{\theta}_{-j}$ from $\tilde{\pi}(\bm{\theta}|\bm{y})$. If needed, the approximated marginal of $f_{i}$ is obtained by first approximating the full conditional of $f_{i}$ with another Laplace approximation. The parameters are then numerically integrated out from $\tilde{\pi}(f_{i}|\bm{\theta},\bm{y})$, which gives

[TABLE]

where $\lambda_{\ell}$ are proportional to $\tilde{\pi}(\bm{\theta}_{\ell}|\bm{y})$. The evaluation points $\bm{\theta}_{\ell}$ can be chosen in different ways, depending on the relative importance of computational efficiency and accuracy in a given setting (Martins et al. 2013).

For the proposed Bayesian GLM model, given $\bm{f}=(\bm{\beta}_{1}^{\prime},\ldots,\bm{\beta}_{K}^{\prime},\bm{b}_{1},\ldots,\bm{b}_{J}^{\prime})^{\prime}$, $\pi(\bm{y}|\bm{f},\bm{\theta})$ is the Gaussian likelihood function defined in (2.2); $\pi(\bm{f}|\bm{\theta})$ is the joint multivariate Gaussian distribution of the independent priors specified on $\bm{b}_{j}$ and $\bm{\beta}_{k}$; and $\pi(\bm{\theta})$ is the joint distribution of the hyperpriors. Since the likelihood is Gaussian for this model, the INLA method is greatly simplified and the only necessary approximation is the numerical integration over $\tilde{\pi}(\bm{\theta}|\bm{y})$.

2.3   Joint PPM for Activation Identification

After fitting model (2.2) with INLA, we may use the resulting estimates of activation amplitudes to identify activated brain regions. There exist a number of threshold adjustment techniques for doing this. Most of these techniques are based on first calculating the marginal probabilities $P(f(\bm{u})>\gamma)$, where $\gamma$ is an activation threshold, then defining the excedence region as $D=\{\bm{u}:P(f(\bm{u})>\gamma)>1-\alpha\}$, where $\alpha$ is some significance level. However, the value of $\alpha$ needs to be adjusted for multiple-comparison issue, which is typically done using Type I error control, false discovery rate thresholding, or posterior probability thresholding (Marchini and Presanis 2004). In this work, we instead employ the joint probabilities using the excursion set method introduced by Bolin and Lindgren (2015).

We define the positive excursion set $E_{{\gamma,\alpha}}^{{+}}$ as the largest set of vertices such that with at least probability $1-\alpha$ the level $\gamma$ is exceeded at all locations in that set, which we can write as

[TABLE]

where $A_{\gamma}^{+}(f)=\{\bm{u}\in\Omega;f(\bm{u})>\gamma\}$. Since the set is defined using the joint distribution of the random field, it should be calculated based on the posterior distribution

[TABLE]

Computing $E_{{\gamma,\alpha}}^{{+}}(f)$ based on this distribution is computationally demanding, but can be done efficiently using the INLA technique. To do so, note that $E_{{\gamma,\alpha}}^{{+}}(f)=A_{\alpha}^{+}(F_{\gamma}^{+})$ where $F_{\gamma}^{+}(\bm{u})=\sup\{1-\alpha,\bm{u}\in E_{{\gamma,\alpha}}^{{+}}(f)\}$ is the so-called excursion function. Using the INLA approximation of $\pi(\bm{f}|\bm{y})$, the excursion function is $F_{\gamma}^{+}(\bm{u})=\sum_{\ell=1}^{L}\lambda_{\ell}F_{\ell}(\bm{u})$ where $F_{\ell}(\bm{u})$ is the excursion function calculated for the conditional posterior $\pi(\bm{f}|\bm{y},{\bm{\mathrm{\theta}}}_{\ell})$ for a fixed parameter configuration ${\bm{\mathrm{\theta}}}_{\ell}$ with corresponding weights $\lambda_{\ell}$, as in (2.2).

For the Bayesian GLM model, $\pi(\bm{f}|\bm{y},{\bm{\mathrm{\theta}}}_{\ell})$ is Gaussian and the computation of $F_{\ell}(\bm{u})$, $\ell=1,\dots,L$ therefore only requires the ability to compute excursion probabilities of multivariate Gaussian distributions. This can be done efficiently using the sequential method described in Bolin and Lindgren (2015). We refer to Bolin and Lindgren (2015) for further details and note that the method is implemented in the R package excursions (Bolin and Lindgren 2016), which has an interface to R-INLA that greatly simplifies the usage of the method for the Bayesian GLM model.

3   Multi-subject Bayesian GLM

So far, we have been concerned with the analysis of a single subject’s data. In practice, researchers often want to estimate population-level effects or pool information across a group of subjects. Multi-subject fMRI data have a hierarchical nature, with lower-level observations (e.g., individual subjects) nested within higher levels (e.g., groups of subjects). It is therefore common to use a two-level model, where in the first level there are a relatively large number of autocorrelated observations (time points), while in the second level there are a relatively small number of independent, identically distributed observations (subjects). Hierarchical models have been proposed as a way to incorporate uncertainty from subject-level models in the group-level estimates of the parameters of interest (Beckmann et al. 2003, Woolrich et al. 2004c, Lindquist et al. 2012, Degras and Lindquist 2014). However, as these models are typically fit separately at each voxel, they fail to account for spatial dependence in the data. On the other hand, Bayesian approaches that account for spatial dependence are usually designed for single-subject analysis. A group-level Bayesian approach that accounts for spatial dependencies for the analysis of fMRI data is therefore needed (Sidén et al. 2017).

Theoretically, it is straightforward to combine data from all subjects to build a single model like (1) to estimate the population effects. Unfortunately, the corresponding design matrix usually has extremely large dimensions, making model fitting computationally prohibitive. Therefore, we first fit model (1) for each subject and then combine the results in a principled manner. Let us write the model for each subject $m=1,\ldots,M$ as

[TABLE]

where $\bm{y}_{m}$ contains all response observations for the $m^{\text{th}}$ subject, $\bm{\beta}_{m}=(\bm{\beta}^{\prime}_{m0},\bm{\beta}^{\prime}_{m1},\bm{\beta}^{\prime}_{m2},\ldots,\bm{\beta}^{\prime}_{mK})^{\prime}$ and $\bm{b}_{m}=(\bm{b}^{\prime}_{m1},\bm{b}^{\prime}_{m2},\ldots,\bm{b}^{\prime}_{mJ})^{\prime}$ contain regression coefficients, and $\bm{X}=[\bm{X}_{0},\bm{X}_{1},\bm{X}_{2}\ldots,\bm{X}_{J}]$ and $\bm{Z}_{m}=[\bm{Z}_{m1},\bm{Z}_{m2}\ldots,\bm{Z}_{mK}]$ are design matrices. Since only $\bm{\beta}_{m}$ is of interest, we may regress $\bm{Z}_{m}$ out of the response $\bm{y}_{m}$ and the task design matrix $\bm{X}_{m}$ to obtain the following model

[TABLE]

where $\tilde{\bm{y}}_{m}$ is the fitted residual in the regression of $\bm{y}_{m}$ against $\bm{Z}_{m}$, and $\tilde{\bm{X}}_{m}$ is the fitted residual in the regression of $\bm{X}_{m}$ against $\bm{Z}_{m}$.

3.1   Two-level Modeling Approach

A natural approach might be to extend the familiar two-level model to a spatial Bayesian setting. The first-level model is then the combination of the $M$ subject-level models in (8), that is

[TABLE]

where $\tilde{\bm{y}}=(\tilde{\bm{y}}_{1}^{\prime},\ldots,\tilde{\bm{y}}_{M}^{\prime})^{\prime}$, $\tilde{\bm{X}}=\mbox{diag}(\bm{X}_{1},\ldots,\bm{X}_{M})$, $\bm{\beta}=(\bm{\beta}_{1}^{\prime}\ldots,\bm{\beta}_{M}^{\prime})^{\prime}$, $\bm{\varepsilon}=(\bm{\varepsilon}_{1},\ldots,\bm{\varepsilon}_{M})^{\prime}$ and $\bm{\Sigma}=\bm{I}_{M}\otimes\bm{V}$. The second-level model is given by

[TABLE]

where $\bm{X}_{G}$ is the $NKM\times NK$ second-level design matrix, $\bm{\beta}_{G}$ is the $NK\times 1$ vector of group-level activations or contrasts, and $\sigma^{2}_{G}$ is the between-subject variance. In the case of a single group of subjects with no covariates, we have $\bm{X}_{G}={\bf 1}_{M}\otimes\bm{I}_{NK}$, but other model specifications are also possible. The two models can be collapsed as

[TABLE]

Unfortunately, it is computationally infeasible to fit this model using a spatial Bayesian framework. A naive solution is to first fit the $M$ subject-level models as described in Section 2 to obtain the posterior mean of $\bm{\beta}$ in (9), denoted $\hat{\bm{\beta}}$. Then, use $\hat{\bm{\beta}}$ as the response variable in model (10), fit this model using INLA with SPDE priors taken on the components of $\bm{\beta}_{G}$, and perform Bayesian inference on $\bm{\beta}_{G}$ as described in Section 2.

Using the posterior means, however, underestimates the uncertainty from estimating $\bm{\beta}$ in the subject-level models. A more accurate but computationally demanding solution is to generate posterior samples of $\bm{\beta}$ in model (9), based on the estimated posterior distribution of $\bm{\beta}_{m}$ in each subject-level model, $m=1,\dots,M$. Let $\bm{\beta}^{(\ell)}$ be the $\ell^{\text{th}}$ such sample for $\ell=1,\ldots,L$. For each sample, we use $\bm{\beta}^{(\ell)}$ as the response variable in model (10) and fit the model using INLA. Finally, we combine the estimates from those models by averaging over the $L$ models. For example, the point estimator for $\bm{\beta}_{G}$ and the excursion function for the posterior distribution of $\bm{\beta}_{G}$ are given by

[TABLE]

where $\hat{\bm{\beta}}_{G}^{(\ell)}$ and $F^{(\ell)}_{\gamma}$ are the posterior mean and the excursion function based on the $\ell^{\text{th}}$ sample.

While this sampling approach accounts for the uncertainty in the posterior distributions of the $\bm{\beta}_{m}$, it is computationally demanding due to the need to fit a separate INLA model to each of $L$ samples. More importantly, assuming SPDE priors on both the subject-level and population-level activation fields may lead to over-smoothing in the second-level model. Therefore, we propose the following joint modeling approach.

3.2   Joint Modeling Approach

For each subject $m$, $\bm{\beta}_{m}=\bm{\Psi}_{m}\bm{w}_{m}$, where $\bm{\Psi}_{m}=\mbox{diag}(\bm{\Psi}_{m0},\ldots,\bm{\Psi}_{mK})$ is block diagonal and $\bm{w}_{m}=(\bm{w}^{\prime}_{m1},\bm{w}^{\prime}_{m2},\ldots,\bm{w}^{\prime}_{mK})^{\prime}$. Taking an SPDE prior on each $\bm{\beta}_{mk}$ implies that $\bm{w}_{m}\stackrel{{\scriptstyle i.i.d.}}{{\sim}}N\left({\bm{0}},\bm{Q}_{\theta}^{-1}\right)$ with $\bm{Q}_{\theta}=\mbox{diag}(\bm{Q}_{\kappa_{0},\tau_{0}},\ldots,\bm{Q}_{\kappa_{K},\tau_{K}})$ and $\bm{Q}_{\kappa_{k},\tau_{k}}$ denoting the precision matrix of the SPDE prior taken on $\bm{\beta}_{mk}$. Then, the full conditional distribution of $\bm{w}_{m}$ is Gaussian with mean vector $\bm{\mu}_{m}$ and precision matrix $\bm{Q}_{m}$, given by

[TABLE]

Since the subject-level models are independent of each other, we let $\bm{w}=(\bm{w}_{1}^{\prime},\ldots,\bm{w}_{M}^{\prime})^{\prime}$ and derive the full conditional distribution of $\bm{w}$ to be Gaussian with mean $\bm{\mu}=\mbox{diag}(\bm{\mu}_{1},\ldots,\bm{\mu}_{M})$ and precision matrix $\bm{Q}=\mbox{diag}(\bm{Q}_{1},\ldots,\bm{Q}_{M})$. Letting $\bm{\beta}=(\bm{\beta}_{1}^{\prime},\ldots,\bm{\beta}_{M}^{\prime})^{\prime}$, we have $\bm{\beta}=\bm{\Psi}\bm{w}$, where $\bm{\Psi}=\mbox{diag}(\bm{\Psi}_{1},\ldots,\bm{\Psi}_{M})$. We then define the group-level $\bm{\beta}_{G}$ to be a linear combination of the subject-level $\bm{\beta}$, i.e., $\bm{\beta}_{G}=\bm{A}\bm{\beta}$, where $\bm{A}$ is a constant matrix determined by the nature of $\bm{\beta}_{G}$. For example, if we are interested in whether or not the $M$ subjects activate on average, we may define $\bm{\beta}_{G}=\sum_{m=1}^{M}\bm{\beta}_{m}/M$ and the corresponding $\bm{A}=(1/M,\ldots,1/M)\otimes\bm{I}_{N}$, an $NK\times NKM$ matrix. Since $\bm{\beta}_{G}=\bm{A}\bm{\Psi}\bm{w}$ and the full conditional of $\bm{w}$ is Gaussian, it is easy to see that the full conditional of $\bm{\beta}_{G}$ is also Gaussian

[TABLE]

where $\bm{\theta}$ denotes the set including all the hyperparameters in $\bm{\mu}$ and $\bm{Q}$. To obtain the posterior density of $\bm{\beta}_{G}$, we need to integrate out $\bm{\theta}$ in (11), i.e.,

[TABLE]

Unfortunately, it is hard (if not impossible) to explicitly solve the integral above. We therefore propose to numerically evaluate it using importance sampling. The marginal posterior density is $\pi(\bm{\theta}\mid\bm{y})\propto\pi(\bm{\theta})\pi(\bm{y}\mid\bm{\theta})$, and since the subjects are assumed to be independent given the parameters, it is easy to show that

[TABLE]

where $\pi(\bm{\theta})$ is the prior density of $\bm{\theta}$ as specified in Section 2, and $\pi(\bm{\theta}\mid\bm{y}_{m})$ is the posterior density of $\bm{\theta}$ from the $m^{\text{th}}$ subject. The proof for (12) is given in Appendix A. For computational ease, we approximate $\pi(\bm{\theta}\mid\bm{y}_{m})$ using the Gaussian approximation with mean $\tilde{\bm{\mu}}_{m}$ and precision $\tilde{\bm{Q}}_{m}$ given by INLA method when fitting the subject-level models. Then, the product term in (12) is proportional to another Gaussian density, denoted by $q(\bm{\theta}|\bm{y})$, with precision matrix $\sum_{m}\tilde{\bm{Q}}_{m}$ and mean $(\sum_{m}\tilde{\bm{Q}}_{m})^{-1}\sum_{m}\tilde{\bm{Q}}_{m}\tilde{\bm{\mu}}_{m}$. As a result, the posterior distribution of $\bm{\beta}_{G}$ can be approximated by

[TABLE]

Let $\bm{\theta}^{(\ell)}$ be the $\ell^{\text{th}}$ sample for $\ell=1,\ldots,L$ simulated from $q(\bm{\theta}|\bm{y})$. Then, $\tilde{\pi}(\bm{\beta}_{G}\mid\bm{y})$ can be approximated numerically as a mixture of Gaussians with weights $\lambda_{\ell}=\pi(\bm{\theta}^{(\ell)})^{1-M}$, and the posterior quantities of $\bm{\beta}_{G}$ can be easily computed. It’s noteworthy that the highly sparse $\bm{Q}$ in (11) makes the computation efficient in spite of its large dimension. Using this approximation, we can also compute the joint PPMs as described in Section 2.3. Specifically, the excursion function $F_{\gamma}$ for the joint posterior $\tilde{\pi}(\bm{\beta}_{G}\mid\bm{y})$ can be written as

[TABLE]

where $F_{\gamma}^{(\ell)}$ is the excursion function based on the Gaussian distribution $\pi(\bm{\beta}_{G}\mid\bm{y},\bm{\theta}_{\ell})$, which can be computed as described in Section 2.3, and $\lambda_{\ell}$ are the same weights as above. Finally, the regions of activation are identified by thresholding $F_{\gamma}$ at a certain level, e.g., 95%.

Since the proposed joint modeling approach is based on the full joint distribution of the population-level activation fields, it is more appropriate than the naive, two-level modeling approach. It also makes few additional assumptions beyond the Gaussian approximation of the posterior densities $\pi(\bm{\theta}\mid\bm{y}_{m})$, which is likely to introduce minimal loss of accuracy in our model. It is therefore expected to result in the most accurate estimates and regions of activation, and fitting the model is quite computationally feasible due to the sparsity of the full joint precision matrix. We compare the joint and two-level modeling approaches using experimental fMRI data in Section 5.

4   Simulation Study

We constructed a $46\times 55$ phantom image using a binary gray matter mask provided by the SPM 8 software package (Welcome Trust, UCL). The selected slice contains $1256$ voxels within the brain. The field-of-view of the image was assumed to be $192$ mm. A dynamic image time series of length $T=200$ was simulated as follows. Two task activation profiles $x_{1}$ and $x_{2}$, depicted in the first column of Figure LABEL:fig:act, were created by convolving a canonical hemodynamic response function (Friston et al. 1998) with two task-specific binary stimulus functions.

Activation amplitude maps $\bm{\beta}_{1}$ and $\bm{\beta}_{2}$ were formed by placing Dirac functions at three separate locations on the phantom image and smoothing each with Gaussian kernels with full width at half maximums (FWHMs) of $10$, $15$ and $20$ mm, respectively. Smoothing was only performed here to generate spatially smooth activation fields; no smoothing was performed on the actual simulated data for the Bayesian GLM. The signal at each voxel $v$ was generated by weighting $(x_{1},x_{2})$ by $(a_{1v},a_{2v})=(1,0)$ for voxels within the first region, $(1,0.5)$ within the second region, and $(0,0.5)$ within the third region. The resulting activation amplitudes and activated regions, consisting of locations with activation amplitude greater than zero, are shown in Figure LABEL:fig:act. The baseline field was based on the gray matter prior map, and AR(1) errors with unit variance and autocorrelation coefficient of 0.3 were independently generated at each voxel. In summary, the data at locations $v=1,\ldots,1256$ and time points $t=1,\ldots,200$ was simulated as

[TABLE]

where $\beta_{0v}$ is the baseline signal, $\beta_{1v}$ and $\beta_{2v}$ are the activation amplitudes, and $\varepsilon_{v}(t)$ is random autoregressive error. Example simulated images at two different time points are shown in Figure S1 in Appendix B.

We fit model (13) using the proposed Bayesian GLM method. Figure LABEL:fig:mesh displays the triangular mesh, where we see regular triangulation inside the brain and two boundary layers of larger triangles surrounding the brain and within the ventricles. Regions of activation were identified using the joint and marginal PPM approaches with excursion level $\gamma=0$ and significance levels $\alpha=0.05$ and $0.01$. For comparison purposes, we also applied the classical GLM method to the simulated data. We first spatially smoothed the data using a Gaussian kernel with $6$ mm FWHM, then fit a linear regression model with AR(1) errors to each voxel separately. The AR coefficient was estimated at each voxel by solving the Yule-Walker equations of the residuals of a linear model with uncorrelated errors, then averaging across all voxels. Regions of activation were identified by performing a $t$-test at every voxel, accounting for multiple comparisons with FDR correction or FWER correction. To control the FDR at $q=0.05$ and $q=0.01$, we used the Benjamini-Yekutieli procedure for dependent observations (Benjamini and Yekutieli 2001); to control the FWER at $\alpha=0.05$ and $\alpha=0.01$ we performed a permutation test by randomly reordering the prewhitened time series of all voxels $1000$ times, re-fitting the model after each reordering, to estimate the null distribution of the maximum $t$-statistic.

Model fitting was done by implementing the R-INLA package on a server with Intel Xeon 32-core CPU at 2.70GHz and 256GB of memory. Computation time to fit the model in (13) was less than 20 minutes and would be reduced to 5 minutes if independent Gaussian errors were assumed instead of AR(1) errors. Therefore, if the temporal pattern is not of interest, it may be removed from the data in the pre-processing stage to improve computational efficiency. This approach is adopted for our experimental data analysis, described in Section 5. Using the excursions package, the computation time to estimate the excursion function based on the joint PPM was approximately 30 seconds.

Estimates of the two amplitude fields based on the Bayesian and classical GLMs are shown in the top row of Figure LABEL:fig:sim:estimates. Both sets of estimates capture the general spatial patterns of the true fields, but the classical estimates are more noisy and also markedly underestimate the activation magnitudes. The difference in performance of the two methods is more clearly illustrated in Figure LABEL:fig:sim_roc, which displays receiver operating characteristic (ROC) curves for the Bayesian and classical GLMs, based on varying the threshold of the excursion function and corrected $p$-value maps, respectively. These curves illustrate that the Bayesian approach achieves both high sensitivity and specificity, with over $0.998$ area under the curve (AUC) for both activations. In contrast, the classical GLM is not able to simultaneously achieve high sensitivity and high specificity, with AUC for each activation of $0.913$ and $0.874$ using FDR correction and $0.942$ and $0.931$ using FWER correction, respectively. Furthermore, the difference in performance is greatest when specificity is close to $1$. As high specificity is often a principal goal in fMRI task activation studies, this suggests that in practice, the proposed Bayesian approach will tend to achieve higher sensitivity than the classical GLM when the false positive rate, and hence specificity, is controlled at a fixed level.

The bottom row of Figure LABEL:fig:sim_roc compares the excursion function values obtained from the joint posterior probabilities and marginal posterior probabilities. The points show that the excursion function based on the joint posterior probabilities is a monotonic transformation of the excursion function based on the marginal posterior probabilities, leading to identical ROC curves for both approaches. This is due to the choice of parametric family implemented in the excursions R package, which has the monotonic property (Bolin and Lindgren 2015). However, the distribution of values obtained from the two approaches is very different. The excursion function values based on the joint PPM (histogram shown in blue) are clustered around 0 and 1, while the values based on the marginal PPM (histogram shown in red) are continuously distributed between 0 and 1. Therefore, besides providing more accurate control of false positives, the joint PPM is also less sensitive to the choice of significance threshold.

Figure LABEL:fig:sim99 displays estimated regions of activation obtained by thresholding the joint PPMs, marginal PPMs, and corrected $p$-values using a significance level of $0.01$ (see Figure S2 in Appendix B for regions of activation based on significance level $0.05$). The false positive rate (FPR) and false negative rate (FNR) are reported for each method and activation, relative to the true active regions shown in Figure LABEL:fig:act. The Bayesian joint PPM method significantly outperforms the classical GLM using either FDR or FWER correction, simultaneously achieving lower FPR and FNR for both activations. This illustrates that the classical GLM approach not only results in a higher rate of false positives, but also suffers from reduced power. Finally, as expected, the marginal PPM approach yields more false positives than the joint PPM approach because it fails to appropriately correct for multiple comparisons.

5   Experimental Data Results

5.1   Data Collection and Processing

We applied the proposed methods to fMRI data from $20$ randomly sampled subjects from the Human Connectome Project (HCP), a database of demographic, behavioral and neuroimaging data collected on over one thousand healthy adult subjects (http://humanconnectome.org). We employed the fully preprocessed cortical surface fMRI data available in the HCP 500-subject data release. A customized Siemens Skyra 3T scanner using a multi-band/multi-slice pulse sequence with an acceleration factor of eight was used to acquire fMRI scans with high temporal and spatial resolution (Moeller et al. 2010, Feinberg et al. 2010, Setsompop et al. 2012, Xu et al. 2012, Uğurbil et al. 2013). High-resolution structural $T_{1}$-weighted images were acquired for anatomical localization and transformation to standard space.

The data first underwent the HCP fMRISurface minimal preprocessing pipeline, which includes removal of spatial distortions, motion realignment, distortion correction, alignment to the structural image, bias field correction, intensity normalization and removal of extra-cerebral voxels (Glasser et al. 2013). These steps are standard practice in fMRI data analysis and are necessary to align the data into a common space and remove major sources of noise. The volumetric fMRI data was then mapped from the volumetric representation to the cortical surface manifold, and subjects registered to the standard template space, as illustrated in Figure LABEL:fig:vol_to_surf. To regularize the mapping process, a small degree of spatial smoothing was applied after transformation using a geodesic Gaussian surface smoothing algorithm with 2mm FWHM (Glasser et al. 2013). While the spherical representation is used for inter-subject registration and fitting the spatial Bayesian model, lesser degrees of inflation, which maintain basic brain shape and structure, are used here for display purposes.

During each subject’s $3.5$-minute fMRI run, a series of $284$ functional volumes were acquired, one every $0.72$ seconds. Subjects performed a series of five motor tasks during the acquisition, each preceded by a three-second visual cue (HCP 2016). The following five tasks were performed twice over the run for twelve seconds each time: move tongue, tap left fingers, tap right fingers, wiggle left toes, and wiggle right toes. For each type of task and the visual cue, a binary stimulus function was defined taking the value $1$ during the periods during which the subject was asked to perform the task or view the cue and [math] otherwise. These six binary stimulus functions were each convolved with a canonical hemodynamic response function (Friston et al. 1998), resulting in a set of stimulus response functions for each subject (see Figure S3 in Appendix C). Three 15-second periods during which the subject was asked to fixate on a crosshair displayed on the screen served as a baseline.

To account for noise due to subject motion, six rigid body realignment parameters estimated in the motion realignment phase of preprocessing and their first-order temporal derivatives were included in the model as nuisance covariates. Linear and quadratic time terms were also included to account for scanner drift. To reduce the computational cost of fitting the Bayesian model, the fMRI time courses were first pre-whitened by assuming an AR($p$) process on the residuals from a classical GLM with uncorrelated errors. A model order of $p=6$ was chosen based on inspection of the partial autocorrelation function. This relatively high model order is not surprising due to the fast temporal resolution of the data. The AR coefficients were estimated by solving the Yule-Walker equations and were allowed to vary spatially, as previous studies have shown that the degree of temporal autocorrelation is not constant across the brain (Worsley et al. 2002, Eklund et al. 2012). To improve estimation efficiency, the AR coefficient estimates at each location were averaged across all subjects in the analysis.

5.2   Model Estimation and Results

The left and right hemisphere each contain 32,492 surface vertices, including approximately 2,700 missing values in the medial wall. To reduce the computational burden while maintaining high spatial resolution, we used the Connectome Workbench to resample each hemisphere to approximately 6,000 vertices (Marcus et al. 2011). A spherical mesh was created based on the radial coordinates of the locations excluding missing values in the medial wall, resulting in a mesh with approximately 5,300 vertices. Figure LABEL:fig:app:mesh displays the lateral (external) and medial (internal) surfaces of the spherical mesh for each hemisphere. After model fitting, coefficient estimates and posterior probabilities for the 32,492 vertices for each hemisphere were obtained by projecting back to the original resolution using the Connectome Workbench tools.

For the $m^{\text{th}}$ ($m=1,\ldots,20$) subject, we considered the following model for each hemisphere:

[TABLE]

where $\bm{y}_{m}$ is a vector containing the $NT$ response values, where $T=284$ is the number of time points and $N\approx 5300$ is the number of locations in the mesh; the $\bm{X}_{mk}$ are the stimulus response functions corresponding to baseline, the visual cue, and the five motor tasks; and the $\bm{Z}_{mj}$ contain the nuisance covariates for motion and drift. Note that here uncorrelated errors are assumed because temporal autocorrelation was removed prior to model fitting through prewhitening. The SPDE spatial priors and diffuse Gaussian priors described in Section 2 were independently assumed on $\bm{\beta}_{mk}$ and $\bm{b}_{mj}$, respectively; other priors were specified as described in Section 2. We used the INLA method to fit model (14) and obtain posterior estimates for each subject. Using the joint and two-level models described in Section 3, we then combined those subject-level estimates to derive group-level posterior estimates. Finally, we used the excursion set method based on the joint PPM described in Section 2.3 to identify group-level regions of activation in response to each task. Two values of the activation threshold $\gamma$ were considered: [math] and 1% of the baseline signal, similar to Quirós et al. (2010).

For comparison purposes, we also applied the classical GLM method to the data. For each subject, we first spatially smoothed the data using the Connectome Workbench geodesic Gaussian surface smoothing algorithm with 6mm FWHM (Glasser et al. 2013). We then prewhitened the data as described above. We fit a classical GLM to each subject’s data including the 14 nuisance terms as in the Bayesian GLM. For each task, we then used the subject-level estimates in a group-level model with weights inversely proportional to the squared standard error of the subject-level coefficients. Finally, we identified group-level regions of activation by performing a $t$-test at every vertex, accounting for multiple comparisons with FDR correction or FWER correction similarly as in the simulation described in Section 4. To control the FWER at the group level, we performed a permutation test by generating $100$ resamplings of each subject’s prewhitened data, refitting the subject-level models for each resampling, then re-estimating the group-level model to obtain the null distribution of the maximum $t$ statistic in the image.

Figure LABEL:fig:app:beta_subj displays the Bayesian estimates of activation amplitude (posterior means) for the visual cue and tongue tasks, along with the classical GLM estimates, for one randomly selected subject. The remaining tasks are shown in Appendix C. Both sets of estimates are noisy, due to the high levels of noise in fMRI data and relatively short duration of the tasks performed for each individual subject ($30$ seconds for the visual cue; $24$ seconds for the tongue task); however, the Bayesian estimates are substantially smoother than the classical estimates, even though the data was smoothed prior to fitting the classical GLM but not the Bayesian GLM. Figure LABEL:fig:app:beta displays group-level estimates of activation amplitude based on the classical and Bayesian GLM approaches. For the Bayesian GLM, the results using the joint modeling approach and the two-level modeling approach with the sampling method described in Section 3.1 are both displayed. In both the subject-level and group-level estimates, the activation fields for each task are generally as expected: during the visual cue, the visual and orbitofrontal (decision making) regions are highly active, while the somatomotor (motor planning) region is activated to a lesser degree; during the tongue movement task, the area of the motor cortex associated with the tongue shows a strong degree of activation; and the remaining motor tasks activate the expected areas of the motor cortex (see Figures S4 and S5 in Appendix C). Compared with the classical GLM, the Bayesian GLM results in smoother estimates of activation, as it accounts for dependence in the degree of activation between neighboring locations. As hypothesized, the two-level modeling approach tends to result in somewhat oversmoothed estimates, compared with the joint modeling approach.

Figure LABEL:fig:app99 shows the regions of activation identified by thresholding the joint PPM excursion function for each task, along with those identified through the classical GLM with FDR and FWER correction, using significance level $0.01$ in each method. For the Bayesian GLM, regions of activation at activation threshold $\gamma=0$ are shown in purple, while those at $\gamma=0.027$, corresponding to $1\%$ of the global baseline signal, are shown in red. Several effects are clearly notable. First, the areas of activation based on controlling the FWER in the classical GLM are very conservative. This is a well-known issue with this approach that worsens as the number of locations increases, which is likely why FWER correction produces reasonable results in the simulation study but very conservative results on real fMRI data. FDR control results in significantly larger areas of activation, reflecting an increase in power while maintaining a known degree of false positive control. Second, the areas of activation based on $\gamma=0$ in the Bayesian GLM are similar to those based on FDR control in the classical GLM, but are smoother and somewhat larger. This likely reflects an increase in power resulting from incorporating spatial dependencies in the Bayesian GLM, as observed in our simulation study. Third, compared with the joint model, the areas of activation based on the Bayesian two-level model are smoother but similar in size. This illustrates that while the two-level approach tends to result in somewhat oversmoothed estimates, the sampling method described in Section 3.1 is able to appropriately account for uncertainty in the subject-level estimates when fitting the group-level model. Finally, the increased power in the joint PPM results in widespread areas of activation based on $\gamma=0$, while the more biologically meaningful threshold of $\gamma=0.027$ ($1\%$ of the baseline signal) results in more conservative areas of activation that correspond well to the highly activated areas shown in yellow in Figure LABEL:fig:app:beta.

These results illustrate the benefits of using a Bayesian framework to account for spatial dependencies in fMRI task activation studies. Compared with the traditional analysis techniques, the proposed Bayesian GLM approach results in smoother estimates of activation for individual subjects and groups of subjects, smoother active regions, and greater power to detect areas of activation, particularly more subtle activations.

6   Discussion

In this article, we have proposed a novel Bayesian GLM approach for analysis of cortical surface fMRI (cs-fMRI) data, which has recently experienced a rise in popularity compared to traditional, volumetric fMRI. While cs-fMRI offers several advantages over volumetric fMRI, most analyses of cs-fMRI still utilize the classical “massive univariate” GLM approach, which suffers from several pitfalls including dependence on naive smoothing methods and multiple comparisons correction techniques that may not provide accurate control over false positive rates. By contrast, the Bayesian GLM approach that we propose achieves appropriate smoothing and control of false positives in a model-based fashion and is, to our knowledge, the first spatial Bayesian method proposed for this type of data. The advantages of the proposed Bayesian approach in comparison with the classical GLM have been demonstrated through a simulation study and a motor task fMRI study from the Human Connectome Project.

In comparison with previously proposed Bayesian methods for volumetric fMRI data, the proposed Bayesian approach for cs-fMRI data offers several advantages. First, the assumptions of isotropy and stationarity are much more reasonable for cs-fMRI data, since geodesic distances along the cortex are more biologically meaningful than Euclidean distances within the volume. Second, a common computational strategy in Bayesian models for volumetric fMRI is to fit a separate model within each slice, an approach that has been shown to result in discontinuities in estimated amplitude maps and regions of activation (Sidén et al. 2017). Since cs-fMRI data is smaller than volumetric fMRI data and can be resampled in a principled way to further reduce dimensionality with minimal loss of information, it is feasible to fit a single model to all the locations within each hemisphere. Third, while most Bayesian methods for volumetric fMRI reduce computational burden by using VB, which is known to underestimate posterior variance, we employ INLA, a computationally efficient but highly accurate approximate Bayesian inference tool. Because INLA is less computationally demanding than MCMC, we are able to fit a complex model based on flexible SPDE spatial processes in order to more accurately capture the unique spatial dependency patterns of each latent field. Fourth, while previously proposed Bayesian methods use the marginal posterior distribution to identify regions of activation, we adopt a computationally efficient excursions set method that uses the joint posterior distribution to account for spatial dependencies in the posterior probabilities of activation. Finally, we introduce a novel approach for multi-subject analysis, which avoids the computational burden of fitting multiple subjects simultaneously, while providing proper group-level spatial Bayesian inference, a gap in the existing literature (Sidén et al. 2017). None of these techniques have, to the best of our knowledge, ever been proposed for analysis of cortical surface fMRI data.

Several limitations of the proposed methods should be noted. First, the SPDE priors we propose in this paper are suitable for modeling stationary and isotropic spatial processes. In volumetric fMRI data, it has been well-established that amplitude fields tend to present some non-stationary features such as varying degrees of activation (e.g., Harrison et al. 2008, Yue et al. 2010). While we largely avoid these issues through the use of cortical surface fMRI data, non-stationary features may exist due to the presence of distinct functional cortical areas and varying degress of smoothness within activated and non-activated areas. Lindgren et al. (2011) introduced non-stationarity to the SPDE models by allowing $\kappa$ and $\tau$ in (2) to depend on location, and Bolin and Lindgren (2011) and Fuglstad et al. (2014) generalized further using nested directional operators and non-isotropic Laplacians, respectively. Future work will investigate the feasibility and performance of non-stationary SPDE models in the Bayesian GLM framework.

Second, while the joint PPM approach described in Section 2.3 allows us to identify a subset of jointly activated locations, this subset may include isolated regions containing only a single voxel or vertex, which are unlikely to be truly activated in isolation. These could be excluded from the excursion set, thereby guaranteeing that the remaining regions have at least the target activation probability. However, this would result in overly conservative regions of activation. Future work should aim to develop methods to identify spatially connected regions of activation based on the joint PPM.

Finally, while the proposed methods assume that spatial dependence is inversely related to distance, a more complex network structure involving long-range connections may exist in the data. Future work will incorporate such long-range dependencies between different brain regions during task activity. We also plan to investigate the potential of an empirical Bayesian framework employing estimates of long-range and short-range dependence to reduce the computational burden compared with a fully Bayesian approach such as the one proposed in this paper. Additionally, the distances utilized in this paper are radial distances along the spherical surface. While the spherical representation of the cortical surface of each hemisphere provides a simple geometric setting while minimizing distance distortions, some level of distortion is inevitable. Future work will focus on utilizing geodesic distances along the true cortical surface of each subject to more accurately capture distance-related spatial dependencies.