Statistical Analysis of fMRI Data

This section covers the fundamentals of statistical analysis in fMRI data, including creating a design matrix, fitting a general linear model (GLM), contrasts, and multiple comparison corrections.

About This Data Set¶

This data set is a single run of a functional MRI (fMRI) experiment. A run is a continuous acquisition of fMRI data while a subject performs a task or rests. Typically runs last from 5 to 15 minutes. In this example, an adult human viewed many small dots shown in white on a black background. In alternating blocks of 20 s each, these dots either moved from the centre of the screen outward (motion condition) or remained stationary (baseline condition). During the motion blocks, the participant was instructed to press a button with their right index finger whenever they saw the dots speed up briefly (which happened at random intervals). During the stationary (baseline) blocks, the participant was instructed to press the button whenever the dots dmiied briefly (which also happened at random intervals).

The purpose of this task was to identify, or localize, brain regions involved in processing visual motion. Localizer tasks are commonly used in fMRI research to identify specific brain regions of interest (ROIs) based on specific response properties (in this case, sensitivity to visual motion). Often this will be used to test smore specific hypotheses about the response properties of a particular area, or how it might interact with other brain regions during more complex tasks. For the present example, we will focus on identifying the well-known human MT+ complex (or V5), a region in the lateral occipital cortex that is strongly selective for visual motion.

The goal of the analysis is to identify brain regions that are more active during the motion blocks compared to the baseline blocks.

fMRI Preprocessing¶

Before performing statistical analysis on fMRI data, it is essential to preprocess the data to correct for various artifacts and prepare it for analysis. Common preprocessing steps include:

• Motion Correction: Aligning all fMRI volumes to a reference volume to correct for head movements during the scan.

• Slice Timing Correction: Adjusting for differences in acquisition time between slices in a volume.

• Spatial Normalization: Transforming the fMRI data into a standard anatomical space (e.g., MNI space) to allow for group analysis and comparison across subjects.

• Smoothing: Applying a spatial filter to the data to increase the signal-to-noise ratio and account for anatomical variability across subjects.

However, for simplicity here we will use preprocessed data that has already undergone motion correction and spatial normalization using established neuroimaging software (FSL). This allows us to focus on the statistical analysis aspects without getting bogged down in preprocessing details.

Introducing BIDS¶

The Brain Imaging Data Structure (BIDS) is a standardized way to organize and describe neuroimaging data. It simplifies data sharing and analysis by providing a consistent structure. BIDS is not a file format; rather it is a set of conventions for organizing and describing data. The intent of BIDS is to make it easier for researchers to share, analyze, and reproduce neuroimaging studies. By organizing data in a standardized way, BIDS facilitates collaboration and enhances the reproducibility of neuroimaging research. It can also facilitate the use of automated analysis pipelines that expect data to be organized in a specific way. This not only makes it easier for invidual labs to manage their data and analysis pipelines, but also promotes best practices in data management across the neuroimaging community. It has become standard practice in the neuroimaging community to share data in BIDS format whenever possible, for example via platforms like OpenNeuro. OpenNeuro requires datasets to be BIDS-compliant before they can be shared on the platform. Nto only does this ensure that datasets are well-organized and easy to understand, but it also makes it easier for other researchers to apply different analysis methods to the same dataset. It also facilitates aggregation of research results; for example it would be quite easy to implement a meta-analysis across multiple studies if they were all BIDS-compliant.

BIDS Directory Structure¶

The data used in this example is organized in BIDS format. The directory structure looks like this:

data/BIDS/
├── dataset_description.json     # Dataset metadata
├── participants.tsv             # Participant information
├── participants.json            # Participant data dictionary
├── README                       # Dataset documentation
├── task-MotionLoc_bold.json    # Task-specific metadata for BOLD
├── task-MotionLoc_events.json  # Task-specific events metadata
├── code/                        # Analysis code
│   └── convert_design_to_events.py
└── sub-001/                     # Subject 001
    └── ses-01/                  # Session 01
        ├── anat/                # Anatomical data
        │   ├── sub-001_ses-01_task-MotionLoc_T1w.nii.gz  # T1-weighted image
        │   └── sub-001_ses-01_task-MotionLoc_T1w.json    # T1 metadata
        └── func/                # Functional data
            ├── sub-001_ses-01_task-MotionLoc_bold.nii.gz    # BOLD fMRI
            ├── sub-001_ses-01_task-MotionLoc_bold.json      # BOLD metadata
            └── sub-001_ses-01_task-MotionLoc_events.tsv     # Stimulus events

There are a few key features of BIDS organization to note:

• Hierarchical structure: Dataset → Subject → Session → Data type

• Standardized naming: Files follow the pattern sub-<label>_ses-<label>_task-<label>_<suffix>.<extension>

• Metadata files: JSON files contain acquisition parameters and other metadata

• Events files: TSV files describe stimulus timing and experimental conditions

• Documentation: README and dataset_description.json provide dataset-level information

In this example, we will use a BIDS-formatted dataset to demonstrate fMRI data analysis. We’ll start with the PyBIDS library for Python.

Read the BIDS Structure¶

First, we need to read the BIDS structure to access our fMRI data. We’ll use the bids library for this purpose.

Read the fMRI Data¶

Next, we will read the fMRI data from the BIDS structure. We will use the .get() method on the BIDS layout object we obtained above, to get the path to the BOLD fMRI NIfTI file, and then use the nibabel library to load the NIfTI files.

Vizualize the fMRI Data¶

It’s always good to visualize your data and do some other validation checks when you first open an image. This way, you can ensure that the data has been read correctly and is in the expected format.

Read fMRI Header Information¶

Nifti files contain header information that provides metadata about the image, such as voxel dimensions, data type, and orientation. We can access this information using the nibabel library.

View dimensions and voxel size¶

The Nifti header contains a lot of information. Let’s pull out specifically the voxel dimensions, and the number of voxels in each plane of the image.

You’ll notice in the output that the fourth dimension corresponds to the TR (repetition time) of the fMRI scan, which is the time between successive volume acquisitions. This is important for understanding the temporal resolution of your fMRI data. It also explains why fMRI data is often represented as a 4D dataset, with three spatial dimensions and one temporal dimension. In the output above, the 4th value listed in the “Number of Voxels in each plane” corresponds to the number of time points (or volumes) acquired during the fMRI scan.

Load Structural (Anatomical) Image¶

Since the anatomical image is strrctural, there is only one volume so it’s a 3D image with no time dimension. You’ll also notice that the voxel dimensions are different between the anatomical and functional images. Anatomical images typically have higher spatial resolution than functional images. Since the voxels are smaller, more of them are required to cover the same brain volume.

Plot Anatomical Image¶

Note that the anatomical image has been skull stripped to remove non-brain tissue. This is a common preprocessing step in neuroimaging to focus analyses on brain tissue only. It also helps anonymize the data by removing identifiable features such as the face and skull. As is sometimes the case, the skull stripping was not perfect, and some non-brain tissue can be seen at the top of the head in the sagittal view.

Read Events File¶

The events file contains information about the timing and type of stimuli presented during the fMRI scan. We will read this file using pandas.

This shows us that there is only one type of event: motion. In this block design, blocks of moving dot stimuli alternatied with blocks of static dot stimuli. However, since there are only two conditions (motion and static), we can model the static condition as the baseline (intercept) in our design matrix. This means we don’t need to explicitly include it as a separate regressor in our model; rather, the baseline is inferred from the absence of any explicitly-defined event occurring at a given time point.

The Design Matrix¶

The design matrix represents the experimental conditions and is used in the GLM analysis. At first it may seem confusing, but it’s essentially a way to encode when different stimuli were presented during the fMRI scan. This is used in the statistical analysis to essentially perform a correlation (regression) between the BOLD time series at each voxel in the image, and the time series that would be predicted for brain areas that are sensitive to the experimental conditions.

