Additional Notes on fMRI

TR and TR in MRI

When deciding a suitable MRI sequence, you have to decide the choice between short or long repetition time (TR) and echo time (TE). Different combinations of these two values produce different types of MR images. See below. To recap, the T2* image has transverse magnetization which decays much faster than would be predicted by natural atomic and molecular mechanisms. In practice, T2* can be considered as an effective/observed T2 image, whereas the T2 image is considered the natural or true T2 (i.e. T2* is always less than or equal to T2). The T2* images result principally from inhomogeneities in the main magnetic field.

Statistical power in MRI

Statistical power is defined as the probability of rejecting the null hypothesis when it is false. In most cases, it is dependent on three main factors: (1) the desired effect size; (2) alpha value or Type I error; (3) the sample size tested. Of these 3 values, the sample size is the one that is within the control of the experimenter (Desmond & Glover, 2002). Unlike in usual behavioural tasks, fMRI data of each voxel is typically normalized as percent signal % change between the Task and Rest conditions. Variability in the fMRI data consists of the within-scan and between-scan variability. Increasing scan time and the number of subjects will potentially reduce these two kinds of variation. 

The concept of statistical power in fMRI data is complex because it involves different expected activation patterns, analysis pipelines, and also scanning parameters. Loosely, some rule of thumb to obtain a decent signal, alpha = 0.05 and 80% power. For example:

a)  Total N = 12 - 20 subjects per condition.

b)  Scan time per subject = 10 - 15 min.

c)  Scan volume (time points) = 400 - 500.

Non-parametric Approach to MRI

Like the usual inferential statistics, we want to make a conclusion in our neuroimaging data that apply to the population. As the name implies, non-parametric tests do not depend on the nature of the distribution of the fMRI data (e.g. normality assumption of the residual component). Through resampling, we construct an empirical distribution from the data itself that can accurately estimate the probability of a particular observance. This is the fundamental idea of non-parametric statistics. Some of these approaches have been clearly discussed in Nichols & Holmes, 2001. The use of non-parametric statistics if the data is not normal has been discussed briefly in this blog previously. A more thorough review article that includes the history and logic behind non-parametric tests is by Winkler A, et al. (2014) in the NeuroImage journal.

Basic Concepts
To see why this makes sense, the task-based block design paradigm will be used, i.e. a subject performs a task in one block and rest in the other. If a specific location of the brain is used during tasks, we should expect a task-related activation during the task, but not during the rest period. Here, we reshuffle or scramble the data across TRs. If a particular voxel is involved in the task, we should not expect any strong correlation between the predicted BOLD response and the actual activation because the 'reshuffled data' is just random noise. Reshuffling does not change the mean and variance of the dataset, that is, we don't tamper with the data.

In analyzing task-based fMRI data, for example, you form a design matrix according to your behavioural task and conduct a GLM (using FEAT, let say) to obtain parametric estimates or PEs of the corresponding regressors. To correct for the familywise error, the multiple comparison problem across voxels, we then use GRF theory by defining cluster-forming threshold (Z = 2.30 in FEAT). According to Echlund et al, GRF performance in controlling Type-I Error is poor. In fact, it is recommended to use non-parametric tests such as permutation testing to control for Type-I Error.

Using randomise on my dataset
For my iMac machine purchased in 2011 with Intel i5 Quad-core, 3.1 GHz, it takes me 8-9 minutes to analyze a decent task-based, 100 x 100 x 63 x 250, preprocessed 4D NIFTI file using 1000x permutation. With FEAT-glm using the same GLM design matrix, it takes me 2-3 minutes.

Registration with ANTs

When I first embarked in the fMRI project, I didn't learn about other software package. Recently, I managed to compare the performance of T1 -> MNI normalization done with 3 different popular software package: ANTs, FSL, and AFNI. The performance was measured using "similarity metric" that I took from Nipype, a python-based neuroimaging package. My data came from MRI scans of stroke patients, as part of my project at the Jewish General Hospital. As a reference image, I took the MNI-152 with 1mm resolution. The left panel below shows what happened when I used MNI-152 2mm template instead. In both FSL_1 and FSL_2, I used the 2 mm template to obtain the non-linear transformation matrix. However, in FSL_1, I explicitly specified the 1mm template when I applied the non-linear warping (with applywarp) process.

There are two main types of cost function: intra-modal (least squares and normalised correlation) and inter-modal (correlation ratio and mutual information-based options).


