A General Statistical Analysis for fMRI Data

K. J. Worsley¹², C. Liao¹, M. Grabove¹, V. Petre², B. Ha², A.C. Evans²
¹Department of Mathematics and Statistics, McGill University
²McConnell Brain Imaging Centre, Montreal Neurological Institute, McGill University

Abstract
Many methods are available for the statistical analysis of fMRI data that range from a simple linear model for the response and a global autoregressive model for the temporal errors (SPM), to a more sophisticated non-linear model for the response with a local state space model for the temporal errors (Purdon, et al., 1998). We have written two Matlab programs, FMRISTAT and MULTISTAT (http://www.bic.mni.mcgill.ca/users/keith) , that seek a compromise between validity, generality, simplicity and execution speed.

Methods

FMRISTAT performs a statistical analysis of a single run of an fMRI data set. The statistical method is based on a linear model with correlated errors. The design matrix of the linear model is specified either directly or via event times, durations and weights. This is first convolved with a hemodynamic response function, chosen by default to be a gamma function with a mean lag of 6 seconds and a standard deviation of 3 seconds, timed to coincide with the acquisition of each slice (Lange and Zeger, 1997). Drift is removed by adding polynomial covariates in the frame times, up to degree 3 by default, to the design matrix. It is important to note that these drift terms are not convolved with the hemodynamic response function. The correlation structure is modeled as an autoregressive process of degree 1 (Bullmore et al., 1996). At each voxel, the autocorrelation parameter was estimated from the least squares residuals using the Yule-Walker equations, after a bias correction for correlation between the residuals induced by the linear model. Without this correction, the autocorrelation parameter is negatively biased by about -0.05. The autocorrelation parameter is then regularized by spatial smoothing with a 15mm FWHM (default) Gaussian filter, then used to `whiten' the data and the design matrix. The linear model is then re-estimated by using least squares in a second run through the whitened data, to produce estimates of effects and their standard errors. This second run requires finding the pseudoinverse of a different whitened design matrix at each voxel. To reduce computations, this was evaluated only for the autocorrelation parameter rounded to the nearest second decimal, then cached.

MULTISTAT combines the output from FMRISTAT across runs within sessions, sessions within subjects, then finally subjects within a population in a hierarchical random effects analysis. Runs are combined with another linear model for the run effects (as data), weighted inversely by the square of their standard errors. A regularized random effects analysis was performed by first estimating the ratio of the random effects variance (VAR_random) to the fixed effects variance (VAR_fixed), then regularizing this ratio by spatial smoothing with a Gaussian filter with a FWHM of FWHM_varatio=15mm (default). The variance of the effect was then estimated by the smoothed variance ratio multiplied by the fixed effects variance to achieve higher degrees of freedom at the expensive of introducing a small bias:

VAR_effect=VAR_fixed*smooth(VAR_random/VAR_fixed).

The degrees of freedom of the smoothed variance ratio was estimated from random field theory by

DF_ratio=DF_random*(2*(FWHM_varatio/FWHM_data)^2+1)^(3/2),

where DF_random is the degrees of freedom of the second linear model for the runs, and FWHM_data is the FWHM of the original data (6mm by default). The final degrees of freedom of the effect is estimated by

DF_effect=1/(1/DF_ratio+1/DF_fixed),

where DF_fixed is the sum of the degrees of freedom of the fixed effect analysis, equal to the sum of the degrees of freedom of the preceding separate FMRISTAT analyses. Thus the FWHM_varatio parameter acts as a convenient way of providing an analysis mid-way between a random effects and a fixed effects analysis; setting FWHM_varatio=0 (no smoothing) produces a random effects analysis; setting FWHM_varatio to infinity, which smoothes the variance ratio to one everywhere, produces a fixed effects analysis. In practice, we choose FWHM_varatio to produce a final DF_effect which is at least 100, so that errors in its estimation do not greatly affect the distribution of test statistics.

The resulting T statistic images are thresholded using the minimum given by a Bonferroni correction and non-isotropic random field theory (Worsley et al, 1996, 1999).

Results

When applied to fMRI data with TR=3s, we have found an autocorrelation parameter of 0.2-0.5 in grey matter, and ~0 elsewhere. We found negligible random effects between runs in the same session, and between sessions on the same subject, but substantial random effects (up to 5 x increase in standard deviation) between subjects within a population.

References

Bullmore, E.T. et al. (1996). Magnetic Resonance in Medicine, 35:261-277.
Lange, N. and Zeger, S.L. (1997). Applied Statistics, 14:1-29.
Purdon, P.L. et al. (1998). NeuroImage, 7:S618.
Worsley, K.J. et al. (1996). Human Brain Mapping, 4:58-73.
Worsley, K.J. et al. (1999). Human Brain Mapping, 8:98-101.