# [

###### Abstract

Statistical modeling of fMRI data is challenging as the data are both spatially and temporally correlated. Spatially, measurements are taken at thousands of contiguous regions, called voxels, and temporally measurements are taken at hundreds of time points at each voxel. Recent advances in Bayesian hierarchical modeling have addressed the challenges of spatiotemproal structure in fMRI data with models incorporating both spatial and temporal priors for signal and noise. While there has been extensive research on modeling the fMRI signal (i. e., the covolution of the experimental design with the functional choice for the hemodynamic response function) and its spatial variability, less attention has been paid to realistic modeling of the temporal dependence that typically exists within the fMRI noise, where a low order autoregressive process is typically adopted. Furthermore, the AR order is held constant across voxels (e.g. AR(1) at each voxel). Motivated by an event-related fMRI experiment, we propose a novel hierarchical Bayesian model with automatic selection of the autoregressive orders of the noise process that vary spatially over the brain. With simulation studies we show that our model has improved accuracy and apply it to our motivating example.

SVARO]Bayesian Analysis of fMRI data with

Spatially-Varying Autoregressive Orders
Teng, M. et al.]Ming Teng
]Farouk S. Nathoo
]Timothy D. Johnson

## 1 Introduction

In the analysis of functional magnetic resonance imaging (fMRI) data a key challenge is dealing with spatial and temporal correlation. The temporal correlation can arise from many sources, including scanner drift at very low frequencies, slow vascular/metabolic oscillations that are typically of moderate to low frequency, and some other sources of noise such as breathing and heartbeat. Simply ignoring these sources of autocorrelation may lead to increased false positive discoveries (Makni et al., 2006). To deal with these issues, a variety of approaches have been proposed. One commonly used approach, namely “prewhitening”, works by estimating the temporal autocorrelation and then de-correlating the noise using the estimates (Bullmore et al., 1996; Locascio et al., 1997). Besides these stationary time series models, non-stationary models have also been proposed (Zarahn et al., 1997; Bullmore et al., 2004). According to Friston et al. (2000), prewhitening can produce an extraneous source of bias. Alternatively, a band-pass filtering technique known as “pre-coloring” can be applied to the data first, followed by statistical modeling that deals with the autocorrelation in the colored data. For a review and discussion of these approaches the reader is referred to Woolrich et al. (2001). While high-pass filtering has proven to be beneficial in increasing the power of the statistical analysis, the low-pass filtering involved in coloring is considered controversial in that it tends to add autocorrelation into the data (Skudlarski et al., 1999; Della-Maggiore et al., 2002).

While accurate temporal modeling is important for estimation of the fMRI signal and its standard error, traditional approaches apply a temporal model at each voxel, independently. That is, they ignore spatial correlation. More specifically, this mass univariate approach, considered to be the classical approach to the analysis of fMRI data, includes a smoothing step involving a spatial Gaussian filter that is applied to the data first (Friston et al., 1995), followed by model estimation at each voxel, and then statistical inference is based on random field theory (Worsley and Friston, 1995) which is applied to adjust for multiplicity in the spatial domain. While this approach remains the most common approach for analyzing fMRI data it has been criticized on a number of grounds. For example, the Gaussian kernel that is used to smooth the data has to be pre-specified and introduces artificial correlation into the data. In addition, this approach does not directly account for spatial correlation in the model.

Partly as a result of these criticisms, Bayesian models with spatial structured priors have been proposed that allow for the calculation of posterior probability maps (PPM) for activation. This Bayesian approach to inference is based on an explicit spatial model and does not require smoothing the data with a Gaussian kernel nor does it require the use of random-field theory-based adjustments for multiplicity. A variety of spatial-temporal Bayesian models have been proposed. One model that is widely used and implemented within the SPM software is the GLM-AR (general linear model, auto-regressive) model (Penny et al., 2003, 2005, 2007), The GLM-AR model assumes that the data can be decomposed into two sources of variability. The first source is the product of the design matrix for the fMRI experiment convolved with a hemodynamic response function (HRF) and experimental factors, and the second source represents temporally correlated noise that is modeled using a low-order AR structure. In addition, the regression coefficients and the autoregressive coefficients vary across space and are assigned spatial smoothing priors. Gössl et al. (2001) proposed a model where the data are decomposed into three sources, a spatial stimulus, a deterministic trend, and a white noise process. However, this modeling approach may not account for some higher frequency stochastic noise components. Woolrich et al. (2004b) assumed that the temporal noise arises from both large scale and small scale variation, and built a space-time simultaneously auto-regressive model that accounts for both scales of variation. Methods focusing on spatial variable selection have also been proposed (see, e.g., Bezener et al. (2016), Lee et al. (2014),Musgrove et al. (2016)); while Kim et al. (2010) proposed a mixture of experts model to represent spatial activation clusters. While these models have a number of different characteristics which make the approaches unique, most of them commonly assume a homogeneous, low order AR or ARMA process for the temporal noise. By homogeneous, we mean that the order of the AR or ARMA process is assumed constant across all voxels. This assumption is also made in Penny et al. (2003); however, as we demonstrate using a simple empirical example in the next section, this homogeneous AR order assumption may be violated in real fMRI data.