The left panel shows the registration accuracy of the subject's T1 image to the standard MNI 1mm template of 3 different software packages. The right panel also shows what happened when the MNI 2mm template was used instead.

As a side note, the MNI template has been widely used in various neuroimaging software packages. It was historically produced in a report for the International Consortium for Brain Mapping (ICBM) back in 2001. The atlas was based on linear mapping of 152 healthy adults ranging from 18-44 years old, mostly white people. It has gone through 2 major revisions lately to increase the accuracy: the first one was in 2006 that made use of a non-linear mapping, they call this the MNI-152 6th generation. FSL uses this template. Another revision was the most recent one in 2009.

Recent development in functional MRI

More recent contributions in functional MRI have been based on graph theory (Bullmore & Sporns, 2009), which underlines the importance of seeing the brain as complex networks. This analysis of MRI data has been widely used in both healthy and dysfunctional brain networks. 

It is important for understanding normal brain development and functions, the networks involved in affect regulation, motor control and execution, and learning and memory (e.g. motor learning study, see Sami and Miall, 2013). In addition, the graph theory has been beneficial to understand the implication resulting from e.g. Alzheimer's disease, schizophrenia, and stroke. For example, Buckner et al. found a correlation between the site of the targeted regions and the location of major hubs in Alzheimer’s disease (2009). He et al. (2007) demonstrated how the frontoparietal network is implicated in the spatial neglect of stroke patients.

In order to establish a brain network using modern graph theory, there are a number of steps to be taken: 

- define the network nodes (usually from resting-state fMRI), 

- estimate association/correlation between nodes, 

- compile pairwise associations between nodes and generate an association matrix, 

- and, calculate the network characteristics. 

Yet, despite this emerging literature, many researchers using graph theoretical approaches for fMRI data have not closely examined whether their data violate the assumptions of graph-theory. One important assumption that has to be met is the assumption of stationarity. Time series or time-varying data is said to fulfill the stationarity principle if the statistical properties of a time series are time-invariant, meaning, the mean and variance do not change across time. Generally speaking, stationarity implies that the statistic parameter of interest does not change over time. Such assumption is also essential for analyzing fMRI time series in the frequency domain, as the Fourier transform is suitable for stationarity. For more than a decade, most scientists have always assumed that fMRI signals are stationary, however, one may refer to some recent works that warrant doing fMRI analyses with caution (e.g. Muhei-aldin et al., 2014; Ou et al, 2014).

Experimental Design - Statistical Tests

Basic inferential statistics.   I decided to separate this from the earlier entry. We shouldn't forget that what we observe is based on a sample, not a population. Inferential statistics has its goal in making a conclusion about a population from just this sample. In doing so, there bound to be some errors as our measure is just an estimate. Suppose we have a population of interest and that we draw all possible samples of size n. And suppose we compute a statistic e.g., a mean, proportion, standard deviation for each sample. The probability distribution of this statistic for all possible samples is called a sampling distribution; the standard deviation of which is called standard error which roughly tells us the spread around the mean of the sampling distribution. Mathematically, statisticians have to make some assumptions in order to model the problem accurately. This will be discussed later on.

How can we infer about a population parameter by using our sample? We use the concept of the sampling distribution (of means) to see where the true population mean lies w.r.t our sample mean. Having obtained the sample mean and standard deviation, we can find a confidence interval within which the true population parameter lies. This interval can be found with a certain confidence level, say 95% or 99%, with the help of a transformation. This data transformation converts our sample into a known sampling distribution, e.g. using a Student t-distribution, by first converting our scores into a t-value. The conclusion of our experiment is obtained by making use of a known sampling distribution and putting that distribution in a certain hypothesis test. The type of statistical test used (e.g. one-sample t-test, ANOVA, or chi-square test) depends on the nature of our variables.

In constructing a hypothesis, there are by definition the null and alternative hypothesis, indicated as H₀ and H₁. In the language of experimental design, the H₀ states that the IV or treatment (marijuana) has no effect on DV (exam score). Conversely, the H₁ states that the treatment has an effect on DV. To do, we first create a sampling distribution assuming the null, H₀ is true. We decide an alpha level α or significant threshold to represent a concept of "very unlikely-ness" or chance that we observe a treatment effect when H₀ is true. We then compare the likelihood of getting our sample statistics p to the alpha level that is usually fixed at 0.05 or 0.01. The final outcome of the test is whether we disapprove or retain the null hypothesis. For example,

