Analyzing fMRI data with GLM

What is GLM?
GLM (general linear model in statistics) forms the core of the statistical pipeline in fMRI data analyses. Suppose we want to localize a brain activity associated with moving the left arm. And suppose we utilize a simple blocked design consisting of two alternate blocks of rest and move. We can predict the parts of the brain associated with moving the left arm through a simple prediction. The area in which BOLD signal changes in the same manner as our task-block can be the focus. In essence, what GLM does is multiple regression, i.e. predicting an outcome based on several predictors (regressors, a.k.a independent variables, or explanatory variables). In fMRI, the outcome Y is the voxel time series and this operation is done for each voxel in the brain.

So, what is GLM? There are three basic meanings:

    a)  It is a model to predict or estimate the actual BOLD response based on a known task design.

    b)  It is 'linear' because our prediction is formed through a linear combination of many factors. 

    c)  It is 'general' because we use assumptions to proceed with the statistical tests (t-test, F-test).

Basic statistics remain the same. The variable that we want to predict is the dependent variable. In GLM, the BOLD response of each voxel in the brain is our dependent variable, or to be more precise, it is the voxel time series. The predictors or regressors are also called EVs or explanatory variables. GLM tries to model each voxel's time series as a linear combination of many EVs. 

In other words:

    a)  GLM is a model-fitting strategy done one by one for every voxel, a mass-univariate statistics.

    b)  The contribution of each regressor is called parameter estimate (PE) or β coefficient.
    c)  As in normal regression, the difference between the fitted and actual data is called the residual.


Fundamentally, GLM is a simple linear regression, that is, a model that consists of a predictor (explanatory variable) and a coefficient (slope or β-value). The following diagram sums up what happens when you construct a GLM when you analyze your fMRI data and how it relates to the subsequent calculation of the t-value. A good source for this is the FSL Course slides.

GLM in a Diagram

Fig-1: An illustration of how General Linear Modeling of functional MRI data analysis works.

Key Steps in GLM:

• Our goal is to predict the voxel time series (y) during the task phase using our predicted BOLD response model. Voxel that best predicts the observed time series should be responsible for task performance.
• We assume we roughly know the shape of the BOLD response based on our experimental paradigm. The approach is to employ a canonical HRF or hemodynamic response function such as the double-gamma function. Convolve this with the boxcar function to obtain the model of the BOLD response that corresponds to our blocked design. Note, this function is also called the basis function because you want to express multivariate data using a certain function.
• Determine the whole set of regressors, EVs, or predictors. These include the canonical model mentioned above, motion-related variables, all other noise signals as noise regressors, etc. A collection of EVs is also known as the design matrix. We include noise regressors to take into account artefacts corrupting the actual BOLD response (hence the term: to regress them out).
• A good design matrix should not contain an EV that is highly correlated or a linear combination of another EVs. If this happens, the model is rank-deficient. Our estimates may then fail to represent the actual BOLD response. A well-designed study/task, hence, should ensure that the regressors in the GLM are orthogonal. This topic needs a separate discussion.
• In GLM, the best prediction occurs when the residual is minimized. Thus, to find the best fit of our data is to estimate β values for each voxel such that the residual is minimum.

The solution for unbiased estimates of betas can be found by using the Gauss-Markov Theorem but with some assumptions involved: 

• The residual errors ~ N(0, σ².I),  are normally distributed and independent of each other, 
• The homogeneity of variance is held, 
• Each regressor cannot be a linear combination of other regressors. 

Unfortunately, fMRI observations contain some sort of temporal autocorrelation denoted by a matrix V, making the dirty error variance Σ_e = σ².V. To ensure the data satisfy the assumption, such autocorrelation shall be handled properly. Data preparation such as prewhitening using FILM in FSL or coloring (temporal smoothing in SPM) are meant to achieve this.

Sources of autocorrelation in fMRI dataset include neural source, scanner drifts, physiological signals (respiration and cardiac pulsation), and the head movement. Not accounting for this temporal autocorrelation results in spuriously high fMRI signal at one time point that bleed over to the subsequent time points, which increases the chance of false positive results in task-based studies. The word 'pre' means the operation is added or applied to the original data. The word 'whitening' is a term that is motivated by the fact that white light contains all visible frequencies in equal amounts, i.e. making the spectrum flatter.