Instead of formulating the model at each voxel and then adopting spatial smoothing priors for parameters across voxels, another branch of research is based on vector autoregressive (VAR) processes, see Harrison et al. (2003). This approach allows for time-lagged dependence across voxels and spatial-temporal interaction but fitting these models across a large number of voxels is computationally intractable and low-rank approximations have to be used. These models are also useful for studying effective brain connectivity, where one time course is used to predict the other (Castruccio et al., 2016; Chang and Glover, 2010). Another line of work chooses to model the temporal noise as a long memory process, and applies discrete wavelet transforms (DWT) towards fitting the model (see, e.g., Jeong et al. (2013); Bullmore et al. (2004); Fadili and Bullmore (2002); Meyer (2003)). While this approach seems promising, our focus in this paper will be with modeling short term memory using the classical AR process and spatial priors. The reason we choose to work with the AR process is because of its mathematical amenability and simplicity, and its wide use in different areas of science. A novel aspect of our work is that we allow the data to determine the order of the AR process at each voxel, borrowing strength from neighboring voxels, using ideas from Bayesian spatial variable selection.

Computation is an important issue when considering Bayesian spatial-temporal models for fMRI data. While the main focus of this paper lies with the development of a new model, another aspect of this work is the comparison of fully Bayesian and approximate Bayesian computation methods. Due to the computational burden associated with fitting models to high-dimensional brain imaging data, approximate Bayesian methods have received considerable attention in the neuroimaging literature. One such method is the variational Bayes (VB) inference (Penny et al., 2003, 2007; Woolrich et al., 2004a). As there are currently no theoretical results quantifying the accuracy of VB methods (in contrast to MCMC which is justified by large sample theory of stationary Markov chains), the evaluation of VB has to be performed on a case-by-case basis. In some cases, the performance of VB can be quite good and in other cases it can be quite poor. In addition to the implementation of our new model based on a suitably designed MCMC sampler, we also develop an MCMC algorithm to sample from the posterior of the original GLM-AR model. We then compare our model with both the VB implementation of the GLM-AR model (using SPM code) and our MCMC implementation of the GLM-AR model. Our studies indicate that under a low signal-to-noise (SNR) ratio the accuracy of MCMC outperforms VB according to several criteria.

### 1.1 Motivating example

Our motivating example comes from a single subject in a fMRI experiment examining a face-repetition stimulus. The experiment involves the presentation of either famous faces (F) or non-famous faces (N) with each type of face presented two times. Convolving this experiment design with the canonical hemodynamic response function and its time and dispersion derivatives leads to a design matrix with twelve columns plus one extra column for an intercept term in the regression model. After performing the necessary pre-processing steps as described in Penny et al. (2005), we fit a simple linear regression at each of the voxels. After obtaining the residuals from each voxel-specific fit, we fit an AR process up to order for each voxel using the “ar” function in R. We then selected the optimal AR orders based on the AIC criterion. Figure 1 displays a pictoral representation of the results.

Figure 1 shows considerable variablity in the estimated AR order across voxels. While most of the estimated optimal AR orders are 4 or less, higher orders up to are selected at some of the voxels. Furthermore, these estimated AR orders tend to show some extent of spatial clustering. If, as is often done, we simply model the data using a homogeneous low-order AR process, then the voxels with higher order estimated AR orders would be incorrectly modeled, and this inaccuracy in the modeling of temporal noise will have an impact on the inference on the covariates of interest (via underestimated standard errors), resulting in potentially false inferences about brain activation. To address this issue, we propose a spatially varying autoregressive order (SVARO) model, where the AR orders vary spatially across the brain. This is made possible by adopting a spike-and-slab prior with a stochastic search variable selection scheme. Spatial clustering of AR orders is incorporated by imposing an Ising prior (Ising, 1925) as the latent indicator for the spike and slab prior. We update the latent indicators using the Swendsen-Wang algorithm (Swendsen and Wang, 1987) alternating with Gibbs sampling in our MCMC algorithm. To prevent the phase transition problem associated with the Ising model, we derive theoretical bounds as in Li et al. (2015) and use these bounds to prevent critical slowing of the algorithm. We compare our model with the GLM-AR model of Penny et al. (2007) (implemented under two schemes: our self written MCMC sampler and the VB algorithm available in the SPM software) in terms of mean squared error (MSE) and sensitivity. We conduct these comparisons using two simulation studies and then compare results on the motivating data set.

The rest of the paper is organized as follows: In Section 2 we present our model and MCMC sampling scheme. We present results from our simulation studies in Section 3, followed by the analysis of the face-repitition data set in Section 4. Lastly, we provide a discussion and outline some possible directions for future work in Section 5.

## 2 The Spatially-Varying Autoregressive Order Model

First, we will define the model likelihood and then specify the spatial and temporal priors. Next, we discuss our posterior sampling scheme, include the construction of bounds for the hyperpriors that we use for the Ising prior on the AR orders. Last, we discuss inference based on posterior probability maps.

### 2.1 Model likelihood

