InTechOpen uses cookies to offer you the best online experience. By continuing to use our site, you agree to our Privacy Policy.

Medicine » Diagnostics » "Assessment of Cellular and Organ Function and Dysfunction using Direct and Derived MRI Methodologies", book edited by Christakis Constantinides, ISBN 978-953-51-2723-9, Print ISBN 978-953-51-2722-2, Published: October 26, 2016 under CC BY 3.0 license. © The Author(s).

Chapter 4

Application of Spatial Bayesian Hierarchical Models to fMRI Data

By Kuo-Jung Lee
DOI: 10.5772/64823

Article top


The stimulus functions and the predicted BOLD signal changes used in the Benchmark example.
Figure 1. The stimulus functions and the predicted BOLD signal changes used in the Benchmark example.
Simulated and estimated binary and ρ images.
Figure 2. Simulated and estimated binary and ρ images.
Four different parcellation schemes. The areas in red are active areas corresponding to γv=1. The blue dashed line is used to partition the area into non-overlapping parts.
Figure 3. Four different parcellation schemes. The areas in red are active areas corresponding to γv=1. The blue dashed line is used to partition the area into non-overlapping parts.
Depiction of temporal occurrences of four different tasks.
Figure 4. Depiction of temporal occurrences of four different tasks.
Activated regions upon representation of different tasks.
Figure 5. Activated regions upon representation of different tasks.
Estimated probability of activation in response to any face task with activation areas in yellow.
Figure 6. Estimated probability of activation in response to any face task with activation areas in yellow.
The activation areas with a posterior probability greater than c, corresponding to FDA = 0.05, shown in red for different tasks.
Figure 7. The activation areas with a posterior probability greater than c, corresponding to FDA = 0.05, shown in red for different tasks.
Areas shown more activity in incongruent than congruent task without consideration of spatial dependence of temporal correlations in the model.
Figure 8. Areas shown more activity in incongruent than congruent task without consideration of spatial dependence of temporal correlations in the model.
Areas shown more activity in incongruent than congruent task with consideration of spatial dependence of temporal correlations in the model.
Figure 9. Areas shown more activity in incongruent than congruent task with consideration of spatial dependence of temporal correlations in the model.

Application of Spatial Bayesian Hierarchical Models to fMRI Data

Kuo-Jung Lee
Show details


Bayesian modelling has attracted great interest in cognitive science and offered a flexible and interpretable way to study cognitive processes using functional magnetic resonance imaging data. In this chapter, a spatial Bayesian hierarchical model is applied to an event-related fMRI study of cognitive control using the Simon test. We consider a sparse spatial generalized linear mixed-effects model to capture the spatial dependence among activated voxels and temporal parameters and to benefit computationally by reducing dimensionality. We demonstrate that the proposed model has the capability of identification of the brain areas related to cognitive tasks. Moreover, the reduction in the false positive rate is observed in the simulation study, and the relevant brain regions involved in processing cognitive control are clearly detected in a real-life fMRI example.

Keywords: Bayesian, functional magnetic resonance imaging, Markov chain Monte Carlo, spatial generalized linear mixed-effects model

1. Introduction

Functional magnetic resonance imaging (fMRI) has increasingly become an important and popular modality that allows researchers to investigate brain activity resulting from a particular stimulus [1]. In an fMRI experiment, a subject is asked to perform a task by responding to a series of stimuli that may involve a motor, sensory or cognitive task, then the MR machine records the changes in the blood oxygen level dependent (BOLD) of the brain across different time points, resulting in three-dimensional fMRI time-series images. Numerous statistical models have been proposed to allow researchers to detect localized regions activated during a task, to describe the networks required for a particular brain function or to assess physical characteristics elicited by cognitive processes in the brain. However, the growth in the complexity, the significant size and the hierarchical structure of fMRI data make it difficult to comprise a fully efficient computationally feasible statistical model to accurately explain the temporal and spatial characteristics of the data.

The standard approach used on fMRI data is known as statistical parametric mapping (SPM) [2], which applied either a voxel-wise t-test or F-test statistics. In order to obtain an activation map of the brain, the next step is to threshold the test statistics at a given overall error rate that leads to a major multiplicity problem. The most common way of solving this problem is to use Gaussian random field theory [3]. This technique is based on the assumption of a stationary Gaussian random field, which may not be satisfied in fMRI settings. Another limitation is that most current methods ignore at least one of the spatial or temporal relationships between observations. Ignoring either spatial or temporal correlations in the model leads to seriously biased conclusions [4].