What are the steps? 

• First, we fit the GLM ignoring the autocorrelation, estimate the residual errors which contain the correlated voxels and have variance Σ_e. 
• Use the obtained residual to estimate the autocorrelation matrix V. FSL performs a local voxelwise estimation (instead of the global in SPM), which can be modelled quite well with an AR(1) function. Then smoothen this estimate V to correct for any biases using a Tukey Taper, hoping to downweight noisy estimates at higher lags. 
• Once V is obtained, we find a new matrix K such that KVK' = I, where I is the identity matrix. Subsequently, prewhitening is achieved by multiplying both sides of the original GLM equation with this matrix K. 
• We then repeat the GLM again to obtain the correct beta weights. The resulting beta estimates are now unbiased! 

Hypothesis Test: Contrast?

Okay, suppose the whitened GLM is completed. What's next? We want to know which voxels are highly activated due to task performance. How? If the resulting β from the GLM is big and the residual is small, highly likely the voxel is significantly activated during the task. As a quick note: in simple linear regression, β can be treated as a slope. It represents the change in the dependent variable Y resulting from a unit change in the predictor. In statistics, we assess how good the predictors are by comparing whether the PEs are significantly different from zero. Why? If the β = 0, it means the predictor does not contribute anything to the outcome Y.

We are now in the position to perform a t-test on every voxel to know which ones carry significant βs. The formula is given below, the denominator denotes the var[c^T β]. This variance depends on the residual and the critical assumption here is that the residual follows a normal distribution ~ N(0, σ²I).

What is the # dof? On each voxel, we have n time points and n is usually large. A few minutes scan can give you 500 volumes, In our t-test, the degree of freedom = n ─ p. With such a huge # dof, our t-statistics is approximately a z-distribution. Our statistical test gives a map called "Statistical Parametric Map". Because of the number of voxels in the brain, and knowing in reality that each voxel isn't that independent, the problem of multiple comparisons is challenging. Another post will discuss how this can be managed, e.g. through Random Field Theory.

The null hypothesis that a particular voxel isn't significantly related to our task is H₀: c^T β  = 0, where c represents a set of contrast. This contrast is to test the β values and serves as a linear combination of regressors. For p-number of regressors (EVs), the contrast can be written as a p × 1 matrix, each is associated with a beta. In our blocked design, the first regressor is associated with the predicted BOLD response. This allows us to put the first contrast as 1, leaving the rest of unwanted regressors as zeros, c = [1 0 0 ... 0]. This does not look like a normal contrast but as long as there is no singularity in the GLM, the hypothesis test is deemed valid.

Another example, say, if you use two different and independent tasks, e.g. visual and auditory tasks, you want to localize voxels related to your visual, auditory, or both on average. You may then set the contrast to be c₁ = [1 0 0 ... 0], then c₂ = [0 1 0 ... 0], c₃ = [1 1 0 ... 0]. If you want to localize voxels that are significantly more active in the visual than in the auditory task, then c₄ = [1 -1 0 ... 0]. The importance of assigning correct contrasts should not be ignored. Note that we can also perform an F-test, just like ANOVA, to test whether any of the contrast is significant.

Fig-2: An example of a rank-deficient GLM. The design matrix contains two task-related regressors from a blocked experiment (no noise regressors for simplicity). We can equally well use either EV1 or EV2 to explain what we see in the data. In other words, Y can be fit with any linear combination of the EVs, as long as β1 + β2 = 0.9. This yields a problem in hypothesis tests with [1 0] or [0 1] contrasts, but we can still get away with that when testing [1 1] contrast. The computation may still be successful, but your inference hereafter won't be accurate. 

For more detailed explanation on GLM, I find the following sources useful:

      [1]   A nice review article by Martin Monti (Front. Hum. Neurosci., 2011).      
      [2]   Chapter 5, in "Statistical Analysis of fMRI Data", a book by Ashby FG (2011).
      [3]   A more classic Chapter 9, in "Functional MRI: An introduction to methods" by K. Worsley (2001).
      [4]   For this post, I consulted Chapter 7, in "fMRI Techniques and Protocols", Woolrich M, et al. (2009).