In our simple moving/not moving block design, this means we want to correlate the data at each voxel with a time series that is high during blocks of moving dots, and low during blocks of static dots. This is often called a boxcar design. Let’s draw this boxcar design using matplotlib and the event timings in events_df:

Modelling the Hemodynamic Response Function (HRF)¶

Although our experimental design is a simple boxcar design, the BOLD response in the brain does not follow this boxcar shape exactly. Instead, the BOLD response is delayed and smoothed due to the hemodynamic response function (HRF), which describes how blood flow changes in response to neural activity. There are several models of the HRF, but a common one is the glover model, originally described by Gary Glover in 1999. This model posits that the HRF in response to a transient (short) neural event can be described by a relatively simple mathematical function that captures the key features of the BOLD response, including its initial rise, peak, and subsequent undershoot. If we plot the Glover HRF it looks like this:

Note that the HRF does not start right at time zero (when the stimulus occurred), but has a delay of about 2 seconds before it starts to rise. This delay reflects the time it takes for capillaries to dilate in response to changes in nneural activity, and allow more blood to flow to the area. The HRF peaks around 5-6 seconds after the neural event and then gradually returns to baseline, often dipping below baseline (undershoot) before stabilizing.

Because of this delayed and smoothed response, we cannot simply use the boxcar design as our regressor in the GLM analysis. Technically we could, but the correlation to the actual BOLD response would be much weaker, since the BOLD response does not follow the boxcar shape exactly. Instead, we need to account for the HRF when creating our design matrix. We do this using a mathematical operation called convolution. Convolution essentially “smears” the boxcar design by the shape of the HRF (imagine the HRF shown above sliding in time from left to right over the boxcar plot), resulting in a regressor that more closely matches the expected BOLD response to the stimuli.

Regressors in the Design Matrix¶

Experimental Regressors¶

Whe we create the desgn matrix, we use the HRF model to convolve our time series this way. In this simple example, we only have one condition (motion), so our design matrix will have one regressor for the motion condition. However, in a more complex design with multiple conditions, we would have multiple regressors (time series), one for each condition.

Nuisance Regressors¶

The design matrix also includes other regressors, such as an intercept term (constant) to model the baseline signal, and potentially other nuisance regressors to account for confounding factors.

The intercept is necessary because the general linear model (GLM), which we use to analyze fMRI data, is an extension of the simple formula for a line that you likely learned in school: y = mx + b. In this formula, ‘b’ represents the y-intercept, which is the value of ‘y’ when ‘x’ is zero. In the context of fMRI data analysis, the intercept represents the baseline level of the BOLD signal when no experimental conditions are present (or in this case, when the dots on the screen were not moving).

In addition to the intercept, we often include other nuisance regressors in the design matrix to account for confounding factors that can influence the BOLD signal but are not of primary interest in our analysis. BOLD data is known to contain drift (meaning that the signal can slowly change over time) due to scanner instability or physiological factors. It may also contain low-frequency (slowly changing over time) noise due to physilogical processes such as respiration and heart rate. To account for these confounding factors, we can include nuisance regressors in our design matrix that model these effects. The simplest way to do this is to include polynomial trends (e.g., linear, quadratic) to model drift, and potentially other regressors derived from physiological recordings (if available) to model physiological noise.

Other examples of nuisance regressors include motion parameters (to account for head movement during the scan), and physiological signals (e.g., heart rate, respiration if they were recorded along with the fMRI data). Including these nuisance regressors helps to improve the accuracy and reliability of our statistical analysis by reducing the influence of these confounding factors on the estimated effects of our experimental conditions. However, in the present example we will keep it simple and only include the intercept and polynomial terms.

All of the regressors in our design matrix look like this:

Create the Design Matrix¶

In practice, we don’t need to generate all of the plots above or write code to implement any of the computations used to create them. This is all done by nilearn’s make_first_level_design_matrix function.