This article introduces a novel Bayesian modelling approach to fMRI data analysis. Bayesian approaches have great potential in applications because they allow a flexible modelling of spatial and temporal correlations in data [5]. We consider a Bayesian spatiotemporal model in a computationally feasible manner to detect brain regions that are activated by the external stimulus in fMRI. Accurate and powerful single-subject task-related activation models are required in order to develop effective imaging biomarkers, which constitutes the primary scientific problem. Moreover, the Bayesian paradigm provides an attractive inferential framework that can directly incorporate the physical characteristics of an experiment. Our goal is to put forward a model and inferential framework by which to investigate task-specific changes in the BOLD signal. Clearly, the model must account for the spatial relationships between the voxels, but there are other possible sources of variation that should not be ignored. Ultimately, we develop a hierarchical Bayesian model, which not only takes into account the spatiotemporal and temporal drift relationships in the data under consideration, but also easily investigates the role of specific regions that integrate brain activity to coordinate cognition and behaviour.

The applied Bayesian hierarchical approach contains several characteristics [4, 6]. Latent binary variables are introduced to indicate activation/inactivation of voxels. The spatial generalized linear mixed model (SGLMM) [7, 8] is considered to capture the spatial dependence of the latent binary variables. In addition, the autoregression (AR) model is used to model the temporal dependence of signal changes. Several studies have found that spatial dependence also appears in the temporal parameters in AR models [4, 9]. Thus, neglecting the spatial dependence of temporal correlations moderates the computational intensity, but the simplification may produce a biased estimation of the temporal coefficients and consequently may result in the spurious detection of brain activities [10]. Therefore, we also consider the spatial linear mixed-effects model to spatially regularize the AR parameters.

In fMRI data, the posterior inferences are based on the estimations of parameters; however, the posterior distribution is typically extremely large, and is unavailable in analytic form. Hence, we employ Markov chain Monte Carlo (MCMC) [11] sampling techniques that combine Metropolis-Hastings [12] schemes to generate samples from the posterior distribution for the purpose of performing inferential tasks. To reduce the computational burden and accelerate the sampling procedure, we parcel the brain so that sampling procedures can be performed in parallel. In addition, prior information on activation in the form of spatially informative variables, e.g. the grey and white areas of the brain, can be incorporated into the model.

2. Statistical modelling

In this section, we introduce the proposed statistical model in the analysis of fMRI data. To make inferences about task-related change in underlying neuronal activity, a general linear model is used to model BOLD signal changes for each voxel. In addition to the essential temporal dependence of BOLD signals in a voxel itself, the BOLD signal changes show spatially contiguous and locally homogenous among voxels [11]. Shmuel et al. [13] in the visual fMRI study have demonstrated BOLD response seems to be well approximated by separable spatiotemporal model. Thus, for computational convenience, we consider a Bayesian separable spatiotemporal model for BODL signal changes to simultaneously account for the temporal dependence in nearby time points and spatial dependence in local neighbouring voxels.

Let yv=(yv1,,yvT), be a T × 1 column vector and denote the observed BOLD signals from a voxel v=1,,V at time t, t=1,,T. Following [14], we model the BOLD response for a particular voxel v with a linear regression model defined as

yv=Xvβv+Lvρv+εv;  εvNT(0, σ2IT),

where Xv is the design matrix, each column of which consists of values obtained from an impulse stimulus function [15, 16] with respect to a task convolved with the hemodynamic response function (HRF); βv, is a p×1 vector and corresponds to the effect of stimuli on the BOLD signal changes. The temporal correlation is modelled by ρv, being an r × 1 autoregression coefficient vector, with Lv, a T×r matrix of lagged prediction errors [2, 14]. We assume that the error terms media/ueq1.png in (1) are independently normally distributed NT(0,σ2IT) with a mean vector 0 of length T and a covariance matrix σ2IT across voxels, where IT is a T×T identity matrix [14]. In the Bayesian framework, the parameters are hierarchically assigned and the corresponding priors are defined, including spatial prior being introduced to capture the spatial dependence of brain activities among voxels.

