The fMRI data was first converted to percentage of whole volume. The statistical analysis of the percentages was based on a linear model with correlated errors. The design matrix of the linear model was first convolved with a hemodynamic response function modelled as a difference of two gamma functions timed to coincide with the acquisition of each slice. Temporal drift was removed by adding a cubic spline in the frame times to the design matrix (one covariate per 2 minutes of scan time), and spatial drift was removed by adding a covariate in the whole volume average. The correlation structure was modelled as an autoregressive process of degree 1. At each voxel, the autocorrelation parameter was estimated from the least squares residuals using the Yule-Walker equations, after a bias correction for correlations induced by the linear model. The autocorrelation parameter was first regularized by spatial smoothing, then used to `whiten' the data and the design matrix. The linear model was then re-estimated using least squares on the whitened data to produce estimates of effects and their standard errors.

In a second step, runs, sessions and subjects were combined using a mixed effects linear model for the effects (as data) with fixed effects standard deviations taken from the previous analysis. This was fitted using ReML implemented by the EM algorithm. A random effects analysis was performed by first estimating the the ratio of the random effects variance to the fixed effects variance, then regularizing this ratio by spatial smoothing with a Gaussian filter. The variance of the effect was then estimated by the smoothed ratio multiplied by the fixed effects variance. The amount of smoothing was chosen to achieve 100 effective degrees of freedom.