We let denote the maximum possible order of the AR process at each voxel and let denote the number of regression coefficients at each voxel. Using similar notation as in Teng et al. (2016), for voxel , we let denote the observed time series of length . For simplicity, our model is specified conditional on the first observations at each voxel so that the likelihood function is constructed based on the model for the remaining observations in the time series. We let denote the design matrix, denote the -dimensional vector of regression coefficients at voxel , and denotes the corresponding error term. Define the vector , the entire time series observed at voxel . The hierarchical model is specified in several stages. The first stage is a general linear model:

(0) |

where we emphasize again the implicit conditioning on . Let denote the embedded error (or lagged prediction) matrix of dimension , with element where the notation denotes the index of the vector. Let denotes a vector of i.i.d mean-zero Gaussian random variables with precision . The second stage is then an AR model at each voxel:

(0) |

where is a vector of autoregressive coefficients for the time series at voxel .

Letting denote a constant term, the log-likelihood for voxel , is

(0) |

where is the th row of the design matrix and is the th row of .

Summing this log-likelihood over the number of voxels, , we obtain the overall likelihood:

(0) |

### 2.2 Spatial prior

At the next level of the model we specify a spatial smoothing prior for the regression coefficients , a matrix, with th row . Following Penny et al. (2005), we assume that the prior for takes the form

(0) | |||||

(0) |

A priori the regression coefficients within voxels are independent (2.2) whilst spatially, the th, , regression coefficients (across voxels) are modeled dependently through an -dimensional multivariate normal distribution (2.2). Here is known as a Laplacian matrix. The th diagonal term of this matrix is equal to the corresponding number of first order neighbors of the voxel . All off-diagonal terms are zero except for in off-diagonal elements and if voxel is a neighbor of voxel , for . This form for the prior accommodates spatial smoothing while also being sparse and convenient to work with computationally. In the SPM12 software this prior is referred to as the “LORETA” prior. Ultimately, what is of primary interest in studies of brain activation is a posterior probability of some function of these regression coefficients, and this posterior probability is computed at each voxel to produce posterior probability maps (PPM, see Penny et al. (2005)). The precision parameter in (2.2), , is assigned a conditionally conjugate hyper-prior:

(0) |

### 2.3 Temporal prior

The key difference between our model and the model of Penny et al. (2007) lies in our modeling of the temporal noise. Rather than assuming AR orders are homogeneous across the brain (we refer the readers to Teng et al. (2016) and Penny et al. (2007) for model details), we allow for variability in the order of the AR processes across voxels. In addition, we adopt a spatial prior for this variability under the assumption that the AR orders of neighboring voxels will be similar. Specifically, for each voxel and order , , we assign the latent indicator variable to the AR coefficient , such that given , will be conditionally independent. will take value if order is present for voxel and otherwise. Conditional on , will either have a normal distribution or unit mass at . This is commonly referred to as the spike-and-slab prior (George and McCulloch, 1993; Mitchell and Beauchamp, 1988), though we note that our formulation is a spatial spike-and-slab prior and that this prior is assigned to the coefficients of the AR process governing the temporal noise.

Here, is the pdf of a normal distribution with mean and variance and is the indicator function that its argument equals 0, and where is the binary indicator. is the precision of the normal component and is again given a Gamma prior .

The advantages of introducing such a prior are three fold: First, the orders in the AR process at each voxel that lack support from the data can be effectively removed from the model as the corresponding AR coefficients can be shrunk exactly to . This allows us to infer which orders are present in which voxels. Second, the number of voxels with high AR orders is non-zero but expected to be small, which is an aspect of this prior that can be controlled by tuning the hyper-parameters. Third, for some of the voxels there might be vacancies in some of the middle orders while there are some non-zero coefficients for higher orders. The proposed model is flexible enough to allow for this behaviour, since we have a total of independent Ising processes, one for each possible order .

There are of course other model selection techniques that could have been considered. For example a type of Bayesian lasso could have bee used as an alternative to the spike-and-slab prior. Wang et al. (2007) have applied the lasso to the selection of AR processes, and for Bayesian lasso we refer to Schmidt and Makalic (2013). A recent alternative prior known as the ”non-local” prior for variable selection has been proposed by Johnson and Rossell (2012) and has been demonstrated to have desirable consistency properties and yield smaller prediction errors in large sample settings. A review of Bayesian priors that can be employed for model selection is presented in O’Hara et al. (2009).

We assume that the indicator processes are independent across different orders, , where . The simplest variable selection model would assume follows a Bernouli distribution (George and McCulloch, 1993). Here, in order to borrow information across neighbors as well as to model the spatial clustering effect of AR orders, we choose to use the Ising prior (Smith and Fahrmeir, 2007) independently for each .

(0) |

where and are two hyper-parameters controlling the sparsity and smoothness of the binary latent field, respectively. Typically, a higher value of results in less sparsity and a higher value of indicates more smoothness. One issue with the Ising model that requires some care is the choice of hyper-parameters. When these parameters take values near what is known as the “phase transition” boundary, the mixing of an MCMC sampler will suffer from critical slow down (Stanley et al., 1987). To avoid the phase transition boundary, we adopt an analytical approach similar to Li et al. (2015) to quantify the value for the bounds for and . An outline of the bound derivation is given in Subsection 2.5.

### 2.4 Posterior sampling scheme