For the purpose of detecting the activation of a voxel, a vector of binary random variables γv=(γv1,,γvp) is introduced to indicate whether the voxel v is in response to a sequence of input stimuli. The voxel v is considered active to the stimulus j if γvj=1 and, on the other hand, inactive if γvj=0. Given that γvj, we assume βvj has a spike and slab mixture prior of two normal distributions [17] given by


where N(μ,σ2) is denoted as a normal distribution with a mean μ and a variance σ2. In Eq. (2), we let cvj2 being fixed and assume τvj2 to have an inverse gamma distribution media/ueq2.png [4]. We consider that media/uequ2.png if the probability density function of X is defined by just below Eq. (3)


where Γ(a)=0ya1eydy. We set cvj2 to be large resulting in the nonzero estimate of βvj so that the stimulus j is considered to activate the voxel v. As George et al. [17] suggested, cvj2 should be taken less than 104 to avoid the computational problems and we find that cvj2=10 is a reasonable choice in our simulation studies and real fMRI example.

Since the response at a particular voxel is likely to be consistent to the responses of neighbouring voxels, we apply SGLMM [7, 35] to capture the spatial relationship. We assume γj=(γ1j,,γVj) are independently distributed in accordance to the Bernoulli distribution γvj|ηvj~Ber(ηvj) with a logistic link logit(ηvj)=αvj+mvϕj, where logit(ηvj)=ln(ηvj1ηvj) and αvj is a constant intended to incorporate the expert knowledge or anatomical information and ϕj is a vector of spatial random effects. The spatial dependence between the binary variables is implicitly captured by ϕj assuming to have


where mv is vth row of M, an V×q matrix consisting of multi-resolutional spatial basis vectors that are able to explain spatial variation sources. The columns of M, consist of the q principal eigenvectors with respect to the first q largest eigenvalues from the adjacency matrix A of voxels, an V×V matrix with the (v,s)th element in Eq. (5), defined as


where v~s indicates s and v are neighbours. In the three-dimensional fMRI data analysis, a 26-adjacent neighbourhood is considered [8]. Typically, q is less than V/2 or equal to the number of eigenvalues greater than 0.05. This choice can reduce the dimensionality significantly but still maintains the spatial structure of the data. The graph Laplacian matrix Q =diag(A1)A relates spatial basis vectors to represent the image data [18] and 1 is a q×1 column vector of ones. We assume a conjugate prior for the smoothing parameter given by


To avoid generating artefactual spatial structure in the posterior distribution, [19] suggested aκ=1 and bκ=4000 in Eq. (6). We refer media/uequ2.png if the probability density function of X is defined in just below Eq. (7)


It is also found that the temporal correlations between voxels tend similar [4]. Such spatial dependence is modelled by a spatial linear mixed-effects model. For computational convenience and simplicity, we make a transformation of ρvr such that ρ˜vr=log1+ρvr1ρvr and assume ρ˜vr~N(mvφr,λr2). For the spatial random effect for the temporal parameter ρ˜s is assumed to be


where φr is the spatial random effect for the rth order of the temporal parameters, and M and Q in Eq. (8) are the same as those listed in Eq. (4). Finally, we assume that media/ueq4.png, media/ueq5.png and media/ueq6.png. The values of a’s and b’s in the gamma or inverse gamma distributions are determined by the user to reflect the strength of one’s prior belief before observing the data.

3. Posterior inference

To explore parallel computation, which allows us to substantially speed-up expensive operations in the MCMC iteration, we partition the brain into non-overlapping areas such as rectangular three-dimensional lattices or Brodmann areas. We then carry out the statistical model in Section 2 to analyse the partitioned data. In this chapter, for simplicity, we use Brodmann’s map to parcel the brain into distinct regions before implementing the proposed model to the real fMRI data. In addition, separate regions divided by a data-drive segmentation procedure using functional clustering can be considered. Next, we introduce the likelihood function and the priors combining the parcellation procedure to form the posterior distribution for inference.

Suppose that a brain can be divided into G parcels, where the gth parcel contains Vg voxels, g=1,,G and where we denote the voxel-level parameters as θv=(βv,ρv,γv,σv2) and the parcel-level parameters as Θg=(φ,κ,τ2,λ2,ϕ,ω). The posterior distribution is obtained by combining the priors π(θv,Θg) and the likelihood L(θv|yv), which are defined in Eqs. (9) and (10) as follows [8, 20]:




Due to the intractability of the posterior, Gibbs and Metropolis-Hastings updates [20, 34] are applied to sample the posterior distribution for estimation and inference of the model parameters.

The posterior quantities of interest are P(γvj=1|y) for the activation map and E(βvj|y) for the magnitude of the effect caused by the stimulus j on the BOLD signal changes for voxel v. They can be directly estimated based on MCMC samples by


where γvj(m) and βvj(m) are the mth sample from a total M MCMC samples for γvj and βvj, respectively. The distribution of P^(γvj=1|y) in a map provides a way to visually inspect brain regions with peak, high, low and practically no activation. In addition, E^(βvj|y) in Eq. (11) offers researchers to view the strength of response in the brain to the stimulus.

The construction of the binary activation map is obtained by thresholding the posterior probability of γvj=1. A threshold value c needs to be used in the detection of brain activity, that is, a voxel v is defined as active to a stimulus j if P^(γvj=1|y)>c. However, threshold determination lacks agreement among investigators. Under certain conditions, a median selection criterion, c=0.5, results in the minimum prediction risk [21]. Smith et al. [6] and Lee et al. [4] defined a threshold by matching a Bayes factor approximation to a likelihood ratio test for activation such that c=0.8722 at a 5% level of significance. Additional to both, Kalus et al. [22] applied the false discovery rate (FDR) [23] to determine a threshold. In practice, it would be reasonable to construct activation maps for a grid of thresholds and to evaluate elicited results with neuroscience experts. In our experience, activation maps seem to be comparably robust against an exact choice of the threshold.

4. Simulation studies

We conduct stimulation studies to investigate the performance of the proposed model on the detection of brain activity. We measure the accuracy of the classification of voxels as either active/inactive and show the proposed model to be robust with regard to different spatial structures. Finally, we illustrate the implementation of the proposed model to a simulated fMRI data that mimics a study in face repetition effects in memory tests [24].

4.1. Benchmark example

Consider a 30×30 binary activation image γv, v1,,900, generated independently from


where logit(ηv)=αv+mvϕ and


We let αv=0, κ=0.5, where M is a 300×900 matrix whose columns are q=300 principal eigenvectors of the adjacency matrix, A, for the image corresponding to the first q=300 largest eigenvalues. Moreover, mv is the vth row of M, and Q=diag(A1)A.

We consider an AR(1) temporal correlation between time points within voxels. For each voxel v, let ρ˜v=log1+ρv1ρv and ρ˜v are generated from


We let λ2=0.1 and ω=2.

We consider a stimulus and the block and event-related designs with a total number of time points spanning 400 and repetition time (TR) equal to 2s. In the block design, the duration time is 20s. In the event-related design, the stimulus is randomly assigned. The stimulus function is convolved with a HRF modelled by a double gamma function [15] to create the design matrix Xv. The stimulus functions and the corresponding design matrices for the block and event-related designs are shown in Figure 1(a) and (b), respectively.


Figure 1.

The stimulus functions and the predicted BOLD signal changes used in the Benchmark example.

Given Xv,γv,ρv, the BOLD signal changes yv are generated from a normal distribution with a covariance σv2Λv, where σv2=1 and the (u,v)th element of Λv is ρ|uv|, and a mean of Xvβv, where


We then apply the proposed model to detect the activation areas, that is, to identify which voxel has γv=1. We ran MCMC chains with 100,000 iterations in order to ensure the largest Monte Carlo standard error (MCSE) [25] of all posterior probabilities of γv=1, less than 0.01.

We obtained posterior activation maps by setting the posterior probability threshold at 0.8722, that is, the voxel is categorized as active if P^(γv=1|y)>0.8722, and it is categorized as inactive otherwise [26]. To measure the performance of our proposed model, we calculate the true classification rate (TCR), the true positive rate (TPR) and the false positive rate (FPR), which are defined as

TCR=number of correctly classified voxelsnumber of voxels;
TPR=number of active voxels correctly claimed as activenumber of active voxels;
FPR=number of inactive voxels correctly claimed as activenumber of inactive voxels.