This function requires a number of inputs: the frame times (time points corresponding to each volume in the fMRI data), the events dataframe, the repetition time (TR), and the type of HRF model to use. It also allows us to specify additional parameters, such as the degree of polynomial trends to include as nuisance regressors. It also allows us to apply a high-pass filter to the data, which is an additional way to remove low-frequency noise.

Visualize the Design Matrix¶

We can plot the design matrix using conventional matplotlib functions.

One thing that sometimes trips people up when first learning about design matrices is the orientation of the matrix. Conentionally, we plot time on the x axis of most plots. However, in the plot above, time is on the y axis: each column represents a regressor (predictor variable), and each row represents a time point (observation). This is the standard orientation for design matrices in fMRI analysis. One advantage is that if the run is long, time is not squeezed to fit the width of the window.

You can think of the colours in each column as representing the height of the corresponding regressor at each time point. Essentially, we’ve taken the plot above and rotated it so that we’re viewing each line from above, rather than from the side. White represents the highest values, black the lowest, and shades of grey represent intermediate values. The last column is the intercept, which we didn’t plot in the previous figure.

Plot design matrix with Nilearn¶

Nilearn also provides a convenient function to plot the design matrix directly.

Run the GLM Analysis¶

To run the statistical analysis,we fit the design matrix to the time series from each voxel of the fMRI scan using the general linear model (GLM). To do this, we create an instance of the FirstLevelModel class from nilearn, specifying parameters including the repetition time (TR), which model of the HRF we want to use, and any smoothing to apply to the data. We then fit the model to the fMRI data and the design matrix, using the .fit() method provided by FirstLevelModel.

Compute and View Contrast¶

Recall from the design matrix that we are really only interested in the first column — the HRF-convolved model of our stimulus time series. All the other columns were nuisance regressors, meant to model variance that might otherwise make it more difficult for us to detect the effect of the stimulation.

To view the results of the stimulation, we need to create a contrast. Contrasts can include comparisons between different columns, but in this case the contrast between moving and static dots was “built in” to the first regressor, since it modelled the motion blocks as “on” periods separated by “off” (static) periods. So in the code below we use a list comprehension to find the column containing the label 'motion' (which was the label we gave it in the design matrix) to create the contrast matrix.

You can think of the contrast matrix as a way to specify which columns of the design matrix we want to focus on for our statistical test. When we define a contrast matrix, it must contain the same number of columns as the design matrix. Any columns with weights of zero are effectively ignored in the contrast — they form the baseline for the test.

In this case, we want to test whether the coefficient (beta weight) for the motion regressor is significantly greater than zero, indicating that there is a positive relationship between the motion condition and the BOLD signal. So we put a 1 in the first column of the contrast matrix (corresponding to the motion regressor), and 0s in all other columns.

If we had more than one experimental condition (for example if in some blocks, the dots appeared only in the left visual field, and in other blocks only in the right visual field), we could create contrasts to compare these conditions directly (e.g., left > right, or right > left) by assigning positive and negative weights to the corresponding columns in the contrast matrix. For example, if we had two conditions, left and right, as well as the same number of nuisance regressors as in the present design matrix, the contrast matrix for left > right would look like this:

[1 -1 0 0 0 0 0]

Plot the Contrast Results¶

We then pass the contrast matrix to the the .compute_contrast() method, asking for z scores as output.

Since z scores are standardized, they are easy to work with. In particular, any value of z is always associated with a specific value of p. Here we define a threshold of z > 3.09, which corresponds to a p < .001 (one-tailed).

Visuzlizing 1-tailed results¶

Although we used the z threshold for a 1-tailed test, the resulting map is 2-tailed, showing both positive and negative z values. The negative z values correspond to voxels where the BOLD signal was significantly lower during the motion blocks compared to the baseline (static) blocks. However, in this particular experiment, we are primarily interested in regions that show increased activity during motion, which is why we used the threshold for a 1-tailed test.

