What it answers. Where in the brain does activity track your task design? The general linear model (GLM) is the mass-univariate workhorse of task fMRI: fit one regression per voxel, then test contrasts between conditions. Use it when you have a known design and want a statistical map of effects.
For the underlying theory, see the GLM chapters in dartbrains. This tutorial is about running the analysis in nltools.
How it works. A GLM analysis runs in two stages:
First level (single subject). Regress each voxel’s timeseries on the design matrix X → one β (and t) per regressor. A contrast is a linear combination of βs, giving a per-subject effect map.
Second level (group). Stack the per-subject contrast effect-size maps and run a one-sample test across subjects.
Feed effect sizes (βs), not first-level t-maps, into the group test: a first-level t is β / SE(β), and SE varies across subjects for reasons unrelated to the effect (scan length, motion). Stacking t-maps would conflate effect magnitude with first-level precision.
import numpy as np
from joblib import Memory
from nltools.data import BrainData, DesignMatrix
from nltools.algorithms import fdr, threshold
from nltools.utils import concatenate
# Memoize per-subject fits to disk (.cache/ is git-ignored) so re-running
# the notebook reloads results instead of refitting every voxel.
memory = Memory(".cache/tutorials", verbose=0)How to do it¶
We use the language localizer demo from nilearn — 10 subjects viewing blocks of sentences (language) vs. consonant strings (string). Each subject’s BIDS derivatives give us three files: the preprocessed BOLD, an events TSV, and a confounds TSV.
import json
from pathlib import Path
from nilearn.datasets import fetch_language_localizer_demo_dataset
from nilearn.interfaces.bids import get_bids_files
DATASET = fetch_language_localizer_demo_dataset(verbose=0)
DATA_DIR = Path(DATASET["data_dir"])
def get_sub_files(sub: str) -> dict:
"""Resolve one subject's BOLD, events, confounds, and TR from BIDS."""
derivatives = DATA_DIR / "derivatives"
sidecar = get_bids_files(
derivatives, file_tag="bold", file_type="json", sub_label=sub
)[0]
return {
"bold": get_bids_files(
derivatives, file_tag="bold", file_type="nii.gz", sub_label=sub
)[0],
"events": get_bids_files(
DATA_DIR, file_tag="events", file_type="tsv", sub_label=sub
)[0],
"confounds": get_bids_files(
derivatives, file_type="tsv", modality_folder="func", sub_label=sub
)[0],
"TR": json.loads(Path(sidecar).read_text())["RepetitionTime"],
}First level (single subject)¶
The recipe for one subject: load the BOLD (BrainData resamples to standard MNI automatically), build the design, and fit. Building a DesignMatrix from a BIDS events file creates boxcar regressors and convolves them with the canonical (Glover) HRF for you — columns come back as language_c0 / string_c0 (pass hrf_model=None for raw boxcars to .convolve() yourself). We append the motion confounds as nuisance columns and add polynomial drift. Wrapping it in memory.cache means each subject is fit once, then reloaded from disk.
@memory.cache
def first_level(sub: str, contrast: str = "language_c0 - string_c0"):
"""Fit one subject's GLM; return its design and the contrast bundle.
We return only the lightweight design and contrast maps (not the
fitted model, which carries residuals and a copy of the data) so the
on-disk cache stays small.
"""
f = get_sub_files(sub)
brain = BrainData(f["bold"])
events = DesignMatrix(f["events"], run_length=brain.shape[0], TR=f["TR"])
confounds = DesignMatrix(f["confounds"], run_length="infer", TR=f["TR"])
brain.fit(X=events.append(confounds, axis=1, as_confounds=True).add_poly(2))
return brain.design_matrix, brain.compute_contrasts(contrast, statistic="all")design, contrasts = first_level("01")
design.plot() # the design we just fit/home/runner/work/nltools/nltools/nltools/data/braindata/io.py:568: UserWarning:
Data resolution (4.500mm) doesn't exactly match template: default 3mm.
data_img = detect_and_update_mask(bd, data_img)
/home/runner/work/nltools/nltools/nltools/data/braindata/__init__.py:987: NearCollinearDesignWarning: Design matrix is nearly collinear (full rank, but ill-conditioned): column pair(s) correlated at |r| >= 0.95: RotX & Y (|r| = 0.96). The OLS betas are estimable but unstable: their variance is inflated, and small changes in the data can flip their signs or magnitudes. Nothing was dropped or modified. Possible fixes: (1) inspect the collinearity with `DesignMatrix.vif()`; (2) consider `DesignMatrix.clean()` to drop near-duplicate columns before fitting (note: which of a correlated pair survives depends on the order the design was built in); (3) try regularization — `fit(model='ridge')` keeps every regressor and shrinks correlated coefficients together.
return fit(
/home/runner/work/nltools/nltools/nltools/models/glm.py:195: RuntimeWarning: [MultiNiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used.
self._glm.fit(
/home/runner/work/_temp/uv-python-dir/cpython-3.12.10-linux-x86_64-gnu/lib/python3.12/functools.py:998: FutureWarning: residuals' is deprecated.
It will be removed in Nilearn 0.16.0.
Use 'residuals_' instead.
val = self.func(instance)
The helper returns the language > string contrast as a bundle — beta, t, z, p, se — computed in one call with statistic="all", so we can threshold the t-map here and reuse the β map for the group analysis below.
contrasts["t"].plot(
method="slices", threshold=3.09, title="sub-01: language > string (t)"
)
Even at one subject the left-lateralized fronto-temporal language network is visible (|t| > 3.09, two-tailed p ≈ 0.001).
Second level (group)¶
The same cached recipe runs per subject, returning one effect-size (β) map each. We loop over eight of the ten demo subjects.
SUBJECTS = ["01", "02", "03", "04", "05", "06", "07", "08"]
beta_maps = []
for sub in SUBJECTS:
_, sub_contrasts = first_level(sub)
beta_maps.append(sub_contrasts["beta"])/home/runner/work/nltools/nltools/nltools/data/braindata/io.py:568: UserWarning:
Data resolution (4.500mm) doesn't exactly match template: default 3mm.
data_img = detect_and_update_mask(bd, data_img)
/home/runner/work/nltools/nltools/nltools/data/braindata/__init__.py:987: NearCollinearDesignWarning: Design matrix is nearly collinear (full rank, but ill-conditioned): column pair(s) correlated at |r| >= 0.95: Z & .nl_poly_1 (|r| = 0.96). The OLS betas are estimable but unstable: their variance is inflated, and small changes in the data can flip their signs or magnitudes. Nothing was dropped or modified. Possible fixes: (1) inspect the collinearity with `DesignMatrix.vif()`; (2) consider `DesignMatrix.clean()` to drop near-duplicate columns before fitting (note: which of a correlated pair survives depends on the order the design was built in); (3) try regularization — `fit(model='ridge')` keeps every regressor and shrinks correlated coefficients together.
return fit(
/home/runner/work/nltools/nltools/nltools/models/glm.py:195: RuntimeWarning: [MultiNiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used.
self._glm.fit(
/home/runner/work/_temp/uv-python-dir/cpython-3.12.10-linux-x86_64-gnu/lib/python3.12/functools.py:998: FutureWarning: residuals' is deprecated.
It will be removed in Nilearn 0.16.0.
Use 'residuals_' instead.
val = self.func(instance)
/home/runner/work/nltools/nltools/nltools/data/braindata/io.py:568: UserWarning:
Data resolution (4.500mm) doesn't exactly match template: default 3mm.
data_img = detect_and_update_mask(bd, data_img)
/home/runner/work/nltools/nltools/nltools/models/glm.py:195: RuntimeWarning: [MultiNiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used.
self._glm.fit(
/home/runner/work/_temp/uv-python-dir/cpython-3.12.10-linux-x86_64-gnu/lib/python3.12/functools.py:998: FutureWarning: residuals' is deprecated.
It will be removed in Nilearn 0.16.0.
Use 'residuals_' instead.
val = self.func(instance)
/home/runner/work/nltools/nltools/nltools/data/braindata/io.py:568: UserWarning:
Data resolution (4.500mm) doesn't exactly match template: default 3mm.
data_img = detect_and_update_mask(bd, data_img)
/home/runner/work/nltools/nltools/nltools/data/braindata/__init__.py:987: NearCollinearDesignWarning: Design matrix is nearly collinear (full rank, but ill-conditioned): column pair(s) correlated at |r| >= 0.95: X & Y (|r| = 0.99), RotY & X (|r| = 0.99), RotY & Y (|r| = 0.97), RotY & RotZ (|r| = 0.97), RotZ & X (|r| = 0.96), ... and 1 more; the condition number of the standardized design is 35 (> 30), indicating near-linear dependence spread across several columns. The OLS betas are estimable but unstable: their variance is inflated, and small changes in the data can flip their signs or magnitudes. Nothing was dropped or modified. Possible fixes: (1) inspect the collinearity with `DesignMatrix.vif()`; (2) consider `DesignMatrix.clean()` to drop near-duplicate columns before fitting (note: which of a correlated pair survives depends on the order the design was built in); (3) try regularization — `fit(model='ridge')` keeps every regressor and shrinks correlated coefficients together.
return fit(
/home/runner/work/nltools/nltools/nltools/models/glm.py:195: RuntimeWarning: [MultiNiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used.
self._glm.fit(
/home/runner/work/_temp/uv-python-dir/cpython-3.12.10-linux-x86_64-gnu/lib/python3.12/functools.py:998: FutureWarning: residuals' is deprecated.
It will be removed in Nilearn 0.16.0.
Use 'residuals_' instead.
val = self.func(instance)
/home/runner/work/nltools/nltools/nltools/data/braindata/io.py:568: UserWarning:
Data resolution (4.500mm) doesn't exactly match template: default 3mm.
data_img = detect_and_update_mask(bd, data_img)
/home/runner/work/nltools/nltools/nltools/models/glm.py:195: RuntimeWarning: [MultiNiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used.
self._glm.fit(
/home/runner/work/_temp/uv-python-dir/cpython-3.12.10-linux-x86_64-gnu/lib/python3.12/functools.py:998: FutureWarning: residuals' is deprecated.
It will be removed in Nilearn 0.16.0.
Use 'residuals_' instead.
val = self.func(instance)
/home/runner/work/nltools/nltools/nltools/data/braindata/io.py:568: UserWarning:
Data resolution (4.500mm) doesn't exactly match template: default 3mm.
data_img = detect_and_update_mask(bd, data_img)
/home/runner/work/nltools/nltools/nltools/models/glm.py:195: RuntimeWarning: [MultiNiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used.
self._glm.fit(
/home/runner/work/_temp/uv-python-dir/cpython-3.12.10-linux-x86_64-gnu/lib/python3.12/functools.py:998: FutureWarning: residuals' is deprecated.
It will be removed in Nilearn 0.16.0.
Use 'residuals_' instead.
val = self.func(instance)
/home/runner/work/nltools/nltools/nltools/data/braindata/io.py:568: UserWarning:
Data resolution (4.500mm) doesn't exactly match template: default 3mm.
data_img = detect_and_update_mask(bd, data_img)
/home/runner/work/nltools/nltools/nltools/models/glm.py:195: RuntimeWarning: [MultiNiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used.
self._glm.fit(
/home/runner/work/_temp/uv-python-dir/cpython-3.12.10-linux-x86_64-gnu/lib/python3.12/functools.py:998: FutureWarning: residuals' is deprecated.
It will be removed in Nilearn 0.16.0.
Use 'residuals_' instead.
val = self.func(instance)
/home/runner/work/nltools/nltools/nltools/data/braindata/io.py:568: UserWarning:
Data resolution (4.500mm) doesn't exactly match template: default 3mm.
data_img = detect_and_update_mask(bd, data_img)
/home/runner/work/nltools/nltools/nltools/data/braindata/__init__.py:987: NearCollinearDesignWarning: Design matrix is nearly collinear (full rank, but ill-conditioned): column pair(s) correlated at |r| >= 0.95: X & .nl_poly_1 (|r| = 1.00), RotZ & X (|r| = 0.99), RotZ & .nl_poly_1 (|r| = 0.99), RotX & Y (|r| = 0.98), Y & .nl_poly_1 (|r| = 0.97), ... and 6 more; the condition number of the standardized design is 54 (> 30), indicating near-linear dependence spread across several columns. The OLS betas are estimable but unstable: their variance is inflated, and small changes in the data can flip their signs or magnitudes. Nothing was dropped or modified. Possible fixes: (1) inspect the collinearity with `DesignMatrix.vif()`; (2) consider `DesignMatrix.clean()` to drop near-duplicate columns before fitting (note: which of a correlated pair survives depends on the order the design was built in); (3) try regularization — `fit(model='ridge')` keeps every regressor and shrinks correlated coefficients together.
return fit(
/home/runner/work/nltools/nltools/nltools/models/glm.py:195: RuntimeWarning: [MultiNiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used.
self._glm.fit(
/home/runner/work/_temp/uv-python-dir/cpython-3.12.10-linux-x86_64-gnu/lib/python3.12/functools.py:998: FutureWarning: residuals' is deprecated.
It will be removed in Nilearn 0.16.0.
Use 'residuals_' instead.
val = self.func(instance)
concatenate stacks the per-subject maps into one (n_subjects, n_voxels) BrainData. BrainData.ttest runs a voxelwise one-sample test, returning the effect-size mean, the parametric t, a signed z, and p. nltools.algorithms.threshold keeps the z values whose p clears a cutoff — here voxelwise p < 0.001.
group = concatenate(beta_maps)
group_result = group.ttest()
group_z = threshold(group_result["z"], group_result["p"], thr=0.001)
group_z.plot(
method="slices", title="Group: language > string (voxelwise p < 0.001)"
)
Multiple-comparisons correction¶
That p < 0.001 map is uncorrected — it ignores that we ran tens of thousands of tests. nltools.algorithms.fdr returns the p-threshold controlling the false-discovery rate. Whole-brain correction is stringent: on a ten-subject demo, far fewer voxels survive than at the uncorrected threshold — exactly the inflation that correction guards against. Restricting the search to an ROI (see the MVPA tutorial) recovers power.
p_values = np.asarray(group_result["p"].data)
n_voxels = p_values.size
fdr_thr = fdr(p_values, q=0.05)
bonf_thr = 0.05 / n_voxels
n_uncorrected = int((p_values < 0.001).sum())
n_fdr = int((p_values <= fdr_thr).sum()) if fdr_thr > 0 else 0
n_bonferroni = int((p_values < bonf_thr).sum())
print(f"voxels surviving, out of {n_voxels}:")
print(f" uncorrected (p < 0.001): {n_uncorrected:5d}")
print(f" FDR (q = 0.05): {n_fdr:5d}")
print(f" Bonferroni (p < 0.05/N): {n_bonferroni:5d}")voxels surviving, out of 71020:
uncorrected (p < 0.001): 127
FDR (q = 0.05): 0
Bonferroni (p < 0.05/N): 0
Recap¶
| Stage | What it does | Key API |
|---|---|---|
| Build design | BIDS events → HRF-convolved regressors + confounds + drift | DesignMatrix(events, run_length=, TR=), .append(confounds, axis=1, as_confounds=True), .add_poly() |
| First level | OLS at every voxel | brain.fit(X=design) |
| Contrast | Linear combination of βs (effect size + inference) | brain.compute_contrasts("A - B", statistic="all") |
| Stack subjects | Concatenate first-level β maps | concatenate([...]) |
| Group test | Voxelwise one-sample t-test → {mean, t, z, p} | group.ttest() |
| Correction | FDR threshold | nltools.algorithms.fdr, nltools.algorithms.threshold |
The per-subject loop is the explicit path; BrainCollection will wrap multi-subject fitting into a single call once it lands on this branch.
Next steps
Encoding models — predict brain activity from stimulus features (GLM vs. Ridge).
Multivariate pattern analysis — decode conditions and compare representational geometry.
Inter-subject correlation — shared responses to naturalistic stimuli.