Table 1 shows the proposed method performs well with regard to detecting the active voxels where only a small number of active voxels are falsely identified as inactive in both designs. Over 10 replications, the model with consideration of spatial dependence performs better on the detection of activations. Especially, FPR reduces around 50%.

Block design (%)Event-related design (%)
TCR 92.44 (90.28–96.67)91.12 (91.82–96.00)
TPR99.97 (99.94–100)99.97 (99.94–100)
FPR6.98 (3.45–8.75)6.54 (3.54–8.66)
Block design (%)Event-related design (%)
TCR96.33 (95.89–97.56)95.78 (95.12–97.50)
TPR99.99 (99.97–100)99.99 (99.97–100)
FPR3.42 (1.12–5.27)3.22 (1.32–4.78)

Table 1.

The percentage of correct classification of voxels with a threshold at 0.8722 for considering the spatial dependence in (a) and without considering spatial dependence in (b) for temporal parameters in the model.

[i] - The values within the parentheses are the range of different rates over 10 replications.

4.2. Structures with spatial dependence

To demonstrate that the proposed model is able to handle the different binary and temporal image spatial dependencies, we generate a dataset that is the same as the benchmark example in the paper of [27]. We create a 20×20 binary image as shown in Figure 2(a), where the areas in red indicate that the voxels are active, and otherwise, they are not. Moreover, the values of ρs are generated from a uniform distribution between –0.3 and 0.3 for the non-active areas. However, we assign the fixed values –0.5, 0.5, and 0.75 to ρ in each active area, respectively, as shown in Figure 2(b). This is designed to investigate the patterns of temporal dependencies within the different areas. Similar to the settings discussed in Section 3.1, we consider both block and event-related designs with one stimulus. First, βv’s generated from U(2, 5) and σv2 from U(1, 3). The design matrix Xv and temporal correlation matrix Λv are defined to be the same as in Section 3.1. Given Xv,γv,ρv,βv, and σv2 the BOLD signal changes are generated from a normal distribution with a mean Xvβv and a covariance σv2Λv. We carry out the proposed approach to detect the activation and estimate the parameters of interest. For ease of visualization of the results, we provide one of the estimated binary images without thresholding and the corresponding estimated image, as shown in Figure 2(c). By visually inspecting estimated values of ρs in Figure 2(d) compared to the simulated ones in Figure 2(b), we are confident that the spatial dependences also can be captured by the model.

After thresholding by a value determined by controlling FDR as 0.05, the accuracy of the classification of voxels measured by TCR is 99.25%, and the FPR is 4.84% with consideration of the spatial dependence of the temporal correlations over 10 replications. On the other hand, when the threshold is taken to be equal to 0.8722, TCR is 99.17% and FPR is 5.01%. Therefore, the proposed model can be applied to detect the activation of brain image data even for different spatial dependence structures.


Figure 2.

Simulated and estimated binary and ρ images.

4.3. Parcellation effect

In the previous simulations, image data was analysed together. In this simulation, we partition the data and then apply the proposed model to analyse each unit of the partitioned data. Our goal is to investigate the effect of the parcellation of an image on the estimation of active voxels. Given the binary image γ shown in Figure 3(a), a unit of data is generated with the parameters being the same as the settings discussed in Section 3.1.

We consider four different parcellations, as shown in Figure 3(a), (b) and (d). Table 2 shows the values of different measurements over 10 replications. Little difference can be found on the measures between different parcellations. However, the computational time is dramatically reduced for case (d), where each parcellation contains fewer voxels. The same results in block and event-related fMRI experiments.


Figure 3.

Four different parcellation schemes. The areas in red are active areas corresponding to γv=1. The blue dashed line is used to partition the area into non-overlapping parts.

Parcellation Scheme(a)(b)(c)(d)
TCR99.12 (98.61–99.84) %99.07 (98.68–99.56) %98.94 (98.03–99.64) %99.08 (98.76–99.51) %
FPR6.35 (1.51–9.97) %6.68 (3.67–9.36) %7.45 (3.12–13.87) %5.82 (2.83–9.44) %
Running time ratio8.752.752.751

Table 2.

The average percentage of correctly classified voxels (TCR), true positives (TPR), and false positives (FPR) over 10 replications for different image parcellations corresponding to Figure 3(a)(d).