Most parameter updates related to the posterior sampling of our model can be accomplished via Gibbs sampling. One exception is the update to the latent indicator . For we use the Swendsen-Wang algorithm alternating with Gibbs updates (Johnson et al., 2013). This strategy, that is, mixing Swendsen-Wang updates with Gibbs updates for has proven successful in improving the mixing of the Markov chain sampler and results in faster block updates in various studies (Higdon, 1998).

### 2.5 Bound construction for the Ising Model Hyper-Parameters

The hyper-parameters in the Ising prior play a vital role in posterior estimation. Without careful selection, we are faced with mixing issues associated with “phase-transition” (Stanley et al., 1987). There are various approaches to sampling such hyper-parameters, Johnson et al. (2013) estimated them using path sampling (Gelman and Meng, 1998), Shu et al. (2015) proposed a Monte Carlo EM algorithm to obtain a point estimate of the hyper-parameters, but these procedures would be too time consuming for our model, considering that we have over independent Ising fields and thousands of voxels. Smith and Fahrmeir (2007) proposed to update the hyper-parameters and binary indicators together, but this approach still suffers from potential possibility of sampling over the phase transition boundary. Here, we adopt an approach similar to that considered in Li et al. (2015) and construct some theoretical bounds to prevent phase transition. The resulting hyper-parameter values are then chosen as fixed values within the estimated bounds. This procedure turns out to work well in our analysis and studies.

To construct the bounds, we first write out the posterior conditional density with respect to ,

In our model where multiple orders exists across space, it is natural to assume that: 1) there are relatively few time series with large AR order; and 2) the posterior density when low AR orders exist is greater than that when low AR orders do not exist, meaning is greater than .

Let denote the candidate voxels selected for order . The maximum number of first order neighbors is 8. Let denote the length of an edge a voxel, then there are neighboring pairs. Based on this we derive

(0) |

According to 1), we know that for high AR orders (typically ) . According to 2), we know that for low AR orders (typically ),

(0) | |||

Reorganizing this by moving the first term on right-hand side to left produces:

(0) | |||

The two terms in the bracket on the left side can be considered as one with and without . Thus, it can be roughly considered as the residual sum of squares of a common linear regression when is included in the model or not. Let denote the coefficient of determination for voxel and order , then we have

(0) |

Combined with Equation ‣ 2.5, we have

(0) |

For a 3-dimensional grid we assume as the number of voxels. Among them, a proportion of are selected for order . So . We assume that 5% of the variation can be explained as a result of order , so . We then have .

Note that the inequality above just gives a range values for the hyperparameters, rather than providing the values directly. In practice, the exact values of hyperparameters are largely determined by the researcher, which should be combined with one’s prior experience and an initial analysis of the data. We suggest obtaining such values based on some exploratory ad-hoc approaches, e.g., a linear regression at each voxel followed by fitting an AR process. Then the the estimated optimal orders can be used as a reference when determining the hyper-parameters in the Ising model. This method has turned out to work well empirically as we demonstrate in Section 5.

### 2.6 Posterior probability maps

A primary emphasis on fMRI data analysis is inference for activation. For Bayesian modeling of fMRI data this is typically achieved via the PPM. Let denote a contrast of interest of the regression coefficients. A PPM is a map of the posterior probability of activation for each voxel: . Here is a pre-specified “activation threshold”, for example, a value that corresponds to of the global mean value. Thus, PPM looks at the probability of the contrast being greater than activation threshold , given the data.

To formally determine activation in the brain, one can look at a thresholded PPM. This is obtained by exerting a second threshold, namely a “probability threshold” , onto the original PPM. Thus, a voxel is “activated” if . This reflects the confidence of the inference and usually takes a value above (e.g. or ). This process discretizes the PPM into “null” and ”activated” voxels and is commonly used in summarizing a Bayesian analysis for brain activation.

## 3 Simulation Study

To evaluate the performance of our model, we make comparisons with the standard GLM-AR spatial model. One implementation of this model that we make comparisons to is the Variational Bayes (VB) method available in SPM12 software. Another implementation is our self-written MCMC sampler for the same model. Although the accuracy of VB has been verified in a setting with high signal-to-noise ratio (SNR) by Teng et al. (2016), under low SNR, MCMC outperforms VB according to certain metrics. This will be illustrated presently. Henceforth, we will refer to the VB implementation and MCMC implementation of the GLM-AR model as PVB and PMCMC, respectively.

### 3.1 Simulation design

Our design matrix consists of two columns (); the first column is the experimental design (fashioned after the face-repetition data set) convolved with the canonical HRF, and the second column is the intercept (a vector of s). The parameters of interest, corresponding to the experimental design, (one at each voxel) are generated under a mean zero multivariate normal distribution: , while the intercepts are generated under a mean multivariate normal distribution: . The model noise will have a precision of . This corresponds to a fairly small signal-to-noise (SNR) ratio, where the temporal noise will play a greater role in the data. In the following, we will carry out two simulations. In the first case the data will be simulated under our model, and in the second case the data will be generated under the standard spatial GLM-AR model. In these two simulations, we investigate the estimation accuracy of the slopes (), intercepts (), and autoregressive coefficients (), and we also examine if the difference in inference for these coefficients will lead to a possible difference in the final inference on brain activation. All simulations are based on replicate data sets, and we perform the simulations on a 2-dimensional axial slice of the brain.

