Correcting for Multiple Comparison in fMRI

I like the statements above (source: FSL training slides). In doing statistical inference or hypothesis testing of the fMRI dataset, we face multiple comparison issues. In the earlier post, we say we adopt the mass univariate approach or GLM but we don't know which voxels in the brain the effect is significant. Under the null hypothesis H₀, we say that all voxel activations are merely due to chance. Each voxel has a random value. For the hypothesis testing, e.g. we are repeating the same statistical test with different voxels until the whole brain is covered. If the number of voxels = 100,000 with α = 0.05, then chances are we may encounter 5,000 false positives, i.e. voxels wrongly concluded as being active or related to our manipulation. This basically means that we may get a false positive (Type-I error) in our inference if we don't do a correction properly.

There are some familiwise error (FWE) correction methods. The most stringent and conservative correction of multiple comparisons is the Bonferroni method. There are other variants such as Sidak-Bonferroni, Holm, and Tukey HSD methods. In general, Bonferroni correction is useful if and only if the tests done are statistically independent, that is, it assumes adjacent points are not correlated. This assumption is violated in fMRI data because a strong activation in one voxel is always accompanied by activity in neighboring voxels. Why? In the real brain, the neuronal activation of a point in space is always correlated with the surrounding activation. 

Nichols & Hayasaka (2003) wrote a nice and a bit technical summary of popular methods used to correct for multiple comparisons in the fMRI dataset.

Random Field Theory (RFT)
A nice summary of RFT pipeline is described by Ashby's "Statistical Analysis of fMRI Data" (2011). RFT is applied here to make a statistical inference or hypothesis testing.

(1). What is a Multivariate Normal Distribution?
A univariate distribution is a probability distribution of one random variable. For example: the distribution of height. When the probability follows a normal or Gaussian distribution, it is called univariate normal distribution. If this is applied to a higher dimension containing n-random variables, it becomes a multivariate problem, e.g. a set of variables defining a person such as height, age, arm spread, etc. A multivariate normal distribution is very useful in fMRI statistics because it is used to model the error or residual portion in GLM and PCA/ICA. In a multivariate normal distribution, the relationship among its random variables is a pairwise linear relationship.

(2). What is a Gaussian random field? 
A Gaussian random field (GRF) is an ordered collection of random variables. In GRF, the random variables in the field have a joint-multivariate normal distribution. From our fMRI stats stage using the GLM, we obtain a set of statistical parametric maps, each voxel has its z-statistic. Hence, we are no longer working with the voxelwise image value obtained from the MRI machine, but on the z-scores obtained from the GLM. Under the null hypothesis, these voxels are then treated as a random field. 

Two assumptions are noted while treating the brain voxels as GRF. The first assumption is that the spatial smoothness of the fMRI signal is constant over the brain, and the second assumption is that the spatial autocorrelation function has a specific shape.

(3). What is the basic intuition of RFT?
Suppose we have a 2D slice of z-statistic with 100×100 voxels under the null hypothesis. To fulfill the assumptions, we smooth the data using the Gaussian kernel, a type of spatial filter performed in the usual fMRI preprocessing steps. As a result, instead of the z-map, what we have is resels of 10 x 10. Worsley said that resels are like the smoothed version of the original voxels. They represent the number of independent datasets. What comes next is to determine a cutoff or threshold that will be applied to the resels. Suppose our 2D image which has 100 resels under the null hypothesis is shown below, the Y-axis is the z-score. 

(source: UCL intro slides to RFT).

Let T be a certain threshold applied to the field, such that there are regions that will be above T called the excursion set. The corresponding thresholded image is the Euler Characteristic (χ_T). Intuitively, Euler Characteristic is roughly equal to the # peaks (blobs) with respect to the # valleys (holes, in between peaks) in the image after the thresholding. Importantly, for sufficiently high threshold T we can only see only a few peaks that exist at voxels, where their largest z-value > T. 

Under the null hypothesis, these peaks correspond to # false positives. Hence, the probability of at least observing one false positive is equal to the probability of seeing at least one peak, i.e. χ_T ≥ 1. It turns out that the α_E  is roughly equivalent to the expected value of χ_T , the formula of which has been derived by (Worsley, et al. 1992, 1996). Subsequently, the corrected p-value can be computed through the probability of χ_T ≥ 1 and this is proportional to the volume size and inversely proportional to the smoothness criterion.