[i] - The values within the parentheses are the ranges of different rates over 10 replications.

4.4. A real-life simulation example

To further demonstrate the capability of the proposed model to accurately identify the activation regions, we simulate an fMRI dataset from an R package: neuRosim [28]. The simulated data follows the example of event-related fMRI study [24]. We consider four different tasks denoted by N1, N2, F1 and F2, which are presented randomly in the experiment. The onsets for each condition are schematically shown in Figure 4. We then specify three areas that are activated by the tasks shown in Figure 5. The simulated four-dimensional fMRI data consists of 351 scans, each containing three-dimensional image of size 53×63×46 and being collected in every 2 s resulting the total time for the experiment is 11 min and 42 s.


Figure 4.

Depiction of temporal occurrences of four different tasks.


Figure 5.

Activated regions upon representation of different tasks.

Now we apply the proposed model to analyse the simulated data. We split the array into disjointed arrays based on the Brodmann level regions. The voxel for all regions ranges from 236 to 969. We collect 1,000,000 samples to estimate the posterior probability of activation in each voxel. We consider the following probability:


which is used to denote whether the voxel v is activated by any of applied stimuli. Once the posterior probability is greater than 0.8722, then the corresponding voxel is considered as activated by a stimulus. The activation maps registered into MNI 152 template (a reference brain map from the Montreal Neurological Institute) and the areas considered active are in yellow as shown in Figure 6. In addition to comparing specified activation regions in Figure 5 to those detected in Figure 6, we also calculate the three measures, TCR, TPR and FPR which are 97.71%, 99.34% and 3.25%, respectively. The results show that proposed model performs well on classifying the active and inactive voxels.


Figure 6.

Estimated probability of activation in response to any face task with activation areas in yellow.

5. Application

We apply the spatial Bayesian variable selection approach discussed in Section 2 to real fMRI data, a study of the Simon effect [29]. In this Simon task, participants responded to one colour with the right hand, and to the other colour with the left hand. Congruence in this aspect means that a colour appears on the default side (congruent condition), or it may appear on the opposite side (incongruent condition). One participant’s brain image data is selected to be analysed for the purpose of illustration. The four-dimensional fMRI BOLD image data was pre-processed using FSL from Oxford Centre for Functional MRI of the Brain [30], including motion correction, realignment, slice timing correction, spatial smoothing and high-pass filtering, before being analysed using our model. Incorrect trials were removed. The four tasks in this experiment were convolved with a double gamma function. We divided the pre-processed data into 48 Brodmann areas. The statistics of interest were registered into a standard template MNI152 for ease of visualization. We used a first-order temporal autocorrelation throughout the study.

We consider a voxel to be active when the posterior probability of γv is greater than c, which is either a deterministic value of 0.8722 [4] or is determined by an FDR equal to 0.05 [22]. However, in this event-related fMRI, 0.8722 seems too conservative, so the value of c is determined corresponding to an FDR equal to 0.05. The activation maps for different tasks are shown in Figure 7(a)(d).

One interesting question is to compare active domains corresponding to congruency. Posterior probabilities helped us find the following inequality between the incongruent and congruent tasks groups:


where βs is the magnitude of corresponding effect on the BOLD signal changes. The posterior activations are shown in Figures 8 and 9 when posterior probabilities are thresholded by a critical value corresponding to FDR equal to 0.05. The active regions include frontal and occipital lobes detected by the models with and without consideration of spatial dependence of temporal correlations. These results are consistent with other studies [31, 32]. However, more active areas were detected by the model without consideration of spatial dependence of temporal correlations. For something like the Simon task, it is expected that the lateral prefrontal cortex to be activated as they are parts of the cognitive control network/multiple demand network, while the medial prefrontal cortex is part of the default mode network [33]. Therefore, it is possible that the activation in the medial part may be some spread of activity and therefore false positives. This is an important example that illustrates considering the spatial dependence of temporal correlations may improve the estimation and reduce the some spurious detection of activation.


Figure 7.

The activation areas with a posterior probability greater than c, corresponding to FDA = 0.05, shown in red for different tasks.


Figure 8.

Areas shown more activity in incongruent than congruent task without consideration of spatial dependence of temporal correlations in the model.


Figure 9.

Areas shown more activity in incongruent than congruent task with consideration of spatial dependence of temporal correlations in the model.