### 3.2 Simulation 1

Here we simulate the AR parameters under our SVARO model, that is from the Ising prior. We assume that the maximum order at each voxel is . The precision parameters are set as . For simplicity, we assume that all AR orders are generated spatially according to the same values for the hyperparameters of the Ising model, i.e., and . The AR order in PMCMC and PVB are set to as is fairly standard practice. We note here that the GLM-AR model is misspecified and we expect that its performance will suffer regardless of which posterior sampling algorithm is used.

Figure 2 shows the true AR orders, the estimated maximum orders from our SVARO model, and the difference of the two. The estimated maximum orders are obtained by averaging the posterior mean at each voxel over their simulation replicates and rounding. We can see that most of the orders match between the two figures indicating good performance. However, there are some negative values in the difference map.

Next we compare SVARO with PVB and PMCMC in estimating the st AR coefficient. As shown in Figure 3, SVARO shows little error compared with the truth, indicating that our model has captured the autoregressive parameter quite well. In contrast, PMCMC and PVB exhibit more bias, indicating a lack of fit for the temporal noise. Note that we are only displaying the SVARO estimates for the st AR coefficient for simplicity and direct comparison, the other AR coefficients are similarly well-estimated.

Table 1 summarizes the average MSE for various parameters. These summaries are obtained by averaging the MSE of the corresponding parameters across all the voxels and over simulation replicates. It is clear that SVARO has the smallest MSE for these three parameters. In addition, PMCMC outperforms PVB in estimating the coefficient of the hemodynamic response, which is the primary parameter upon which inference is based. This finding is in line with our previous findings in Teng et al. (2016), where we found that a low SNR is one setting where MCMC outperforms VB for this particular model (GLM-AR). We also calculate and present the log-pseudo marginal likelihood (LPML, Gelfand and Dey (1994)) in Table 1. SVARO has a larger LPML than the MCMC implementation of GLM-AR, indicating a better model fit. Note that LPML cannot be obtained from the VB algorithm. In terms of timing, SVARO takes 108 minutes with 10,000 iterations following 10,000 burn-in iterations, PMCMC takes 11 minutes with the same number of iterations, PVB is the fastest at 1 minute computation time.

We next investigate how the differences observed for the individual parameters will impact the overall inference of interest. A sensitivity plot is presented in Figure 4. This figure is obtained by plotting the average sensitivity against a range of marginal posterior probability thresholds from to . We choose this range because it covers those values most often used in practice.

MSE | LPML | Timing | |||
---|---|---|---|---|---|

(min) | |||||

SVARO | 0.478 | 0.030 | 0.001 | -1842902 | 108 |

PMCMC | 113% | 135% | 509% | -1926620 | 11 |

PVB | 199% | 138% | 510% | NA | 1 |

In terms of the underlying activation threshold, we use two thresholds: the true value of the contrast that corresponds to the top and top of all voxels. Thus, corresponding to a certain activation threshold and a certain probability threshold, the higher the sensitivity, the better the model is in terms of capturing activation. Again, a notable difference is observed when comparing the three methods, with SVARO giving the uniformly highest sensitivity across the entire range of probability thresholds and PVB resulting in the lowest sensitivity. PMCMC is better than PVB but still under-performs relative to SVARO.

We plot the posterior probability maps (PPM) in Figure 5. This figure depicts the locations of the true activations and the posterior probability maps from SVARO. In addition, differences in the probability maps comparing SVARO with PMCMC and PVB are also depicted. Again, SVARO appears to perform the best in producing the highest posterior probabilities for regions that are truly activated. PMCMC is similar to SVARO but its probability on those activated regions are slightly lower than those from SVARO, especially on the boundary. PVB under performs compared with the other two approaches by providing greater posterior probability on null locations while providing smaller posterior probability on actived locations.

### 3.3 Simulation 2

Although the real data examined earlier suggest the existence of heterogeneous AR orders we want to compare the performance of SVARO under a homogeneous AR order assumption, where now the GLM-AR model is correctly specified. To do so, we simulate under the competing GLM-AR model. The AR coefficients are simulated under the LORETA prior and the AR order is set to for every voxel, with prior precision . We set the maximum AR order to when applying SVARO. Thus, PMCMC and PVB are working under the true model while SVARO is working under a more general model.

Table 2 shows the MSE summaries of the estimators. When data are simulated under the competing model, SVARO still exhibits good performance in terms of MSE for the two regression parameters. Its MSEs are only slightly larger the those under the correctly specified model using MCMC posterior simulation. It is worth mentioning that PVB again under performs relative to PMCMC in terms of the hemodynamic response parameter and the AR coefficient.

MSE | LPML | Timing | |||
---|---|---|---|---|---|

(min) | |||||

SVARO | 0.502 | 0.031 | 0.003 | -1817287 | 206 |

PMCMC | 99% | 97% | 54% | -1875900 | 11 |

PVB | 167% | 98% | 49% | NA | 1 |

