fMRI Tutorials

5.1 NiLearn fMRI Basic Analysis

6 min read

Analysis Overview

This tutorial implements the standard 3-level fMRI analysis hierarchy:

  1. Run-level (1st level): GLM for each individual run → COPEs and VARCOPEs
  2. Subject-level (fixed-effects): Combine runs within subjects → Fixed-effects COPEs
  3. Group-level (2nd level/random-effects): Model population using subject-level COPEs → Group statistics

Key principle: We propagate effect sizes (COPEs), not standardized statistics (z-maps), up the hierarchy. This preserves effect magnitudes and allows proper random-effects modeling across subjects.

fixed effects vs random effects — when to use which

  • Within subject (runs → subject): combine runs of the same task within one subject using fixed effects. You’re not generalising beyond this person, so pooling like this maximises statistical power. Use Nilearn’s compute_fixed_effects or stack runs into one design matrix.
  • Across subjects (subjects → group): use random / mixed effects. The goal here is generalisation to the population, so you must include between-subject variance. Nilearn’s SecondLevelModel performs OLS at the group level (equivalent to summary-statistics RFX). FSL FLAME 1 does proper mixed-effects MCMC (see 4.2 High-level Analysis with FSL). Treat group inference results as conditional on the sample size — n ≥ 20 is a common pragmatic minimum.

0  Prerequisites

RequirementWhy you need it
Python ≥ 3.10 with nilearn 0.10.x (or newer stable), numpy, pandas, scipy, nibabel, joblibCore libraries
fMRIPrep-pre-processed data (BIDS derivatives)Gives you: *_desc-preproc_bold.nii.gz, *_desc-brain_mask.nii.gz, *_events.tsv, *_confounds_timeseries.tsv
An MNI brain mask (e.g. MNI152_T1_2mm_brain_mask.nii.gz)For group-level masking

Install everything (one-off):

pip install nilearn nibabel scipy pandas joblib

1  Per-run first-level GLM

import os, glob, json, pandas as pd
from nilearn.glm.first_level import FirstLevelModel
 
# Paths come from environment variables — set BIDS_DIR / DERIVS_DIR / FMRIPREP_DIR
# before running, e.g.:  export BIDS_DIR=/path/to/ds000003
bids_dir     = os.environ["BIDS_DIR"]
derivs_dir   = os.getenv("DERIVS_DIR",   f"{bids_dir}/derivatives")
fmriprep_dir = os.getenv("FMRIPREP_DIR", f"{derivs_dir}/fmriprep")
results_dir  = os.getenv("RESULTS_DIR",  f"{derivs_dir}/level1_nilearn")
 
subjects    = sorted([d.split('-')[1] for d in glob.glob(f'{bids_dir}/sub-*')])
runs        = ['1', '2']
contrasts   = {
    'Congruent'   : 'congruent',
    'Incongruent' : 'incongruent',
    'Incongr_gt_Cong': 'incongruent - congruent'
}
motion_cols = ['rot_x','rot_y','rot_z','trans_x','trans_y','trans_z']
 
