Basic Independent Component Analysis (ICA)

What is ICA?
As an exploratory method, Independent Component Analysis (ICA) provides an alternative to seed-based resting-state fMRI analysis. Traditionally, ICA attempts to solve the cocktail party problem. The story is this. Suppose I have five people with a microphone talking at the same time. How can I separate my speaker output into voices of the first person, second, third, and so on? Indeed, ICA is a powerful method to discover hidden independent patterns or features from a set of data, thus, an exploratory analysis. Unlike GLM, it is model-free because no assumption is made on the shape/pattern of the actual BOLD response.

There are two types of ICA, Temporal and Spatial ICA. The cocktail party problem is an example of Temporal ICA. The data to be separated contains temporal information, the mixing coefficients or weights of which vary uniquely across different locations. In contrast, the dataset in the Spatial ICA contains spatial maps or networks that get activated altogether at a given time. The weights vary uniquely across time for different networks. Spatial ICA is more popular for fMRI data analysis simply because there are more voxels than time points. An image of 64 x 64 x 64 acquired for 250 volumes or TRs contains 262,144 voxels but 250 time points.

Running ICA on the resting-state data will ideally yield to a set of ICs, some of which are clearly related to activation or network, physiological processes (heart rate, breathing), or even some artefacts (e.g. motion, ghosting, slice dropout, noise, etc). Since rs-fMRI is not task-based, we don't have prior knowledge of the temporal waveform of our resting-state networks. This is where data exploration is useful. Spatial ICA attempts to split the data into a set of spatial maps, each with an associated time course. No assumption is made, no experimental paradigm to be specified. The diagram from the FSL website below summarizes this. The content below is thus based on FSL-MELODIC tool.



ICA for Resting-state fMRI
Our observation is the fMRI data. Like standard fMRI analysis, the pipeline starts with the common fMRI data preprocessing steps such as motion and slice timing correction, registration, spatial smoothing, etc. We have to also prepare the preprocessed dataset before ICA. First, simplify the 4D dataset into time vs space. Consider Y_t´m as the whole-brain fMRI image of m voxels and t number of volumes, where t < m. We then mean-center the data by removing the mean spatial map from each row of Y. We normalize the time series variance such that each column of Y will have unit variance, a step called variance normalization. Next, we remove the mean time series from each column of Y. Without variance normalisation, the PCA step below will be biased towards tissue exhibiting high temporal variability.

Fig-1: Temporal std deviation of a sample dataset together with its PCA components (FSL slides).

Two typical yet strong assumptions are made when doing ICA in FSL: (a) the observed source signals are statistically independent; and (b) they are non-Gaussian or not normally distributed. So, ICA is performed on observations that are assumed to be a linear combination of independent sources. Our fMRI image Y is a mixture of independent sources such that Y = A.s , where the A is a t ´ t mixing matrix, and s is a t ´ m spatially independent components.

The first part of ICA is carrying out dimensionality reduction to simplify the problem or computation. This can be achieved with the help of PCA decomposition using SVD. The criterion is that only components that represent a large amount of the dataset remain. The original dataset is now whitened with unit variance and reduced dimension. Rewrite the equation, Y_w = P.Y = P. A.s ; so we have the new mixing matrix A_w = P.A. Note: removing autocorrelation in the fMRI dataset is called whitening, a necessary step in the GLM parametric analyses to make the results more accurate! (Woolrich, et al., 2001). Different statistical software tools (AFNI, FSL, and SPM) have their own steps to accurately model this temporal correlation.

Because Y_w is whitened, its variance-covariance matrix Σ_w is equal to identity matrix I. In other words:  Σ_w = (Y_w Y_w^T) / (n ─ 1) =  (A_ws).(A_ws)^T/ (n ─ 1) = (A_ws s^TA_w^T) / (n ─ 1) = I. This is because the variance-covariance matrix of s.s^T = I as the components are independent to each other. Note that for any orthogonal matrix, A_w^-1 = A_w^T , making the problem solving even simpler.