Figure 6 presents the sensitivity curves. Despite the data being simulated under a constant order AR assumption, SVARO demonstrates similar sensitivity to that of PMCMC. The sensitivity curve for PVB is uniformly smaller than the other two because of the inaccurate estimation of the hemodynamic response parameter.

Overall, these simulation studies show that our SVARO model is more efficient from a statistical point of view than the GLM-AR model when the AR orders of voxels differ—a situation we believe to be more common than not. They also show that our model does not suffer much in terms of statistical efficiency when the AR order is constant across voxels.

## 4 Analysis of Face Perception fMRI Data

We turn our focus back to the face repetition data set that originally motivated our model development and compare results from the two models and three algorithms. In this analysis we use the complete experimental design consisting of famous faces, repeated twice and unfamiliar faces, repeated twice. These four design vectors are then convolved with the canonical HRF as well as its time and dispersion derivatives. An intercept term is also added to the final design matrix for a total of 13 covariates. We allow an AR order up to a maximum of when fitting the SVARO model, and an for the GLM-AR model using both the PMCMC and PVB algorithms. While the choice of an for the latter two approaches may seem arbitrary, this is exactly the justification for the use of the SVARO model where such an arbitrary assumption need not be made. We consider data collected on a single subject in what follows.

Pre-processing steps are applied to the data prior to fitting the Bayesian models. All functional images are aligned to the first image using a six-parameter rigid-body transformation. Then slice-timing correction is performed to set the standard acquisition time as the 12th slice. Images are spatially normalized to a standard EPI image. The global mean, , is computed and each time series is divided by to represent a percentage of . Finally, a high-pass filter with cut-off frequency of Hz is apply to the data and design matrix to remove low frequency signals that arise through scanner drift.

Figure 7 presents the distribution of optimal AR orders estimated from SVARO across voxels. The most frequent order is the zero order, or no autocorrelation in the time series, accounting for approximately 35% of the voxels. Interestingly, the next highest is order 8, with 9.4% of the voxels. Overall, roughly 51% of the voxels exhibit an AR order greater than 3. The existence of these higher orders and the variability in the orders is in general agreement with our exploratory analysis of the face repetition data set (see Section 1.1) and indicates the necessity of our proposed model.

Next, we compare models/algorithms with respect to a particular contrast of interest, namely the effect of “fame”. That is, famous faces versus unfamiliar faces. The estimated posterior mean and standard deviation (SD) maps for the fame contrast, as well as the estimated posterior mean for the order autoregressive coefficient, are displayed in Figure 8. While the posterior SD estimated from SVARO and PMCMC are very close, the estimated posterior means are different between the two approaches. Also, the estimated SD obtained from PVB shows apparent discretization. This is due to a graph-partitioning that is incorporated in the algorithm for the sake of computational speed (Penny et al., 2007). It is clear that the boundaries of these graph-partitioned regions have substantially larger estimated SD than the interior locations. The posterior mean of PVB also seems to exhibit artifacts at the partition boundaries, thought the effect is not as pronounced. The estimated LPML is under our proposed model and is under GLM-AR with MCMC sampling. According to this model selection criterion, our proposed model is preferred.

Finally, we look at the effect of fame using thresholded PPMs. The activation threshold is set to of the global mean value, and the probability threshold is set to . Figure 9 shows the activation regions projected onto the pial surface. We can see that there is a match in terms of a majority of activation regions inferred from SVARO and PMCMC. A closer look reveals that PMCMC tends to make more scattered predictions across the posterior regions of the brain. The number of activation regions from PVB are far greater than the number obtained from the other two approaches, and are more widely dispersed across the brain. From our simulation study results, we suspect that these scattered activated regions are likely false positives.

In terms of computation time, PVB took about hour, PMCMC took day, while SVARO took about week of computation. Much of our MCMC algorithm is amenable to parallel programming, which is an avenue for further development.

SVARO | PMCMC | PVB |
---|---|---|

## 5 Discussion

In this paper, we have developed a Bayesian hierarchical model, SVARO, that allows the AR order to vary spatially across the brain, with the orders themselves displaying a certain level of spatial clustering based on an Ising model. We compared our proposed model with a self-written MCMC sampler for the standard GLM-AR model and the VB implementation for the same model. The results are interesting, under a low SNR ratio, VB seems to suffer from variance overestimation, leading to a much bigger MSE than the other two methods. It is likely that as temporal noise increases, a more vital role is played by the AR correlation that increases the posterior correlation between different parameters and this makes the mean field assumption underlying the VB approximation less accurate.

We have shown that our model outperforms the GLM-AR model not only in terms of accuracy and sensitivity, but also according to formal model selection based on the LPML criterion. Through an application of our proposed model and through an exploratory analysis, we have shown that the typical constant low-order AR assumption can be violated with real fMRI data. It is very likely that this issue, seen in the face repetition data set, is also present in other fMRI data sets.

There is a computational price to be paid for gaining the flexibility we have proposed in our model. Our model takes a longer time to run than either implementation of GLM-AR. This is mostly due to the estimation of the varying AR orders and the associated greater number of parameters to estimate. However, as previously mentioned, there are elements of our MCMC algorithm that are amenable to parallel programming. This will be investigated in future work.

