A comparison of random field theory and permutation methods for the statistical analysis of MEG data

Abstract

We describe the use of random field and permutation methods to detect activation in cortically constrained maps of current density computed from MEG data. The methods are applicable to any inverse imaging method that maps event-related MEG to a coregistered cortical surface. These approaches also extend directly to images computed from event-related EEG data. We determine statistical thresholds that control the familywise error rate (FWER) across space or across both space and time. Both random field and permutation methods use the distribution of the maximum statistic under the null hypothesis to find FWER thresholds. The former methods make assumptions on the distribution and smoothness of the data and use approximate analytical solutions, the latter resample the data and rely on empirical distributions. Both methods account for spatial and temporal correlation in the cortical maps. Unlike previous nonparametric work in neuroimaging, we address the problem of nonuniform specificity that can arise without a Gaussianity assumption. We compare and evaluate the methods on simulated data and experimental data from a somatosensory-evoked response study. We find that the random field methods are conservative with or without smoothing, though with a 5 vertex FWHM smoothness, they are close to exact. Our permutation methods demonstrated exact specificity in simulation studies. In real data, the permutation method was not as sensitive as the RF method, although this could be due to violations of the random field theory assumptions.

Introduction

Magnetoencephalography (MEG) is used to image electrical activity in the brain. Clusters of thousands of synchronously activated pyramidal cortical neurons are believed to be the main generators of MEG signals. In particular, the currents associated with their large dendritic trunks, which are locally oriented in parallel and perpendicular to the cortical surface, are the primary source of the neuromagnetic fields outside the head (Hämäläinen et al., 1993). Imaging approaches to the MEG inverse problem exploit this concept by restricting the reconstruction to elemental sources (dipoles) oriented normally to the cortical surface (Dale and Serano, 1993). Consequently, a commonly used approach extracts a tessellated representation of the cerebral cortex from a coregistered MR image and solves the inverse problem for a current dipole located at each vertex of the tessellated surface. Since the position and orientation of the dipoles is fixed, image reconstruction is a linear problem and can be solved using standard techniques (Baillet et al., 2001, Hämäläinen et al., 1993, Katila and Karp, 1983, Phillips et al., 1997). However, the highly convoluted nature of the human cortex requires the use of many thousands of dipoles for an accurate representation of the cortical surface. The inverse problem becomes hugely underdetermined and the resulting current density maps (CDMs) are of low resolution; interpretation is further confounded by the presence of additive noise exhibiting a highly nonuniform spatial correlation.

As with fMRI images, objective assessment of CDMs requires a principled approach to identifying regions of activation. The analysis of CDMs involves testing thousands of hypothesis (one per surface element) for statistically significant experimental effects. This raises the possibility of large numbers of false-positives simply as a result of multiple hypothesis testing. To effectively control the number of false-positives over all tests, we must therefore consider the multiple hypothesis-testing problem. Many false-positive measures have been proposed in this context, including familywise error rate, expected false discovery rate, per-comparison error rate and per-family error rate (Nichols and Hayasaka, 2003). The standard approach, and the one investigated in this paper, is to control the Familywise Error Rate (FWER), i.e., the chance of one or more false-positives under the null hypothesis.

The simplest approach to controlling the FWER is the Bonferroni correction method (Hochberg and Tamhane, 1987, Nichols and Hayasaka, 2003). This method produces conservative thresholds unless the tests are independent, a case that is rarely true in neuroimaging experiments and certainly not for the smooth images reconstructed from MEG data. Other methods that consider the spatial dependence of the data make inferences based on the global maximum distribution. The FWER is directly related to the maximum statistic; one or more voxels T_i will exceed the threshold u under the null hypothesis H₀ only if the maximum exceeds the threshold:

$$P\left(\text{FWER}\right)=P({\cup }_i{T}_i>u|H_o)=P({\max }_i{T}_i>u\left|H_o\right)=1-F_{\max T|H_o}\left(u\right)=1-(1-\alpha )=\alpha$$

Therefore, we can control the FWER at level α, if we choose the threshold u to be in the (1 − α) 100th percentile of the maximum distribution.

Random Field (RF) theory methods approximate the upper tail of the maximum distribution F_max using the expected value of the Euler characteristic of the thresholded image (Worsley et al., 1996). They are implemented in various software packages (SPM-http://www.fil.ion.ucl.ac.uk, VoxBo-http://www.voxbo.org, and FSL-http://www.fmrib.ox.ac.uk/fsl among others), and are typically used in PET and fMRI studies. However, RF theory relies on several assumptions including the following: the image has the same parametric distribution at each spatial location, the point spread function has two derivatives at the origin, sufficient smoothness to justify application of continuous RF theory, and a sufficiently high threshold for the asymptotic results to be accurate.

Resampling methods are a different approach to controlling the FWER that exploit the information contained in the data to estimate the empirical distribution of the maximum statistic. They do not assume parametric distributions, they adapt to underlying correlation patterns in the data, and are now computationally feasible. The two main categories are bootstrap-based, which allow for a general modeling framework, and permutation-based, which require some knowledge of exchangeability conditions under the null hypothesis. Here, we only consider permutation methods since they are exact, that is, they give precise control of the FWER, while bootstrap methods are only asymptotically exact. Furthermore, the permutation approach relies on a less restrictive exchangeability condition than the requirement in the bootstrap that samples are independent and identically distributed.

Permutation and RF theory methods have been applied widely in functional (Andrade et al., 2000, Nichols and Holmes, 2001, Worsley et al., 1992, Worsley et al., 1996) and structural (Bullmore et al., 1999, Chung, 2001, Pantazis et al., 2004, Sowell et al., 1999a, Sowell et al., 1999b, Thompson et al., 2001, Thompson et al., 2003) brain imaging. However, until recently, error rate control in MEG experiments has drawn little attention. Dale et al. (2000) normalized the CDMs using an estimate of the background noise variance at each cortical element. These normalized images follow a t distribution under the null hypothesis of Gaussian background noise. Thresholding of the resulting statistical maps was then used to detect significant activation, however, the multiple comparisons problem was not addressed. Barnes and Hillbrand (2003) presented an application of RF theory to MEG data but their method is specifically tailored to beamforming solutions rather that the general linear inverse methods. Carbonell et al. (2004) used Hotelling's T² random fields to localize significant MEG/EEG activation in time, and then t statistics to achieve spatial localization. Permutation tests were applied by Blair and Karnisky (1994) for the analysis of EEG data as recorded on an array of electrodes, and by Pantazis et al. (2003) for the analysis of MEG data in reconstructed cortical maps of brain activation. An alternative permutation scheme, proposed by Singh et al. (2003), detects event-related synchronization or desynchronization components in an MEG study involving visual stimulation and a Linearly Constraint Minimum Variance (LCMV) beamformer applied to data decomposed into multiple frequency bands. The current work presents a novel general RF theory-based method to control FWER in MEG. Also, it extends the results in Pantazis et al. (2003) to extract thresholds for each time-point. Finally, we compare RF theory and permutation methods in terms of specificity, sensitivity, and possible limitations.