So now, our unmixing matrix is the inverse of the new orthogonal matrix A_w. The job now is to predict or estimate the inverse of Â_w, such that ŝ = Â_w^-1 Y_w is maximally independent. Note that ŝ and Â_w are estimates since they are both unknown. The criteria of independence can be empirically found by using a certain algorithm, e.g. mutual information, negentropy, etc.

In dealing with fMRI data, specifically, FMRIB Oxford team proposed a probabilistic ICA model in the noisy dataset, a more robust independent component estimation where the noise properties follow is assumed to be Gaussian (Beckmann, et al., 2004). For a more detailed but readable explanation to write this post, I consult Chapter 10 in Ashby, 2011, textbook. I have been using FSL-MELODIC to perform a single-subject ICA and their website is also informative to read.

Theoretically speaking, when all ICs are added together (each one being a 4D signal formed by the outer product of the spatial map and timecourse) they equal the original data. Unlike PCA, ICA enforces independence between the components spatially, while PCA enforces orthogonality both spatially and temporally. Sometimes, the ICs may also share similar time courses. Once all the independent components have been identified, they are ordered according to the degree of importance, that is, % total variance explained.

The last step in the MELODIC is the thresholding stage that produces thresholded ICs overlayed on a background image of your choice. A threshold level of 0.5 (set by default) means that a voxel 'survives' as soon as the probability of being in the 'active' class exceeds the probability of being in the 'background' noise class. This 0.5 assumes we set an equal loss on false-positives and false-negatives. The user guide says that if instead we consider e.g. false-positives as being twice as bad as false-negatives you should change this value to 0.66.

Fig-2: ICA method to rs-fMRI data (TR = 450 msec) displays an independent component associated with the physiological signal. The thresholded spatial map illustrates the location of brain regions correlated with the heart rate (~1 Hz). The files are generated by FSL MELODIC with image cropping for practical reasons.

Fig-3: ICA method to rs-fMRI data (TR = 450 msec) displays an independent component associated with the visual and posterior parietal area (0.01 - 0.1 Hz). The files are generated by FSL MELODIC with image cropping for practical reasons.

By default, FSL MELODIC produces the following files: 

- HTML report with a logfile that collates all steps performed by melodic.

- The mask, the list of ICs, and other statistical information (smoothest : estimated smoothness, PPCA : estimated intrinsic dimensionality estimated from PPCA).

- melodic_mix: an ASCII text file that contains the estimated mixing matrix (in the noise-free case the ICA decomposition is typically written as X = A*C, where X is the original data, A is the mixing matrix, and C is the matrix containing the estimated independent components as its rows). melodic_mix contains #ICs time courses as its columns. Each time course is plotted in the IC report that melodic produces.

- melodic_FTmix: a matrix containing the power spectrum at different freq for the time courses contained in melodic_mix (plotted in the IC report under the time courses)

- Eigenvalues_adjusted : the set of eigenvalues from the initial PCA decomposition (after variance normalisation and adjusting for the dimensionality of X).
 

Group ICA Method
To assess resting-state network with ICA at higher-level or group application, group ICA (gICA) is introduced. One can apply gICA to a set of resting-state data through temporal concatenation gICA; or a set of task-based fMRI using tensor gICA. In temporal concatenation gICA, we are finding common spatial patterns across subjects with unpredictable or different time series. In tensor gICA, we assume that each subject carries a similar time-series pattern, e.g. in task-based fMRI.

Temporal concatenation is recommended for the group-level resting-state pipeline as an alternative method to seed-based analysis. Often, researchers are faced with ambiguity in seed selection and gICA can be useful, especially with the use of dual-regression or back-projection method. For a more detailed explanation, refer to the review paper by V. Calhoun's group (NeuroImage, 2011).

Smith et al., (PNAS, 2009) show that the functional connectivity networks obtained at rest from gICA method are similar to different task-based activation networks obtained from a separate analysis of thousands of subjects. This is quite an important study as the authors show that resting-state networks are not simply random fluctuation, but are recruited when subjects perform the actual task. Several functional networks are identified from resting-state data, e.g. the visual network, sensorimotor network, attention network, limbic network, cerebellar network, the default mode network (DMN), and auditory network. This so-called functional segregation shows that the brain is active even at rest.