While we have based our model specification on a set of independent Ising processes, one for each possible order of the AR process, another approach would be to assume a Potts model for the orders of the AR coefficients. A Potts model, combined with a Dirichlet process prior for parameters has been investigated for selecting covariates of interest in brain imaging (Johnson et al., 2013). Here we can also apply it to the selection of autoregressive orders to yield a more parsimonious, yet still flexible, model. Investigation of hyper-parameter estimation in the Ising model and the use of alternative spatial models is also of interest, as is increasing the scope of our comparison of methods to include wavelet approaches that focus on long memory errors, or VAR models (Harrison et al., 2003).

## References

- Bezener et al. (2016) Bezener, M., Eberly, L. E., Hughes, J., Jones, G. and Musgrove, D. R. (2016) Bayesian spatiotemporal modeling for detecting neuronal activation via functional magnetic resonance imaging. Handbook of Big Data Analytics, Springer Handbooks of Computational Statistics. Springer.
- Bullmore et al. (1996) Bullmore, E., Brammer, M., Williams, S. C., Rabe-Hesketh, S., Janot, N., David, A., Mellers, J., Howard, R. and Sham, P. (1996) Statistical methods of estimation and inference for functional mr image analysis. Magnetic Resonance in Medicine, 35, 261–277.
- Bullmore et al. (2004) Bullmore, E., Fadili, J., Maxim, V., Şendur, L., Whitcher, B., Suckling, J., Brammer, M. and Breakspear, M. (2004) Wavelets and functional magnetic resonance imaging of the human brain. Neuroimage, 23, S234–S249.
- Castruccio et al. (2016) Castruccio, S., Ombao, H. and Genton, M. G. (2016) A multi-resolution spatio-temporal model for brain activation and connectivity in fmri data. arXiv preprint arXiv:1602.02435.
- Chang and Glover (2010) Chang, C. and Glover, G. H. (2010) Time–frequency dynamics of resting-state brain connectivity measured with fmri. Neuroimage, 50, 81–98.
- Della-Maggiore et al. (2002) Della-Maggiore, V., Chau, W., Peres-Neto, P. R. and McIntosh, A. R. (2002) An empirical comparison of spm preprocessing parameters to the analysis of fmri data. Neuroimage, 17, 19–28.
- Fadili and Bullmore (2002) Fadili, M. and Bullmore, E. (2002) Wavelet-generalized least squares: a new blu estimator of linear regression models with 1/f errors. NeuroImage, 15, 217–232.
- Friston et al. (2000) Friston, K., Josephs, O., Zarahn, E., Holmes, A., Rouquette, S. and Poline, J.-B. (2000) To smooth or not to smooth?: Bias and efficiency in fmri time-series analysis. NeuroImage, 12, 196–208.
- Friston et al. (1995) Friston, K. J., Holmes, A. P., Poline, J., Grasby, P., Williams, S., Frackowiak, R. S. and Turner, R. (1995) Analysis of fmri time-series revisited. Neuroimage, 2, 45–53.
- Gelfand and Dey (1994) Gelfand, A. E. and Dey, D. K. (1994) Bayesian model choice: asymptotics and exact calculations. Journal of the Royal Statistical Society. Series B (Methodological), 501–514.
- Gelman and Meng (1998) Gelman, A. and Meng, X.-L. (1998) Simulating normalizing constants: From importance sampling to bridge sampling to path sampling. Statistical science, 163–185.
- George and McCulloch (1993) George, E. I. and McCulloch, R. E. (1993) Variable selection via gibbs sampling. Journal of the American Statistical Association, 88, 881–889.
- Gössl et al. (2001) Gössl, C., Auer, D. P. and Fahrmeir, L. (2001) Bayesian spatiotemporal inference in functional magnetic resonance imaging. Biometrics, 57, 554–562.
- Harrison et al. (2003) Harrison, L., Penny, W. D. and Friston, K. (2003) Multivariate autoregressive modeling of fmri time series. Neuroimage, 19, 1477–1491.
- Higdon (1998) Higdon, D. M. (1998) Auxiliary variable methods for markov chain monte carlo with applications. Journal of the American Statistical Association, 93, 585–595.
- Ising (1925) Ising, E. (1925) Beitrag zur theorie des ferromagnetismus. Zeitschrift für Physik, 31, 253–258.
- Jeong et al. (2013) Jeong, J., Vannucci, M. and Ko, K. (2013) A wavelet-based bayesian approach to regression models with long memory errors and its application to fmri data. Biometrics, 69, 184–196.
- Johnson et al. (2013) Johnson, T. D., Liu, Z., Bartsch, A. J. and Nichols, T. E. (2013) A bayesian non-parametric potts model with application to pre-surgical fmri data. Statistical methods in medical research, 22, 364–381.
- Johnson and Rossell (2012) Johnson, V. E. and Rossell, D. (2012) Bayesian model selection in high-dimensional settings. Journal of the American Statistical Association, 107, 649–660.
- Kim et al. (2010) Kim, S., Smyth, P. and Stern, H. (2010) A bayesian mixture approach to modeling spatial activation patterns in multisite fmri data. IEEE transactions on medical imaging, 29, 1260–1274.
- Lee et al. (2014) Lee, K.-J., Jones, G. L., Caffo, B. S. and Bassett, S. S. (2014) Spatial bayesian variable selection models on functional magnetic resonance imaging time-series data. Bayesian Analysis (Online), 9, 699.
- Li et al. (2015) Li, F., Zhang, T., Wang, Q., Gonzalez, M. Z., Maresh, E. L., Coan, J. A. et al. (2015) Spatial bayesian variable selection and grouping for high-dimensional scalar-on-image regression. The Annals of Applied Statistics, 9, 687–713.
- Locascio et al. (1997) Locascio, J. J., Jennings, P. J., Moore, C. I. and Corkin, S. (1997) Time series analysis in the time domain and resampling methods for studies of functional magnetic resonance brain imaging. Human brain mapping, 5, 168–193.
- Makni et al. (2006) Makni, S., Ciuciu, P., Idier, J. and Poline, J.-B. (2006) Joint detection-estimation of brain activity in fmri using an autoregressive noise model. In Biomedical Imaging: Nano to Macro, 2006. 3rd IEEE International Symposium on, 1048–1051. IEEE.
- Meyer (2003) Meyer, F. G. (2003) Wavelet-based estimation of a semiparametric generalized linear model of fmri time-series. IEEE transactions on medical imaging, 22, 315–322.
- Mitchell and Beauchamp (1988) Mitchell, T. J. and Beauchamp, J. J. (1988) Bayesian variable selection in linear regression. Journal of the American Statistical Association, 83, 1023–1032.
- Musgrove et al. (2016) Musgrove, D. R., Hughes, J. and Eberly, L. E. (2016) Fast, fully bayesian spatiotemporal inference for fmri data. Biostatistics, 17, 291–303.
- O’Hara et al. (2009) O’Hara, R. B., Sillanpää, M. J. et al. (2009) A review of bayesian variable selection methods: what, how and which. Bayesian analysis, 4, 85–117.
- Penny et al. (2007) Penny, W., Flandin, G. and Trujillo-Barreto, N. (2007) Bayesian comparison of spatially regularised general linear models. Human brain mapping, 28, 275–293.
- Penny et al. (2003) Penny, W., Kiebel, S. and Friston, K. (2003) Variational Bayesian inference for fMRI time series. NeuroImage, 19, 727–741.
- Penny et al. (2005) Penny, W. D., Trujillo-Barreto, N. J. and Friston, K. J. (2005) Bayesian fMRI time series analysis with spatial priors. NeuroImage, 24, 350–362.
- Schmidt and Makalic (2013) Schmidt, D. F. and Makalic, E. (2013) Estimation of stationary autoregressive models with the bayesian lasso. Journal of Time Series Analysis, 34, 517–531.
- Shu et al. (2015) Shu, H., Nan, B. and Koeppe, R. (2015) Multiple testing for neuroimaging via hidden markov random field. Biometrics, 71, 741–750.
- Skudlarski et al. (1999) Skudlarski, P., Constable, R. T. and Gore, J. C. (1999) Roc analysis of statistical methods used in functional mri: individual subjects. Neuroimage, 9, 311–329.
- Smith and Fahrmeir (2007) Smith, M. and Fahrmeir, L. (2007) Spatial bayesian variable selection with application to functional magnetic resonance imaging. Journal of the American Statistical Association, 102, 417–431.
- Stanley et al. (1987) Stanley, H. E., Stauffer, D., Kertesz, J. and Herrmann, H. J. (1987) Dynamics of spreading phenomena in two-dimensional ising models. Physical review letters, 59, 2326.
- Swendsen and Wang (1987) Swendsen, R. H. and Wang, J.-S. (1987) Nonuniversal critical dynamics in monte carlo simulations. Physical review letters, 58, 86.
- Teng et al. (2016) Teng, M., Johnson, T. and Nathoo, F. (2016) A comparison of variational bayes and hamiltonian monte carlo for bayesian fmri time series analysis with spatial priors. arXiv preprint arXiv:1609.02123.
- Wang et al. (2007) Wang, H., Li, G. and Tsai, C.-L. (2007) Regression coefficient and autoregressive order shrinkage and selection via the lasso. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 69, 63–78.
- Woolrich et al. (2004a) Woolrich, M. W., Behrens, T. E. and Smith, S. M. (2004a) Constrained linear basis sets for hrf modelling using variational bayes. NeuroImage, 21, 1748–1761.
- Woolrich et al. (2004b) Woolrich, M. W., Jenkinson, M., Brady, J. M. and Smith, S. M. (2004b) Fully bayesian spatio-temporal modeling of fmri data. IEEE transactions on medical imaging, 23, 213–231.
- Woolrich et al. (2001) Woolrich, M. W., Ripley, B. D., Brady, M. and Smith, S. M. (2001) Temporal autocorrelation in univariate linear modeling of fmri data. Neuroimage, 14, 1370–1386.
- Worsley and Friston (1995) Worsley, K. J. and Friston, K. J. (1995) Analysis of fmri time-series revisitedâagain. Neuroimage, 2, 173–181.
- Zarahn et al. (1997) Zarahn, E., Aguirre, G. K. and D’Esposito, M. (1997) Empirical analyses of bold fmri statistics. NeuroImage, 5, 179–197.