Voxel vs Cluster-based Thresholding
In general, the main statistical analysis involves performing GLM at every voxel of the brain. That's why it is also called the mass univariate analysis. Once the statistical test is carried out, we get a z-map or t-map with a large #dof. Each voxel is constructed under the null hypothesis that nothing interesting happens because of our behavioral task. In statistics, the most common (frequentist) approach of a hypothesis testing is to obtain a p-value which is then compared against a certain threshold α. This α also represents the chance of a false positive which is capped at 0.05. Simply put, when we have data with 100,000 voxels, we have 5,000 voxels deemed false positive when we use α = 0.05. The error or bias in drawing a conclusion due to multiple comparison problems is also called familywise error because we are essentially repeating the same t-test for all voxels.

In statistics, the most common correction method for familywise error is Sidak-Bonferroni correction. In fMRI, however, the method is highly conservative. Two more popular methods are widely used in the neuroimaging community: the classic Random Field Theory or GRF (Worsley, K. 1995/1996), the False-discovery Rate or FDR (Benjamini & Hochberg, 1995). A more recent development based on non-parametric statistics using the permutation method is also popular. Detailed discussions about these methods are not presented here.

Post-stats thresholding is the last step in the pipeline to draw the conclusion about our neuroimaging data. This step is done with two objectives in mind. First, we want to know whether a certain voxel or a set of voxels is really active due to the task. Second, we draw that conclusion after having a proper correction for multiple comparisons. Technically, there are two different ways of thresholding:

(1) Voxel-wise thresholding
We can correct for multiple comparisons and do thresholding at every voxel by showing which part of the brain is active at a particular significance level. We reject the null hypothesis that there is no activation if t > u_v. This is called voxel-wise thresholding with threshold u_v. The merit of this method is the high specificity but it faces a serious multiple comparison problem. In earlier days when the scientific community was overly excited about fMRI, the results were reported using voxel-wise thresholding using a more stringent Bonferroni correction, say p < 0.0001, focusing on an area of interest, say frontal motor cortex or visual cortex only. Theoretically, Bonferroni correction is not suitable as it assumes adjacent points are not correlated.

(2) Cluster-based thresholding
Cluster-based thresholding is the more accepted way in the community currently. If a brain region is activated, we expect that not only one, but the surrounding set of voxels gets activated too. Instead of dealing with an individual voxel, we look at a set of voxels called a cluster. Cluster-based thresholding is done in two steps. We first create a binary image (pass/fail) of voxels passing the cluster-forming threshold u_k. Following this, we find another thresholding parameter that is based on either cluster size or peak voxel. For example, we set the cluster parameter with a size greater than a certain value p_k. Here, the sensitivity is typically better but worse specificity. Correction for multiple comparisons is more complicated because it involves a series of contagious voxels. This will be discussed in more detailed in my other post (Random Field Theory, RFT).

Case study: My own data
I'd like to write a simple example of how one can perform statistical analysis (GLM) in FSL. A participant performs a very simple motor task inside the scanner with the right-dominant arm. This is a blocked design with two conditions (rest | move | rest | move |.....) lasting for about 7 minutes (TR=1.69 sec). The scan recorded 250 volumes, meaning, the data would have 250 time points.

The following graphs illustrate what we expect from the GLM analysis to identify voxels that are significantly associated with the task. The red line in the top panel is our observed data, the time series of a particular voxel at a certain coordinate. Ideally, the best fit would see a smooth HRF curve that spans over the 'move' block. The purple plot of a full model fit of the actual time series given by the GLM analysis, that is, given all of the regressors, whether any activation in this voxel can be attributed to any of the task conditions. Lastly, the green line represents only the contrast of interest, and is usually only meaningful when having 3-4 task conditions or taking a simple main effect. Another way of looking at full model vs partial mode fit is this. The former shows what happens when all EVs are used in the fit, whereas the latter shows what happens when only some EVs are used (those involved in the contrast).

The resulting z-map shown at the bottom panel is thresholded at Z = 3.0 with a corrected cluster p-value < 0.05. This simply means that among the voxels that satisfy Z > 3.0, we apply the cluster-based thresholding with a correction for multiple comparisons. In FSL, GRF theory is used as the default correction method.

Fig-3: Outputs of the task-based analysis using GLM. The color map overlayed on the brain image is sometimes called the statistical parameter mapping or SPM (after Friston, et al). Voxelwise temporal autocorrelation must be taken care of by default (prewhitening stage) to ensure that the GLM outputs are valid.