18. The General Linear Model and Foundations of Analysis

The General Linear Model and Foundations of Analysis¶

Overview¶

The General Linear Model is arguably the most widely used statistical method for relating task manipulations to brain activity. Applied to fMRI, it identifies regions that respond to an event type, compares responses across conditions, and links brain activity to behavioral or clinical variables. In a typical analysis the outcome variable is the time series from a single voxel (or region of interest), and the analysis is massively univariate: a separate GLM is fit at every voxel, and the resulting parameter estimates are assembled into statistical maps.

The GLM is really a family of models in which one continuously distributed outcome ($y$) is predicted by a linear combination of predictors (the columns of a design matrix $X$). This family includes t-tests, multiple regression, ANOVA, and ANCOVA — the number, distribution, and grouping of predictors determines which label applies. GLMs have three big advantages: they have a closed-form solution so estimation is fast and reproducible; they are easy to interpret; and they are flexible, since nonlinear relationships can be captured with transformed predictors, piecewise terms, or splines. Extensions handle situations where the basic assumptions fail: generalized least squares for correlated or unequal-variance errors, linear mixed models for grouped observations such as participants, and generalized linear models (via link functions) for binary or count outcomes.

Modeling task-evoked BOLD signal begins with two ingredients: a neural stimulus function describing when neural activity is assumed to occur — brief events or sustained epochs (blocks) — and a hemodynamic response function describing the sluggish vascular response to a burst of neural activity. Assuming the brain behaves as a linear time invariant (LTI) system, the predicted signal is the convolution of the two:

LTI systems have three key properties: scaling (doubling neural activity doubles the BOLD response, so amplitude differences between conditions can be interpreted as neural differences), superposition (responses to nearby events sum), and time invariance (shifting a stimulus shifts its response). Linearity is a good approximation when events are spaced at least ~5 seconds apart, though nonlinearities (e.g., refractory effects) can be substantial for stimuli spaced under ~2 seconds.

A popular canonical HRF is the difference of two gamma functions, which captures the initial peak around 5–6 seconds and the subsequent undershoot:

with common choices $\alpha_1 = 6$, $\alpha_2 = 16$, $\beta_1 = \beta_2 = 1$, and $c = 1/6$.

To formulate the GLM, each condition’s indicator vector (1 during hypothesized neural activity, 0 elsewhere) is convolved with the HRF, and the resulting predicted time courses become columns of the design matrix $X$, alongside an intercept and any nuisance covariates (head motion estimates, drift terms). With $n$ time points and $p$ predictors, the model is

Estimation seeks the $\beta$ values that minimize the sum of squared errors $(y - X\beta)^T(y - X\beta)$. A closed-form solution exists — the ordinary least squares (OLS) estimate:

with error variance estimated from the residuals $r = y - X\hat{\beta}$ as $\hat{\sigma}^2 = \frac{r^Tr}{n - p}$. Under the model assumptions, $\hat{\beta}$ is the best linear unbiased estimate (BLUE), distributed as $\hat{\beta} \sim N(\beta,\; \sigma^2 (X^TX)^{-1})$.

Inference on each $\hat{\beta}_i$ compares it with a null value of zero. Its standard error is the square root of the $i$th diagonal of $\hat{\sigma}^2 (X^TX)^{-1}$, and the ratio $t = \hat{\beta}_i / SE(\hat{\beta}_i)$ follows a Student’s t distribution with $df_e = n - p$ error degrees of freedom. A low P value lets us declare a region activated ($\hat{\beta}_i > 0$) or deactivated ($\hat{\beta}_i < 0$) — always relative to the model’s intercept, which reflects unmodeled baseline periods. The same machinery extends to contrasts across conditions (Chapter 20).

Two refinements matter in practice. First, fMRI errors are not independent — they are autocorrelated — so OLS standard errors are too small and t-statistics inflated. The remedy is generalized least squares (GLS): assume a structured error covariance $V$ and prewhiten, giving $\hat{\beta} = (X^TV^{-1}X)^{-1}X^TV^{-1}y$, with $V$ and $\beta$ estimated iteratively (Chapter 19). Second, the canonical HRF’s shape varies across brain regions and individuals, and mis-modeling it reduces power and can produce false positives. Basis sets — such as the canonical HRF plus temporal and dispersion derivatives, or the highly flexible finite impulse response (FIR) model with one parameter per post-stimulus time point — let the fitted response adapt, at the cost of a fundamental tradeoff between flexibility and power: too much flexibility risks fitting noise.

Hands-on tutorial¶

In this tutorial you will build a design matrix from event onsets and fit a GLM to a simulated voxel time series — the core of every first-level fMRI analysis. The figure below shows the recipe: indicator functions marking event onsets are convolved with an assumed HRF, and the resulting predictors become columns of $X$.

Step 1 — Build a design matrix from event onsets. We specify onset times (in seconds) for two event types, A and B, and convolve them with a canonical HRF.

% Requires CanlabCore + SPM12 on your MATLAB path
% Adapted from CANlab tutorials (github.com/canlab)
TR = 2;                                       % repetition time (s)
n_scans = 180;                                % number of volumes (6-min run)

ons = {};
ons{1} = [10 60 110 160 210 260 310]';        % Condition A onsets (sec)
ons{2} = [35 85 135 185 235 285 335]';        % Condition B onsets (sec)