6. Discussion and conclusions

In this work, we applied a new approach to performing Bayesian variable selection with consideration of spatial dependencies from both regression and temporal coefficients in single-subject event-related fMRI data. Through simulations, improvement in the detection of brain activity was observed for a wide range of different spatial structures, parcellations and experimental designs. We found that the proposed approach potentially decreases false positive rates as shown in Table 1 in the simulation and Figure 9 in the real example. In addition, prior information from brain structure and function for subject-level inference can be incorporated into the analysis.

It is worth noting that the proposed Bayesian approach partitions a brain into several regions before implementing the model to the fMRI data. In addition to Brodmann areas, several parcellations of the brain into distinct regions are available, such as the separate regions divided by a data-drive segmentation procedure using functional clustering. From the simulation study, we found the computational burden to be greatly reduced by the joint use of parcellation and an SGLMM. In particular, the probability of activation and activation magnitudes were readily computed without requiring an adjustment for multiple comparisons in a post-processing step. For broader goals, it would be of interest to extend the model to group studies. In addition, we are confident that Bayesian approaches represent an important direction in fMRI and in high-dimensional research in general.


This research was, in part, supported by the Ministry of Education, Taiwan, R.O.C. and National Cheng Kung University (the Aim for the Top University Project to National Cheng Kung University, NCKU). We thank the Mind Research and Imaging Center (MRIC) at NCKU for consultation and instrument availability.