Functional Connectivity - a brief overview

Background
During its initial application in scientific research, fMRI was used primarily with task-based behavioral studies. It was used to locate brain activation associated with a certain task or stimulus. For example, we want to study V1 in response to seeing patterns of a checkerboard, or S1 in response to a tactile stimulus on our left forearm. Since Biswal et al., (1995) found a temporal correlation in the sensorimotor area in the absence of tasks, the interest in the so-called spontaneous fluctuations in the brain has accelerated. This type of fMRI, the resting-state fMRI, is unique because the frequency content is low < 0.1 Hz. That is why it was initially thought of as mere physiological noise.

In another study, Fox et al. (2005) reveal a set of interesting brain regions where they get deactivated after a person engages in a task. These regions are known to routinely exhibit task-negative responses (deactivations) during attention-demanding tasks. For terminology sake, brain regions that get deactivated during task performance are called task-positive. When talking about spontaneous fluctuations, we have the terms correlation and anticorrelation regions. Anticorrelation signifies a form of temporal negative correlation. In that study, the fluctuation that is negatively correlated with those in the attention and executive systems is later called the default mode network.

Statistical maps produced by this spontaneous BOLD acquisition is thought to reflect functional connectivity network or resting-state network (RSN). Some authors (Seekey, 2007) used the word intrinsic connectivity to denote two distinct and dissociable networks, the "salience network," (dACC and orbital frontoinsular cortices, with connectivity to subcortical and limbic structures), and "executive-control network" that links dorsolateral frontal (dPFC) and parietal neocortices of the brain. These intrinsic networks are correlated with task-based MRI measured.

Fig-1: (A) BOLD activity during finger tapping and spontaneous fluctuation with a seed located at point a. (B) Spontaneous deactivation observed during task performance using PET. Regions that are negatively correlated with the seed a is called the default mode network.

Methodological Perspective
As mentioned, the only difference between rs-fMRI and task-based fMRI is that the acquisition is performed at rest. The subjects typically lie down in the scanner with their eyes open and looking at a cross-hair, or with eyes closed. Some differences in brain activations exist between the two variants. The minimal requirements of an rs-fMRI study make it easy to scan a wide variety of populations, subjects may get bored and fall asleep. A recent study indicates that falling asleep is the major challenge for subjects and that fixating subjects rarely fall asleep.

The way we analyze resting-state fMRI data is different from the method used in the task-based fMRI. This is because we do not have any clear task model for our design matrix. Having said that, some steps in the preprocessing pipeline are still valid here. For example, the removal of non-brain tissues and skulls, spatial smoothing, high pass filtering, image registration, and motion correction.

Resting-state fMRI suffers two essential drawbacks. First, separating wanted signals from the image artifacts produced by magnetic field distortion (e.g. near the sinuses), signal dropout, and motion can be challenging. Unless we can correct or model these events, the presence of artifacts potentially reduces the temporal correlation between two brain areas. Image distortion due to sinuses can be corrected using field map correction. FSL/FreeSurfer software package has a subroutine called BBR that helps to correct this when registering the EPI to T1-image.

Second, physiological signals such as cardiac and respiratory cycles can also modulate the wanted signals via several mechanisms (Murphy et al., 2013). More recent works have shown that, with a very high sampling rate or TR < 400 msec, we can segregate the frequency contents of spontaneous brain signals. But such MRI sequences are still work-in-progress. There are, however, a few methods to correct for this:

• Regress out the average time series obtained within the white matter region and ventricles, which is mostly of no interest.
• Regress out respiratory and heart-rate signals acquired from the actual physiological recording.
• Regress out global (whole-brain) average time series. This method is still debatable as it allegedly introduces the systematic anticorrelation areas.
• Use ICA to identify artifactual components around the perimeter, cardiac-related, or anything with the frequency > 0.1 Hz.
• Use band-pass filtering to limit the frequency content to be 0.01 - 0.1 Hz.