for sub in subjects:
    for run in runs:
        # ----- file paths ---------------------------------------------------
        func   = f'{fmriprep_dir}/sub-{sub}/func/' \
                 f'sub-{sub}_task-flanker_run-{run}_space-MNI152NLin2009cAsym_desc-preproc_bold.nii.gz'
        mask   = func.replace('preproc_bold', 'brain_mask')
        events = f'{bids_dir}/sub-{sub}/func/' \
                 f'sub-{sub}_task-flanker_run-{run}_events.tsv'
        confounds_tsv = f'{fmriprep_dir}/sub-{sub}/func/' \
                        f'sub-{sub}_task-flanker_run-{run}_desc-confounds_timeseries.tsv'
 
        events_df    = pd.read_csv(events, sep='\t')
        confounds_df = pd.read_csv(confounds_tsv, sep='\t')[motion_cols]
 
        # ----- read TR from the BIDS sidecar (portable across datasets) -----
        sidecar = f'{bids_dir}/sub-{sub}/func/' \
                  f'sub-{sub}_task-flanker_run-{run}_bold.json'
        with open(sidecar) as f:
            t_r = json.load(f)['RepetitionTime']
 
        # ----- model --------------------------------------------------------
        model = FirstLevelModel(
            t_r=t_r, smoothing_fwhm=6.0, mask_img=mask,
            drift_model='cosine', high_pass=1/128, noise_model='ar1'
        ).fit(func, events=events_df, confounds=confounds_df)
 
        out_dir = f'{results_dir}/sub-{sub}/run-{run}'
        os.makedirs(out_dir, exist_ok=True)
 
        for name, expr in contrasts.items():
            model.compute_contrast(expr, output_type='effect_size'   ).to_filename(f'{out_dir}/{name}_cope.nii.gz')
            model.compute_contrast(expr, output_type='effect_variance').to_filename(f'{out_dir}/{name}_varcope.nii.gz')
            model.compute_contrast(expr, output_type='z_score'       ).to_filename(f'{out_dir}/{name}_zmap.nii.gz')

read TR from the sidecar, not from memory

Reading RepetitionTime from the BIDS *_bold.json sidecar makes this script portable across datasets without manual edits. The moment you hard-code t_r=2.0, the next dataset (TR = 1.5 s, 0.8 s, multi-band, …) silently misfits the HRF and your design matrix. One line, no excuse.

What happened?

  1. Design matrix construction: Events are convolved with HRF, combined with motion regressors
  2. Model fitting: One call to fit() performs design-matrix creation, HRF-convolution, and GLS estimation
  3. Contrast computation: For each contrast, we compute:
    • COPE (effect_size): The contrast estimate (c’β) - this is what we need for group analysis
    • VARCOPE (effect_variance): The variance of the contrast estimate - needed for fixed-effects
    • Z-map (z_score): Standardized statistic for this run only (COPE/sqrt(VARCOPE))

2  Within-subject fixed-effects (combine runs)

from nilearn.glm import compute_fixed_effects
 
subject_copes = {}    # {sub: {contrast: fixed-effects cope img}}
 
for sub in subjects:
    subject_copes[sub] = {}
    for c in contrasts:
        copes   = [f'{results_dir}/sub-{sub}/run-{r}/{c}_cope.nii.gz'    for r in runs]
        varcps  = [f'{results_dir}/sub-{sub}/run-{r}/{c}_varcope.nii.gz' for r in runs]
 
        fx_cope, fx_var, _, fx_z = compute_fixed_effects(
            copes, varcps, precision_weighted=True, return_z_score=True
        )
        
        # Save the fixed-effects contrast estimate (effect size) for group analysis
        cope_fout = f'{results_dir}/sub-{sub}/{c}_fixed_cope.nii.gz'
        fx_cope.to_filename(cope_fout)
        subject_copes[sub][c] = cope_fout
        
        # Also save z-map and variance for completeness/diagnostic purposes
        fx_z.to_filename(f'{results_dir}/sub-{sub}/{c}_fixed_zmap.nii.gz')
        fx_var.to_filename(f'{results_dir}/sub-{sub}/{c}_fixed_varcope.nii.gz')

What happened?

  • compute_fixed_effects does an inverse-variance weighted average of runs
  • Key point: We save the fixed-effects contrast estimates (COPEs) for group analysis, not z-maps
  • Z-maps represent standardized effects specific to each subject’s noise level and are not appropriate for group-level random-effects analysis

3  Group (second-level) random-effects

from nilearn.glm.second_level import SecondLevelModel
from nilearn.datasets import load_mni152_brain_mask
import pandas as pd
 
group_dir = f'{results_dir}/group'
os.makedirs(group_dir, exist_ok=True)
 
# Our fMRIPrep derivatives are in MNI152NLin2009cAsym. Passing mask_img=None lets
# nilearn compute the group mask from the data itself; here we use an explicit
# MNI brain mask (nilearn resamples it to the data grid) so results are reproducible.
brain_mask = load_mni152_brain_mask()
 