1 - J. A. Detre and T. Floyd. Functional MRI and its applications to the clinical neurosciences. The Neuroscientist. 2001;7:64–79.
2 - W. D. Penny, K. J. Friston, J. T. Ashburner, S. J. Kiebel, and T. E. Nichols, editors. Statistical Parametric Mapping: The Analysis of Functional Brain Images. London: Academic Press; 2006.
3 - K. J. Worsley, S. Marrett, P. Neelin, A. C. Vandal, K. J. Friston, and A. C. Evans. A unified statistical approach for determining significant voxels in images of cerebral activation. Human Brain Mapping. 1996;4:58–73.
4 - K. Lee, G. L. Jones, B. Caffo, and S. Bassett. Spatial Bayesian variable selection models on functional magnetic resonance imaging time-series data. Bayesian Analysis. 2014;9:699–732.
5 - L. Zhang, M. Guindani, and M. Vannucci. Bayesian models for functional magnetic resonance imaging data analysis. Wiley Interdisciplinary Reviews: Computational Statistics. 2015;7:21–41.
6 - M. Smith and L. Fahrmeir. Spatial Bayesian variable selection with application to functional magnetic resonance imaging. Journal of the American Statistical Association. 2007;102:417–431.
7 - J. Hughes and M. Haran. Dimension reduction and alleviation of confounding for spatial generalized linear mixed models. Journal of the Royal Statistical Society Series B. 2013;75:139–159.
8 - D. R. Musgrove, J. Hughes, and L. E. Eberly. Fast, fully Bayesian spatiotemporal inference for fMRI data. Biostatistics. 2016;17:291–303.
9 - M. Woolrich, M. Jenkinson, J. Brady, and S. Smith. Fully Bayesian spatio-temporal modeling for fMRI data. IEEE Transactions on Medical Imaging. 2004;23:213–231.
10 - A. Eklund, T. Nichols, M. Andersson, and H. Knutsson. Empirically investigating the statistical validity of SPM, FSL and AFNI for single subject fMRI analysis. In: 2015 IEEE 12th International Symposium on Biomedical Imaging (ISBI); New York, NY; 2015; pp. 1376–1380.
11 - A. Quirós, R. M. Diez, and D. Gamerman. Bayesian spatiotemporal model of fMRI data. NeuroImage. 2010;49:442–456.
12 - W.K. Hastings. Monte Carlo sampling methods using Markov chains and their applications. Biometrika. 1970;57:97–109.
13 - S. Brooks, A. Gelman, G. Jones, and X. Meng. Handbook of Markov Chain Monte Carlo. Chapman and Hall/CRC; 2011 Taylor & Francis Group, 6000 Broken Sound Parkway NW, Suite 300 Boca Raton, FL 33487–2742.
14 - W. D. Penny, N. J. Trujillo-Barreto, and K. J. Friston. Bayesian fMRI time series analysis with spatial priors. NeuroImage. 2005;24:350–362.
15 - M. A. Lindquist, J. M. Loh, L. Y. Atlas, and T. D. Wager. Modeling the hemodynamic response function in fMRI efficiency, bias and mis-modeling. NeuroImage. 2009;45:S187–S198.
16 - J. C. Rajapakse, F. Kruggel, and J. M. M. D. Y. Cramon. Modeling hemodynamic response for analysis of functional MRI time-series. Human Brain Mapping. 1998;6:283–300.
17 - E. I. George and R. E. McCulloch. Approaches for Bayesian variable selection. Statistica Sinica. 1997;7:339–374.
18 - D. M. Cvetkovć, M. Doob, and H. Sachs. Spectra of Graphs: Theory and Applications. 3rd ed. New York, NY: Wiley; 1998.
19 - N. G. Best, L. Waller, A. Thomas, E. M. Conlon, and R. Arnold. Discussion of Bayesian models for spatially correlated disease and exposure data. In: J. O. Berger, A. P. Dawid, and A. F. M. Smith, editors. Bayesian Statistics 6. Oxford University Press; 1999. pp. 147–150.
20 - A. Shmuel, E. Yacoub, D. Chaimow, N. K. Logothetis, and K. Ugurbil. Spatio-temporal point-spread function of fMRI signal in human gray matter at 7 Tesla. NeuroImage. Oxford 2007;35:539–552.
21 - M. Barbieri and J. O. Berger. Optimal predictive model selection. Annals of Statistics. 2004;32:870–897.
22 - S. Kalus, P. G. Sämann, and L. Fahrmeir. Classification of brain activation via spatial Bayesian variable selection in fMRI regression. Advances in Data Analysis and Classification. 2014;8:63–83.
23 - Y. Benjamini and Y. Hochberg. Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society, Series B. 1995;57:289–300.
24 - R. Henson, T. Shallice, M. Gorno-Tempini, and R. Dolan. Face repetition effects in implicit and explicit memory tests as measured by fMRI. Cerebral Cortex. 2002;12:178–186.
25 - G. Jones, M. Haran, B. Caffo, and R. Neath. Fixed-width output analysis for Markov chain Monte Carlo. Journal of the American Statistical Association. 2006;101:1537–1547.
26 - O. Josephs and R. Henson. Event-related fMRI: modelling, inference and optimisation. The Royal Society: Philosophical Transactions B. 1999;354:1215–1228.
27 - M. Bezener, J. Hughes, and G. L. Jones. Bayesian spatiotemporal modeling using hierarchical spatial priors with applications to functional magnetic resonance imaging. Available from:
28 - M. Welvaert, J. Durnez, B. Moerkerke, G. Berdoolaege, and Y. Rosseel. neuRosim: An R package for generating fMRI data. Journal of Statistical Software. 2011;44:1–18.
29 - T. Wen and S. Hsieh. Neuroimaging of the joint Simon effect with believed biological and non-biological co-actors. Frontiers in Human Neuroscience. 2015;9:1–13.
30 - FMRIB, Oxford, UK. FMRIB Software Library [Internet]. Available from:
31 - X. Liu, M. Banich, M. Jacobson, and J. Tanabe. Common and distinct neural substrates of attentional control in an integrated Simon and spatial Stroop task as assessed by event-related fMRI. NeuroImage. 2004;22:1097–1106.
32 - B.U. Forstmann, W. P. M van den Wildenberg, and K. R. Riderinkhof. Neural mechanisms, temporal dynamics, and individual differences in interference control. Neuroscience. 2008;20:1854–1865.
33 - T. A. Niendam, A. R. Laird, K. L. Ray, Y. M. Dean, D. C. Glahn, and C. S. Carter. Meta-analytic evidence for a superordinate cognitive control network subserving diverse executive functions. Cognitive, Affective, and Behavioural Neuroscience. 2012;12:241–268.
34 - A. Gelman, J. C. Carlin, H. Stern, D. Dunson, A. Vehtari, and D. Rubin. Bayesian Data Analysis. 3rd ed. Chapman and Hall/CRC; 2013.
35 - C. Berrett and C. A. Calder. Bayesian spatial binary classification. Spatial Statistics. 2016;16:72–102.