Use of the Gibbs sampler for error covariance estimation in Statistical Parametric Mapping

Wellcome Dept. of Imaging Neuroscience, Institute of Neurology, WC1N 3BG, London, UK

Subject: Modeling & Analysis

Abstract
Introduction:
In Statistical Parametric Mapping (SPM), the estimation of non-sphericity or the error covariance matrix is a crucial step in the formation of statistical maps. An important example of non-sphericity is serial correlations in fMRI. In the current development version of SPM, these serial correlations are modelled using an autoregressive moving average (ARMA(1, 1)) model. The hyperparameters of this model are estimated using linear expansions of the covariance matrix and the Restricted Maximum Likelihood (ReML) method (1). The hyperparameters are estimated over all cerebral voxels and used to form a serial correlation matrix which enters the next step as a known variable. A suitable test-statistic is then formed using a voxel-wise variance estimate and a Satterthwaite approximation (2). We refer to this process as ReML-S. An alternative approach is to compute the error covariance matrix from the ARMA or AR coefficients estimated by Gibbs sampling (3). Like ReML-S, the Gibbs sampler provides unbiased estimates but eschews the linearization of the covariance components. Importantly, the output of the Gibbs sampler provides the marginal posterior distribution of the contrast weighted parameters, (e.g. activations). This makes two things possible. First, it allows one to make inferences at single voxels while parameterizing the error covariance matrix with multiple hyperparameters. Secondly, it enables Bayesian inference, i.e. posterior probability maps (PPMs, 1).

Methods:
We assume a (generalized) linear model based on an AR(p) process. This is Y X , _k a_k (Eq. 1), where Y is the data, X the design matrix, the parameters of the linear model, _k are lagged errors, a_k are the p AR coefficients, and the innovation N(0, ² I). We use non-informative priors. The Gibbs sampler is based on drawing iteratively from the full conditional distributions for the parameters (, a_k, ²). From these we obtain marginal posterior distributions of all .

Results:
As an example, we consider simulated data from two subjects at a single voxel assuming an AR(1) process. We use two regressors for each subject and are interested in a difference between subjects. Typically, the probability density function (PDF) of the test statistic as approximated by ReML-S (without pooling over voxels) is too narrow and ReML-S slightly overestimates p-values in the (interesting) upper tail of the distribution. This overestimation results in an increased false positive rate, i.e. an invalid test. This is confirmed by observing (over 1000 realizations) that at a 0.05 test level, the false positive rate is 0.066 for ReML-S and 0.048 for the Gibbs sampler.

Conclusion:
We have implemented the Gibbs sampler for both AR and ARMA models. The approach can be either used to improve classical inference or to generate PPMs. The Gibbs sampler can be applied to any neuroimaging time series, in particular to fMRI and EEG/MEG data.

References:
(1) K.J. Friston et al., Neuroimage, submitted.
(2) K.J. Worsley and K.J. Friston, Neuroimage (2): 173-181.
(3) J.B. Carlin et al., Bayesian Data Analysis, 1995.