for c in contrasts:
    # Use contrast estimates (effect sizes), NOT z-maps, for group analysis
    inputs = [subject_copes[sub][c] for sub in subjects]
    design = pd.DataFrame({'intercept': [1]*len(inputs)})
 
    model  = SecondLevelModel(smoothing_fwhm=6.0, mask_img=brain_mask)
    model  = model.fit(inputs, design)
 
    # Generate group-level statistics
    cope   = model.compute_contrast('intercept', output_type='effect_size')
    zmap   = model.compute_contrast('intercept', output_type='z_score')
    tmap   = model.compute_contrast('intercept', output_type='stat')
    
    # Save outputs
    cope.to_filename(f'{group_dir}/{c}_group_cope.nii.gz')
    zmap_fn= f'{group_dir}/{c}_group_zmap.nii.gz'
    zmap.to_filename(zmap_fn)
    tmap.to_filename(f'{group_dir}/{c}_group_tmap.nii.gz')

Why use contrast estimates (COPEs) instead of z-maps?

  • Random-effects analysis requires effect sizes that can be averaged across subjects
  • Z-maps are standardized by each subject’s individual error variance - they don’t preserve the actual magnitude of effects
  • COPEs represent the actual c’β (contrast times beta) values that can be meaningfully compared across subjects

4  Cluster-level inference + table

from nilearn.glm import cluster_level_inference
from nilearn.reporting import get_clusters_table
from scipy.stats import norm
 
z_thresh  = norm.isf(0.001)                        # voxel p<0.001
zmap_file = f'{group_dir}/Incongr_gt_Cong_group_zmap.nii.gz'
 
# A. All-Resolutions Inference (Rosenblatt et al. 2018): for each supra-threshold
#    cluster it returns a map of the *proportion of voxels that are truly active*
#    (FWE-controlled) — NOT a cluster-FWE-thresholded z-map.
proportion_active = cluster_level_inference(
    zmap_file, mask_img=brain_mask,
    threshold=z_thresh, alpha=0.05
)
proportion_active.to_filename(f'{group_dir}/Incongr_gt_Cong_ARI_propmap.nii.gz')
 
# B. cluster table (≥20 vox @ voxel p<0.001)
table = get_clusters_table(zmap_file, stat_threshold=z_thresh, cluster_threshold=20)
print(table.head())

What cluster_level_inference actually returns

This is All-Resolutions Inference (ARI), not classical cluster-FWE thresholding. Each cluster’s voxels are assigned the guaranteed proportion of true positives it contains (e.g. 0.6 = “at least 60% of this cluster is genuine signal”), with FWE control across the whole analysis. For a conventional cluster-extent-thresholded map, use threshold_stats_img(..., height_control='fpr', cluster_threshold=...) — see 5.3 Visualization and Statistical Thresholding with Nilearn.

5  Visualisation

from nilearn import plotting
from nilearn.image import index_img
 
plotting.plot_stat_map(
	proportion_active,
	threshold=0.0, display_mode="z", vmax=1,
	cmap="inferno",
	title="Incongruent > Congruent: proportion of true positives (ARI)",
)
 
plotting.plot_stat_map(
	zmap_file,	
	threshold=z_thresh,
	display_mode="z",
	title="Incongruent > Congruent (uncorrected p < 0.001)",
)

6  Parallelising (optional)

from joblib import Parallel, delayed
Parallel(n_jobs=6)(
    delayed(process_run)(sub, run, ...)  # wrap the first-level code in a function
    for sub in subjects for run in runs
)

7  Summary

StepNilearn function(s)What it replaces in FSL/SPM/AFNI
1FirstLevelModelFEAT (FSL), 3dDeconvolve (AFNI), SPM first-level
2compute_fixed_effectsFEAT “higher-level” fixed-effects
3SecondLevelModelFLAME, 3dttest++, SPM second-level
4cluster_level_inference, get_clusters_tablecluster, randomise, 3dClustSim
5plot_glass_brain, view_img_on_surfFSLeyes, AFNI viewer, SPM GUI

Further reading

Graph View

Backlinks