We can plot only the positive z values, thresholded at 3.09, by using the vmin and vmax parameters (kwargs) of the plot_stat_map function. vmin is the lowest value to plot, and vmax is the highest value to plot. Any values below vmin or above vmax are not shown in the plot. By setting vmin=0, we ensure that only positive z values above this threshold are displayed. As well, we can set a maximum value for the colour scale using vmax. This helps to enhance the contrast in the plot (rather than having it scale to the largest z value in the entire map, which is likely an outlier), making it easier to see the regions of interest. Setting vmax is also important if you’re comparing multiple statistical maps, as it ensures that the colour scales are consistent across the maps. If the same colour means different things in different plots, it can be misleading.

Visualize MT+ Activation Maps¶

The crosshairs in the image above are automatically positioned, which is not always helpful since the plot_stat_map function has no idea what brain area(s) we’re interested in. We can manually set the crosshair positions to be over left and right MT+ using known coordinates from the literature. Since it ends up being a fair bit of code, and we’re going to generate several plots like this, let’s put the code into a function called plot_stat_maps(), which takes as input the z-map we computed above.

Now we can use that to plot the z-map we computed above.

We can see that there is strong activation in both left and right MT+, as expected. The activation is a bit more extensive in the right hemisphere for this particular subject, but that can vary across individuals. There’s also a large area. of activation medially, possibly corresponding to the another motion-sensitive area, V3A.

Multiple Comparisons Correction¶

Recall that we fit the design matrix to the time series from each voxel of the fMRI scan using the general linear model (GLM). There are thousands or even millions of voxels in a typical fMRI scan, as we can see by computing the product (multplying) the dimensions of the z_map_motion image:

# compute the number of voxels in z_map_motion
n_voxels = np.prod(z_map_motion.shape)
print(f"Number of voxels in z_map_motion: {n_voxels}")

Number of voxels in z_map_motion: 10485760

So we’ve performed over ten million statistical tests! This raises an important issue known as the multiple comparisons problem. When we apply a statistical threshold such as p < .001 to the results from a single voxel, we are estimating that there is only a 1/1000 chance of incorrectly rejecting the null hypothesis (i.e., identifying that voxel as significantly active when it is not). However, we are performing far more than 1000 tests simultaneously, as we have over ten million voxels in this example. This means that a simple threshold of p < .001 is not sufficient to control the overall error rate across all voxels. Indeed, if you look closely at the image above you can see small areas of purple (which are the weakest “significant” activations that just barely pass the threshold) scattered throughout the brain, many of which are likely to be false positives.

When viewing the results of voxel-wise statistical tests across the entire brain, it is important to apply multiple comparisons correction methods to control the overall error rate, often refrred to as the family-wise error rate. There are many different ways to do this. Common methods include Bonferroni correction, False Discovery Rate (FDR), and cluster-based thresholding. These methods help to reduce the likelihood of false positives and improve the reliability of the results.

Bonferroni Correction¶

Bonferroni correction is a straightforward method for controlling the family-wise error rate when performing multiple statistical tests. The basic idea is to adjust the significance threshold for each individual test based on the total number of tests being performed. Specifically, if you want to maintain an overall significance level of α (e.g., 0.05) across all tests, you would divide α by the number of tests (N) to obtain a new, more stringent threshold for each individual test: α_corrected = α / N.

In the present data set, we have over ten million voxels, so if we wanted to maintain an overall significance level of 0.05 across all voxels, we would use a corrected threshold of 0.05 / 10,485,760 ≈ 4.77e-09 for each individual voxel. This is a very stringent threshold, and it may result in very few voxels being identified as significantly active. Thus Bonferroni correction is often considered too conservative for fMRI data, where the number of tests (voxels) is very large and the tests are not independent due to spatial correlations in the data.

In the present case, our activation is very robust and so as you can see below, even with Bonferroni correction we still see strong activation in MT+. But if you compare the plot below to the previous one, you can see that the clusters got smaller and more distinct from each other.