• The probability of our sample statistics of 0.01 means there is only 1 out of 100 chance you find no treatment effect in your data when the null H₀ is true  ..... We reject H₀.
• The probability of our sample statistics of 0.80 means there is 8 out of 10 chance you find no treatment effect in your data when the null H₀ is true  ..... We retain H₀.
• Rejecting the null while in reality it's true is a Type I error with a chance α. Accepting the null while in reality it's false is a Type II error with a chance β. Note that these errors are not measurement errors made by mistakes but rather a measure of chance or probability.

Fig 1: Summary of hypothesis test with its associated Type I/II errors. 

Behavioral science is unique in the sense that the test must first assume the null hypothesis to be correct. Also, the test result simply makes a case without even saying which hypothesis, H₀ or H₁ is indeed true or factual. Even if H₀ is indeed true, we don't know how much the difference is with the untreated sample group. For such reason, the effect size d comes into play. The effect size is defined as the difference in means between the experimental group and control group, divided by the standard deviation of the control group.

Lastly, the ability to correctly reject the H₀ or find something significant in your data while H₀ is indeed wrong is called power, and it is simply 1 – β. A powerful design is able to detect a small effect size < 0.20 unit. Power increases or β is reduced when the alpha level increases, the mean difference between two conditions under the null hypothesis increases, and sample size n increases. The concept of power is crucial because it determines the replicability of our study. According to Cohen (1966), most of the behavioral studies have low power, i.e. the chance of replicating the same experiment yielding the same result is 50-50.

A technique called power analysis is used to determine the appropriate sample size for the study. Generally, power = 0.80 is an acceptable value. Tools such as "pwr" in R is able to help us in power analysis, e.g. find the minimum sample size to achieve certain p-value and power of 0.80.

Okay, we move on to talk about assumptions in statistics, a crucial component to make our problems mathematically solvable. Most statistical procedures such as t-test, linear regression, and ANOVA are based on assumptions. For example, a linear regression is valid when the DV and IVs are linearly related to begin with (for categorical variables such as "Yes/No", one uses logistic regression instead). Regression has other important assumptions in general, e.g. the residuals are independent, normally distributed, and their variance is homogenous. The type of test that requires certain assumptions to be met is called parametric statistics. Violating these assumptions causes incorrect inference or decision. The assumptions for parametric tests are presented below.