Once we have more or less 'clean up' our dataset, there are two most common and essential methods in resting-state data analysis, i.e. seed-based and ICA (Independent Component Analysis) methods. The methodological aspects of rs-fMRI are massive and it is always evolving. But some key points I'd like to write about RSN analyses:

• The frequency spectrum is usually between 0.01 - 0.11 Hz, i.e. the 'low frequency' spontaneous fluctuation, although a wider range up to 0.25 Hz was recently accepted.
• There are two most commonly used methods to determine RSN: the seed-based correlation analysis (SCA) and Independent Component Analysis (ICA).
• The SCA method is a model-driven analysis that requires a strong apriori region of interest selection. This method allows us to use the usual framework of general linear modeling (GLM). Example: if we are studying the effect of playing tennis on the brain, we should select M1 as one of our seeds.
• ICA is a data-driven analysis and does not require apriori assumption. ICA is also very attractive in identifying various neural activities at rest, collectively reflected as independent components, e.g. physiological artifacts and jerky head movements.
• It is important that changes in RSN do not constitute the causality of our behavioral paradigm.

Exploring Brain Organization
In several studies, ICA has been used as a robust method to segregate different functional networks of the brain (RSNs). With improvements in the MRI sequences, more and more details can be picked up using rs-fMRI techniques. The number of wanted components was found to be 10 of 25 in (Damoiseaux et al., 2006), 45 of 70 in (Smith et al., 2009), and ~23 of ~200 in (Marcus et al., 2013). Another method is more like a seed-based approach, which involves calculating pairwise correlations between all voxels and using clustering algorithms to identify groups of highly correlated voxels. Those studies identified the brain as being partitioned into somewhere around 15–20 large-scale clusters (e.g. Power et al., 2011; Yeo et al., 2011). Refer to Fig-2.

Fig-2: Various resting state networks revealing different functional areas of the brain (motor, sensory, visual etc).

Some of these networks have been known anatomically for long ago but others, such as attentional networks, are new. Interestingly, in some cases, functional connectivity grouped brain areas that were not so widely recognized as connected. With sophisticated data-driven approaches, research has been spent to refine boundaries between adjacent functional brain regions.

... Resting-state correlations grouped a set of regions and made it easier to recognize that they shared a variety of specific characteristics, bolstering the case that these regions form a functional system. The large-scale (system-level) patterns in resting-state activity, therefore serve as a useful organizing framework for interpreting results and patterns in other modalities.

Recently, a new paradigm of analyzing rs-fMRI data called the network-based approach emerges. When viewing the data as a network, less focus is given on the properties of a single brain area, but more within the larger neural system. However, it is wise to note that a network in this sense is not based on anatomical or physical neuronal networks. Instead, it is from time-varying signals arising of either BOLD (functional networks) or DTI (diffusion), which are essentially mathematical models with some underlying assumptions, that may not fully reflect the physical entity.


Source: Power, JD, Schlaggar BL, Petersen, SE (2014). "Studying Brain Organization via Spontaneous fMRI Signal", Neuron-Primer. Some other notable papers: (Raichle, 2010) for a historical and metabolic perspective; (Deco et al., 2011; Hutchison et al., 2013) for dynamical perspectives; (Bullmore and Sporns, 2012; Sporns, 2014) for network perspectives; (Murphy et al., 2013) for a methods perspective; (Lee et al., 2013) for a clinical perspective; and (Buckner et al., 2013; Craddock et al., 2013) for general perspectives.

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.

How should I do this?

I cannot catch up! My work requires me to master stuffs from three different domains: neuroscience, motor behavior, and statistics. It is simply too tiring to catch up with the amount of information that I have to learn with the time spent in writing this down online. From now on, I will just write a point of interest that I recently came across. It will be updated and revisited again and again in the future if I have more findings.

This is especially true when, for example, someone asked me about a few questions that I'm not too sure, but there is no sufficient time for me to find all of the answers. Furthermore, it is better not to retype the summary or basic theories of a certain topic. Citing the appropriate source here directly is more practical and time-efficient!

Okie, the 3rd year of my PhD officially has started :(