% Convolve onsets with canonical HRF; intercept is the last column
X = onsets2fmridesign(ons, TR, n_scans * TR);

plotDesign(ons, [], TR);                       % plot regressors and onsets
figure; imagesc(X); colormap gray;             % view X as an image
xlabel('Regressor'); ylabel('Time (TRs)');

import numpy as np, pandas as pd
from nilearn.glm.first_level import make_first_level_design_matrix
from nilearn.plotting import plot_design_matrix

t_r, n_scans = 2.0, 180        # t_r = repetition time (s); n_scans = volumes (6-min run)
frame_times = np.arange(n_scans) * t_r

events = pd.DataFrame({
    "onset":      [10, 60, 110, 160, 210, 260, 310,     # Condition A
                   35, 85, 135, 185, 235, 285, 335],    # Condition B
    "duration":   [1.0] * 14,                          # every event lasts 1 s
    "trial_type": ["A"] * 7 + ["B"] * 7,
})

# Convolve onsets with canonical HRF; 'constant' (intercept) is last column
X = make_first_level_design_matrix(frame_times, events,
                                   hrf_model="spm", drift_model=None)
plot_design_matrix(X)                          # view X as an image

Example output: running the Python tab draws the classic design-matrix image — 180 rows (one per volume) by 3 columns (regressors A and B, plus the constant/intercept):

Step 2 — Fit the GLM to a voxel time series. We simulate a voxel with known true effects ($\beta_A = 0.8$, $\beta_B = 0.4$, intercept 100) plus modest Gaussian noise ($\sigma = 0.5$, where $\sigma$ is the noise standard deviation) and recover them with the OLS formula $\hat{\beta} = (X^TX)^{-1}X^Ty$.

% Simulate a voxel: y = X * beta_true + noise
rng(9);                                        % seed for reproducible noise
beta_true = [0.8 0.4 100]';                    % true effects: A, B, intercept (last col)
sigma_noise = 0.5;                             % noise SD; modest so this short 7-event demo recovers betas clearly
y = X * beta_true + sigma_noise * randn(n_scans, 1);

% Ordinary least squares: beta_hat = (X'X)^{-1} X'y
beta_hat = (X' * X) \ (X' * y);

r      = y - X * beta_hat;                     % residuals
dfe    = n_scans - size(X, 2);                 % error degrees of freedom (n - p)
sigma2 = (r' * r) / dfe;                       % error variance estimate
t_vals = beta_hat ./ sqrt(sigma2 * diag(inv(X' * X)));

disp(table(beta_true, beta_hat, t_vals))

# Simulate a voxel: y = X @ beta_true + noise
Xm = X.to_numpy().copy()                       # design as array; columns: A, B, constant
Xm[:, :2] /= Xm[:, :2].max(axis=0)             # rescale task regressors to unit peak, so betas are in signal units
rng = np.random.default_rng(1)                 # seed for reproducible noise
beta_true = np.array([0.8, 0.4, 100.0])        # true effects: A, B, intercept
sigma_noise = 0.5                              # noise SD; modest so this short 7-event demo recovers betas clearly
y = Xm @ beta_true + sigma_noise * rng.standard_normal(n_scans)

# Ordinary least squares: beta_hat = (X'X)^{-1} X'y
beta_hat = np.linalg.solve(Xm.T @ Xm, Xm.T @ y)

resid  = y - Xm @ beta_hat                     # residuals
dfe    = n_scans - Xm.shape[1]                 # error degrees of freedom (n - p)
sigma2 = resid @ resid / dfe                   # error variance estimate
t_vals = beta_hat / np.sqrt(sigma2 * np.diag(np.linalg.inv(Xm.T @ Xm)))

print(pd.DataFrame({"true": beta_true, "estimate": beta_hat, "t": t_vals},
                   index=X.columns))

Example output: the Python tab prints the recovered parameters — close to the true values, with clear t-statistics for both conditions (the MATLAB tab produces an analogous table):

           true   estimate            t
A           0.8   0.787951     6.314865
B           0.4   0.399651     3.223820
constant  100.0  99.964533  2647.254402

The estimates land close to the true values, with reliable t-statistics for both conditions. The full labs continue the arc: inspecting residuals and model fit ($R^2$), computing standard errors and P values, and testing a simple A − B contrast — the bridge to Chapter 20.

Thought questions¶

1. Suppose your experiment presents stimuli in rapid sequences, spaced 1–2 seconds apart. Which LTI assumptions are most at risk, how might your $\hat{\beta}$ estimates be distorted, and what design or modeling choices could mitigate the problem?

2. The HRF varies across brain regions and individuals — for instance, patients with vascular disease may show delayed or blunted responses. If you compare patients and controls using a fixed canonical HRF, how could HRF mis-modeling masquerade as group differences in “activation”?

3. Adding nuisance regressors (motion, drift) can both improve and worsen your inferences. Using the concepts of residual variance and error degrees of freedom, explain when each outcome would occur.

4. Activation in the GLM is always relative to the intercept (baseline). How does your choice of what to leave unmodeled change the interpretation of a significant positive beta? Give a concrete example where two labs modeling the same data differently would reach different verbal conclusions.

5. The FIR basis set can fit almost any response shape, yet most researchers use a single canonical HRF. Drawing on the power–flexibility tradeoff, when would you recommend each approach, and what would you look at in the data to decide?

Quiz yourself¶

⌂ Back to home · Table of contents · How to use this site