(4). Cluster-forming threshold
As with the usual hypothesis testing, we assume that the null hypothesis has to be true for an image voxel to apply GRF to it. Here, we start by forming a so-called cluster step-by-step until we obtain a supervoxel in 3D. We can iterate the process to obtain the distribution of χ_T and z-values. What determines this relationship is the resels only, not the original z-statistic values. If T is the threshold that corresponds to E(χ_T) of 0.05, then by using T, we can expect that any remaining peaks have a probability of 0.05 that they have occurred by chance given that the null hypothesis is true. This α_E value is akin to the Type-I error in the usual hypothesis testing. This Z- threshold value is called the cluster-forming threshold. 

Example: for our 2D resels of 100 resels, the equation gives E(χ_T) = 0.049 for a threshold Z = 3.80, i.e. the probability of getting one or more blobs where z > 3.80 is 0.049. What does it mean? If we threshold our image at Z = 3.80, we can conclude that any blobs that remain have a probability of less than or equal to 0.05 that they have occurred by chance. In other words, this is the case when we have to reject the null. See the illustration above.

For more technical details, this is a good online summary. In SPM and FSL, once we set the cluster-forming threshold (Z = 2.30, with p < 0.05), GRF is automatically applied to obtain the contiguous cluster that survives the threshold while at the same time controlling the Type-I error at 0.05.

(5). How conservative is GRF?
It is good to note that # resel counts is not the same as the number of independent tests. GRF is conservative when the # dof is small, but less conservative than the usual FWE corrections. Without smoothing, the GRF assumes that voxels are independent and the threshold T is equal to the Sidak-Bonferroni threshold, which is a more stringent way of correction. As the kernel size increases, smoothing and spatial correlation increase, resels size increases, and so T decreases. The choice of the most optimal smoothness places a serious problem in applying GRF to the fMRI dataset. Another thing to note that in a 3D brain volume, the resel takes into account the shape of the volume. A more sophisticated method using search volume can be used (Worsley, et al. 1996).

False Discovery Rate
FDR approach is totally different from controlling familywise error. Rather than limiting the errors V themselves to occur, we limit the proportion of significant results that are false positive, that is, V/R. Thus, FDR is the probability of V/R and the corrected p-value guarantees that FDR < q. The goal of FDR is finally to find regions with p-values associated with q smaller than 0.05, say.

What is the algorithm in finding q?
1. Convert N number of z-values to p-values. Let the value of q = 0.05.
2. In the case of activation, p- is the area under the null distribution to the right of z
3. In the case of no activation, p- is the area under the null distribution to the left of z.
4. Rank or order all p-values from the smallest to the largest. Let the k-th smallest p-value as P[k].
5. The z-value of the k-th smallest p-value is significant if P[k] < qk/N;  otherwise it isn't significant.

Compared to Bonferroni and GRF thresholding, FDR almost always indicates more and larger significant regions in the brain. But this is also because FDR does not control for Type-I error α_E, but rather, the rate of its occurrence as a whole.

Permutation-based Method
The permutation-based method is considered non-parametric as it doesn't depend on the nature of the distribution and correlation that exists in our image data. Recent opinions support the permutation-based approach as an alternative to GLM-based statistics. To see why this makes sense, the task-based block design will be used. A participant performs a task in one block and rest in the other. If a specific location of the brain is required for task performance, we should expect a task-related activation during task block as reflected in the BOLD response. However, this activation does not appear during the rest block.

Now, we reshuffle or scramble the data across time points or TRs. If a particular region is involved in the task, we should not expect any strong correlation between the predicted BOLD response and the task performance because the relationship has collapsed. The "reshuffled data" behaves just like a random noise. Reshuffling does not change the mean and variance of the dataset, that is, we don't tamper the data. Although attractive, the method has one drawback: the computation time is much longer as you have to keep randomizing the time points. Although reshuffling the data for 500 or 1000 times is quite a norm, there is no consensus of how many times this should be repeated.

How can this method correct for multiple comparison problem? It sounds simple, but I think I have to read more. Therefore, a separate blog entry will attempt to discuss this method with my own stuff to play. Anyways... what it does in correcting for multiple comparison is mentioned below:
1. Randomly shuffle the data point (TRs). Reorder every voxel time series in the ROI so that the
    TRs agree with the shuffle.
2. Run the GLM to produce a statistical parametric map. Find the maximum z-value (z_max) of the map.
3. Repeat Step 1 & 2 so that eventually we have M maps and thus z_max. Usually M = 500 or 1000.
4. Rank the order of maximum z-values, the largest being z_max[1] down to z_max[M].
5. The corrected threshold T = z_max.(Mα_E + 1), where α_E is the familywise error.