[Reference #1]: Among some relevant textbooks, I chose the materials from the classic JL Myers' "Fundamentals of Experimental Design" (1979).

The normality assumption.   As mentioned, statisticians make assumptions on their methodology to simplify the maths. Both t-test and ANOVA have to satisfy the following main assumptions:
(a) Samples are independent of each other; the score of Tom is independent of Ni;
(b) The normality assumptions (for ANOVA, it means the distribution within groups);
(c) The homogeneity of variance.
The normality assumption means that the sampling distribution comes from a population that is normally distributed or Gaussian. This is hard to know as we don’t have the luxury of knowing all possible samples related to our population of interest. Out of faith, we assume that if the data we take is normally distributed then the sampling distribution will also be normally distributed. Every stats guru also tells us about the central limit theorem that is, the distribution of the sample means approaches a normal distribution as the sample size N increases, regardless of the shape of the population of interest. As a rule of thumb, N > 30 is generally accepted for the assumption to be valid. The average of our sample means is itself the population mean, and the standard deviation is the standard error.

Fig 2: Comparison between normal data and its QQ-plot (left) and "not-really" normal data (right)

This can be tricky if N is as few as 15 or 20. Now, how do we test for normality?
(1) Visual comparison. Normal data has a bell-shaped distribution. It has a specific kurtosis and skewness characteristic. It is recommended to transform the kurtosis and skew values to z-scores. Another useful visual test is through using Q-Q plot that compares theoretical values against the actual or observed values. An ideal case will be a straight line with 45-deg w.r.t the horizontal. The following figure shows an example how the data that is less normal has the tendency to have Q-Q plot that is not straight, but rather curved away from the perfect straight line. Most of the times, the presence of outliers distort the normality, the skewness or kurtosis of your data. It is quite common to treat outliers by removing the particular data point or by transforming the data.

(2) The Shapiro-Wilk test. In R, you can simply write shapiro.test(variable name); or in Matlab as swtest(x, alpha). This test is super cool as it directly assesses whether a certain distribution is significantly different from the normal distribution with the same mean and standard deviation. If the p is significant (p < 0.01, say), then the data is different from the normal distribution.

Homogeneity of variance.   Suppose you have three groups in your experiment. Homogeneity of variance says that the variability of the data in all groups should be equal, because strictly speaking, they come from the same population. The same notion applies, for example, in the case of repeated measures where you test the same group in two different experimental manipulations. If the homogeneity principle holds true, the variance of the data in either case remains rightfully the same. In most cases, we usually ignore this principle and immediately perform t-test or ANOVA, etc. Refer to the following illustration.

There are two tests that you can use to assess this: Levine’s Test and Hartley’s F-test (good for large sample). Both methods test the null hypothesis that the variance of the two datasets or groups are equal. The Hartley's F-test is also called the equality of variance test, or sometimes variance ratio test since you're actually taking the ratio of. The computed F-statistics will be compared against the standard F-table according to their respective degree of freedom. In general, people have shown that ANOVA is robust against a mild non-normality and heterogenous variance (only for equal sample size among groups!).

Fig 3: Comparison between dataset with homogenous (left) and non-homogeneous variance (right)

Non-parametric statistics.    Non-parametric statistics means the methods don't rely on any assumption discussed above! In fact, the information on the distribution is obtained directly from the dataset itself. This is the most fundamental idea of non-parametric statistics! Bootstrapping, for example, is one popular non-parametric test that relies on "random sampling with replacement". This test can be used to assess stability of your statistics and it is pretty straightforward. Another test called permutation test builds, rather than assumes, the sampling distribution by 'shuffling' the observed data. Unlike bootstrapping, we use only the existing data set and the shuffles are without replacement.

I took this from Matlab website. Suppose you have a set of LSAT scores and  GPA. You can easily plot the data and compute the correlation. It is found that there is a positive relationship between LSAT and GPA, r = +0.78. Now although it may seem large, it is unknown whether this observation is statistically significant or just happens by chance. Using the bootstrp function you can resample the two vectors as many times as you like and consider the variation in the resulting correlation coefficients. If the finding is true, you would expect that bootstrapping process yields more or less the same correlation value. Refer to the following histogram. Bootstrapping gives a distribution of correlation coefficients that is centered near r = +0.78. Using permutation strategy, we can also show how the correlation found is not by chance. We permute the LSAT-GPA pairs for 1000 times, each time is followed by computing the correlation. The correlation values of the permuted data sets should collapse since the relationship has been "broken" following shuffling.

Fig 4: Distribution of correlation coefficients with bootstrapping and the location of the original r = +0.78 (red line).

During the period when I first learnt statistics, we didn't care much about such assumptions and treated them as dogma. Now, I'm learning to accept the presence of non-parametric statistical tests used when the assumptions are not met, that is, to replace the usual t-test etc. The principle behind those tests is assigning ranking. High scores will be represented by large ranks, low scores by small ranks. The analysis is then performed on the ranks rather than the original raw data. In other words, non-parametric tests are useful when dealing with ordinal data. Example of ordinal dataset includes: the satisfaction rating from 1 - 5, the order of athletes based on speed, and the school ranking.

When we are interested in comparing two independent datasets, we need to use the Mann-Whitney test or the Wilcoxon's rank-sum test. This is the counterpart of the independent t-test. In R, you can do this by:  wilcox.test(x,y,____), some parameters inside the brackets have to be defined. In Matlab, use the following instead ranksum(x,y). When we deal with repeated measures or paired sample, then we have to use the Wilcoxon signed-rank test. For example: testing the performance before and after drug consumption. You can use the same command as in (1) in R but by defining "paired = TRUE". In Matlab, however, the command will be different now, signrank(x,y). If we are interested in testing more independent groups, we have to do non-parametric version of the regular ANOVA, e.g. Kruskal-Wallis test and Friedman's ANOVA for repeated-measure designs. Other non-parametric analysis such as the Spearman's correlation coefficient (the counterpart of Pearson's r) and Kendall tau coefficient.



Final note: ANOVA can be understood from the perspective of Linear Modeling and in itself is a very huge topic, covering multiple regression, mixed-effects or hierarchical modeling, and generalized linear models. The mathematical discussions on this theme will be beyond this blog post.....

[Reference #2]: Andy Field et al, "Discovering Statistics using R" (2012).