False Discovery Rate (FDR)¶

The False Discovery Rate (FDR) is a statistical method used to correct for multiple comparisons when performing voxel-wise analyses in fMRI data. It was first described by Benjamini and Hochberg (1995) and was introduced to the fMRI community by Genovese, Lazar, and Nichols (2002). FDR has since become a widely used approach in neuroimaging research.

Unlike Bonferroni correction, which controls the family-wise error rate (the probability of making at least one Type I error across all tests), FDR controls the expected proportion of false positives among the voxels that are declared significant. This makes FDR less conservative than Bonferroni correction, allowing for (in theory) more true positives to be detected while still controlling the rate of false discoveries. Indeed, comparing the FDR-corrected results below to the Bonferroni-corrected results above, we can see that more voxels are identified as significantly active with FDR correction.

View a table of the results¶

While the thresholded statistical maps above may look pretty, we often want to report the results in a more systamtic way, such as in a table. This makes it easier to compare results between different studies rather than “eyeballing” similarities. Obviously, however, we don’t want to make a table listing details of every voxel that survived multiple comparisons correction! Instead, it’s useful to apply clustering to find groups of contigous (adjacent) voxels all above threshold, which could be inferred to represent a specific functional area in the brain.

Nilearn provides a convenient function to extract and view clusters of activation from the statistical maps. This function identifies clusters of contiguous voxels that exceed a specified statistical threshold and exceed a specified cluster size. This is an alternative approach to multiple comparisons correction, which estimates the probability that a cluster (a set of voxels all touching each other and all above the z threshold) of a certain size would be observed by chance. This is based on the recognition that it’s less probable that 1000 voxels all z > 3.09 would be clutered together by chance, than 10 voxels all above that threshold.

However, the implementation of this in Nilearn is relatively simple, in that the user must provide a cluster size threshold. This makes it somewhat arbitrary as to how the clusters are determined, rather than defining the cluster size threshold in a principled mathematical way. Other software package use methods such as Gaussian random field theory or nonparametric permutation to estimate cluster thresholds in principled ways.

Regardless, we can use this function to extract clusters from our thresholded z map and view them in a table, below.

The table shows that the clustering algorithm only identified two distinct clusters. In the maps above, it looked like there were at least three distinct areas of activation (left MT+, right MT+, and medial occipital cortex). However, it sometimes happens in fMRI analysis that some areas of activation “bleed into” others. That is likely the case here — the medial cluster is probably connected to one of the MT+ clusters. Note that each cluster has a number of sub-clusters listed (e.g., 1a, 1b, etc.) within this larger cluster. These represent smaller peaks within the larger connected blob of significatnt voxels. Based on our simple explanation of MNI coordinates above, we can infer that cluster 1 is probably located medially (since the x values are small numbers, close to the centre of the brain), while cluster 2 is probably the left MT+ (since its x coordinates are negative). Cluster 1c has an x value of 44 which suggests it may correspond to right MT+. However, without anatomicl labels it’s a bit hard to interpret this table.

Get Anatomical Labels for Clusters¶

We can use the atlasquery Python library to get anatomical labels for the clusters identified above. This library allows us to query various brain atlases to obtain anatomical labels based on MNI coordinates. We will use the Harvard-Oxford atlas, which is a commonly used anatomical atlas in neuroimaging.

The last column of the table now contains anatomical labels. We can see from these labels that the largest cluster is centered in the intracalcarine cortex. This is the large blob we observed earlier in the middle of the occipital lobe, likely corresponding to area V3A. However, as we guessed earlier, sub-cluster 1c is located in the right lateral occipital cortex, which likely corresponds to right MT+. Cluster 2 is located in the left lateral occipital cortex, corresponding to left MT+. Thus, our analysis has successfully identified both left and right MT+ as being significantly more active during blocks of moving dots compared to static dots, as expected based on prior literature.