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.

Group Analysis in Functional MRI - a brief overview

Previously, I mentioned subject-level or first-level GLM analysis that is performed on individual subject data. This step gives us within-subject contrasts across parameter estimates, e.g. estimates of brain activity due to visual stimulation vs. rest. In this blog entry, a second-level analysis will be discussed which provides group-level inference on whether they are significantly different from zero, e.g., group status, behavioral performance, etc. Statistically speaking, the goal of the standard group analyses is to do inference of a set of acquired sample to a larger population. How can we perform group analysis in neuroimaging data? There are two ways we can model the "effects" of a behavioral study on the brain, that is, fixed effects and random effects. A statistic that takes both variations into account is called mixed-effects modeling in psychology or behavioral sciences. 

Consider the following scenario. You want to know the distribution of hair length between adult men and women, but doing this for all world population is just impossible. You start to select four men and women. For every single person, there is bound to be differences in hair length which represents a within-subject variability. In addition, different people in our study will have different hair cut or style, a form of between-subject variability. You can sample a strain of hair of each person or multiple hairs.

1. Fixed effects model: we assume that any systematic variability introduced in our measurement is due to within-subject variability only. In other words, the four men (or women) are more or less equal in terms of hairstyle.
2. Random effects model: this model assumes a more realistic scenario where a man (or woman) is randomly sampled from a population. The model takes into account another factor called "subject", and any unpredictable variation in the data is contributed by an individual or between-subject differences. 

In functional MRI, we employ mixed-effects modeling as accurate statistics when analyzing group data. This is done by taking sample voxel's time series per subject. Theoretically speaking, a fixed-effects model is used when a single subject has multiple runs. In order to assess multiple subjects that have multiple runs, we have to use mixed-effects modeling because every person is understood to be unique and we want to make a statement that can generalize to all people.

Fixed Effects Analysis
Suppose Y_i = X_ib_i + ε_i  is the fixed effects model of our neuroimaging data of N-subjects, each is acquired with n time points. The subject-level GLM consists of p parameter estimates. The within-subject residual ε_i has a multivariate normal distribution with a mean vector 0 and variance-covariance matrix of Σ_i = σ_w².I. Under the null hypothesis, H₀: c^Tb_G = 0, we can then compute:

which obeys a t-distribution with # dof = N(n – p). Note how we incorporate the group-level PE from the average of individual subject's PEs, and within-subject error variance from the pooled individual error variance.

Random Effects Analysis
Our group modeling will now have both fixed and random effects. Suppose Y_i = X_ib_i + ε_i  is the model of our neuroimaging data of N-subjects, each is acquired with n time points. The within-subject parameter estimates and residual are b_i and ε_i respectively. In the random-effects model, we further treat the individual subject as coming from a set of the population with its own between-subject variability. Thus, b_i can be further decomposed into b_G + εGi. the group parameter estimates and group error of which the subject belongs. Under the null hypothesis, H₀: c^Tb_G = 0, we can then compute:

which obeys a t-distribution with # dof = N – 1. Note that the numerator is the same as the fixed effects model, but the denominator contains more terms. The denominator in random effects is comparatively larger than that in fixed-effects model. So, if you use fixed effects for group analysis, you tend to get unnecessary activation, that is, a false positive.

Notes on Linear Mixed Model
The statistics mentioned in this post are also known by linear mixed model (LMM). With mixed modeling, the model will fit the average intercept and slope as a fixed effect, and each subject is assigned a different intercept/slope. Hence, LMM is an extension of GLM, but this is unlike a standard linear regression model which has only fixed effects. LMM approach is powerful to model data where there is a hierarchical or nested structure, or we deal with multilevel modeling. LMM can also be used in the case of repeated-measures design, or cases when there are some missing data.

Let's take another classic example of LMM: modeling distribution of performance of college students in New York City. Variation in the performance can primarily be due to a random variation at the individual (each student) level. But suppose the city has many colleges. Each college can also be seen as a contributing factor to the performance of each of the individuals at that school, say, due to the good teachers and attitude of the school. Hence those observations cannot be treated as fully independent of each other, but dependent on a higher level grouping called (which college?) - breaking the major assumptions of more traditional linear models. This example reflects different sources of randomness which are in a hierarchy, i.e. individuals-classes-college. This resembles our fMRI dataset, as we move up from an individual subject data to a group-based inference.

LMM is also common for two explanatory variables that are not independent. For example: suppose we want to predict the pitch level as a function of two independent variables: age and gender. Following GLM notation, we can write down:

pitch ~ nationality + gender + error

Notice how age and gender are dependent as coming from the same person or subject. Taking into consideration fixed and random effects with LMM, we rewrite the equation to be:

pitch ~ nationality + gender + (1 | subject) + error