Nuisance Regression in the fMRI GLM

Overview

This assignment extends a completed first-level FSL FEAT analysis by rerunning the same task model in three versions. The baseline model includes task regressors only. The second model adds motion parameters and motion derivatives. The third model adds motion parameters, motion derivatives, and five aCompCor components.

Throughout the tutorial, the example assumes three FEAT directories:

M1.feat  # Baseline task model
M2.feat  # Task model + motion
M3.feat  # Task model + motion + aCompCor

If your actual output folders have different names, replace the directory names in the code blocks.

Part 1. Three model variants

The task EVs should be identical across the three models. The only difference should be the nuisance regressors. This makes the comparison interpretable.

Prepare the confound matrices

You need two nuisance-regressor text files. One should contain the 12 motion columns. The other should contain those same 12 motion columns plus five aCompCor columns.

Expected columns

Generic Python template

Use this if your confounds are in a tab-separated file. Adjust the column names to match your actual confounds.tsv.

import pandas as pd

conf = pd.read_csv("confounds.tsv", sep="\t")

motion_cols = [
    "trans_x", "trans_y", "trans_z",
    "rot_x", "rot_y", "rot_z",
]

motion_derivative_cols = [
    "trans_x_derivative1", "trans_y_derivative1", "trans_z_derivative1",
    "rot_x_derivative1", "rot_y_derivative1", "rot_z_derivative1",
]

acompcor_cols = [
    "a_comp_cor_00", "a_comp_cor_01", "a_comp_cor_02",
    "a_comp_cor_03", "a_comp_cor_04",
]

m2_cols = motion_cols + motion_derivative_cols
m3_cols = motion_cols + motion_derivative_cols + acompcor_cols

conf[m2_cols].fillna(0).to_csv(
    "motion_12.txt",
    sep=" ",
    header=False,
    index=False
)

conf[m3_cols].fillna(0).to_csv(
    "motion_acompcor_17.txt",
    sep=" ",
    header=False,
    index=False
)

Check row count

wc -l motion_12.txt
wc -l motion_acompcor_17.txt
fslnvols M1.feat/filtered_func_data.nii.gz

Add confounds in FEAT

1. Create M2 from M1

   Open the FEAT design used for M1. Save a copy as the M2 design, then add motion_12.txt as the confound EVs text file if your FSL version supports it.

2. Create M3 from M1 or M2

   Save another copy of the design and add motion_acompcor_17.txt as the confound EVs text file.

3. If using the EV interface manually

   Add each nuisance column as a separate EV. Set the waveform to Custom (1 entry per volume), set convolution to None, and leave Add temporal derivative unchecked.

4. Run all three models

   feat M1.fsf
   feat M2.fsf
   feat M3.fsf

   You may also run them from the FEAT GUI. Keep careful notes on the output directories.

Write-up question 1

10 pts Report the total number of regressors in each model. Explain why adding regressors reduces residual degrees of freedom and when adding too many nuisance regressors becomes counterproductive.

Answer key

If FEAT includes additional columns for filtering, temporal derivatives, or other automatically modeled terms, your design matrix may show more columns. For the assignment logic, the expected counts are 3, 15, and 20 when counting task regressors, nuisance regressors, and intercept.

Suggested written response

Adding nuisance regressors increases the number of columns in the design matrix. Each additional estimated beta consumes one degree of freedom, so the residual degrees of freedom decrease as the model becomes larger. This can be beneficial when the nuisance regressors explain non-neural variance, because residual variance decreases and the task estimates become more precise. However, adding too many regressors can become counterproductive if they remove meaningful task-related variance, introduce collinearity, or overfit idiosyncratic noise. The best confound model should improve fit while preserving enough degrees of freedom to estimate the task effects reliably.

Write-up question 2

10 pts Compute whole-brain mean R² for M1, M2, and M3. Then compute R² by tissue class using masks.

Whole-brain R²

fslstats M1.feat/stats/r2.nii.gz -M
fslstats M2.feat/stats/r2.nii.gz -M
fslstats M3.feat/stats/r2.nii.gz -M

R² by tissue class

for MODEL in M1 M2 M3; do
  echo "----- ${MODEL} -----"
  echo "Whole brain:"
  fslstats ${MODEL}.feat/stats/r2.nii.gz -M

  echo "Grey matter:"
  fslstats ${MODEL}.feat/stats/r2.nii.gz -k gm_mask.nii.gz -M

  echo "White matter:"
  fslstats ${MODEL}.feat/stats/r2.nii.gz -k wm_mask.nii.gz -M

  echo "CSF:"
  fslstats ${MODEL}.feat/stats/r2.nii.gz -k csf_mask.nii.gz -M
done

Answer key interpretation

R² should usually increase from M1 to M2 and from M2 to M3, because each added set of nuisance regressors explains additional variance in the BOLD signal. The increase is often especially visible in white matter and CSF, because aCompCor components are derived from those tissues and are designed to capture non-neural physiological variance. Grey matter may also show increased R² if physiological noise is shared across tissue compartments. A larger R² is not automatically better, because R² always tends to increase when predictors are added. The key question is whether the added regressors reduce non-neural residual variance without distorting the task contrast.

Part 2. Beta weight comparison

FEAT stores contrast estimates in stats/cope*.nii.gz and variance estimates in stats/varcope*.nii.gz. This section compares the estimated task effect across the three model variants.

Write-up question 3

10 pts Create a difference map between M1 and M3 for the primary contrast of interest.

Code

fslmaths M3.feat/stats/cope1.nii.gz \
  -sub M1.feat/stats/cope1.nii.gz \
  beta_diff_M3_minus_M1.nii.gz

fsleyes M1.feat/mean_func.nii.gz beta_diff_M3_minus_M1.nii.gz &

What to look for

• Are the largest differences near tissue boundaries?
• Are differences concentrated near ventricles, large vessels, or edges of the brain?
• Do differences appear in homogeneous grey matter regions?
• Are M3 estimates generally higher, lower, or spatially mixed relative to M1?

Answer key interpretation

Differences near tissue boundaries, ventricles, and edge regions often suggest that M3 is removing variance related to motion or physiological artifacts. If differences occur in task-relevant grey matter regions, the interpretation depends on direction. M3 estimates may decrease if M1 had absorbed motion or physiological variance into the task contrast. M3 estimates may increase if removing nuisance variance reduces residual noise and reveals a cleaner task effect.

Write-up question 4

10 pts Extract the cope value from a 5 mm sphere around the peak activation coordinate and compute Cohen's d for each model.

Step 1. Find the peak activation voxel in M1

fslstats M1.feat/stats/zstat1.nii.gz -x

Record the three voxel coordinates returned by FSL. The coordinates are zero-indexed voxel coordinates.

Optional cluster table

cluster \
  --in=M1.feat/stats/zstat1.nii.gz \
  --thresh=3.1 \
  --mm \
  > M1_clusters.txt

cat M1_clusters.txt

Step 2. Create a single-voxel seed image

Replace X, Y, and Z with your peak coordinates.

X=45
Y=37
Z=28

fslmaths M1.feat/stats/zstat1.nii.gz \
  -mul 0 -add 1 \
  -roi ${X} 1 ${Y} 1 ${Z} 1 0 1 \
  peak_seed.nii.gz

fsleyes M1.feat/stats/zstat1.nii.gz peak_seed.nii.gz &

Step 3. Expand the seed to a 5 mm sphere

fslmaths peak_seed.nii.gz \
  -kernel sphere 5 \
  -fmean \
  -bin \
  peak_sphere_5mm.nii.gz

fslstats peak_sphere_5mm.nii.gz -V

Step 4. Extract cope and varcope values

COPE_M1=$(fslstats M1.feat/stats/cope1.nii.gz -k peak_sphere_5mm.nii.gz -m)
COPE_M2=$(fslstats M2.feat/stats/cope1.nii.gz -k peak_sphere_5mm.nii.gz -m)
COPE_M3=$(fslstats M3.feat/stats/cope1.nii.gz -k peak_sphere_5mm.nii.gz -m)

VARCOPE_M1=$(fslstats M1.feat/stats/varcope1.nii.gz -k peak_sphere_5mm.nii.gz -m)
VARCOPE_M2=$(fslstats M2.feat/stats/varcope1.nii.gz -k peak_sphere_5mm.nii.gz -m)
VARCOPE_M3=$(fslstats M3.feat/stats/varcope1.nii.gz -k peak_sphere_5mm.nii.gz -m)

echo "M1 cope=${COPE_M1} varcope=${VARCOPE_M1}"
echo "M2 cope=${COPE_M2} varcope=${VARCOPE_M2}"
echo "M3 cope=${COPE_M3} varcope=${VARCOPE_M3}"

Step 5. Compute Cohen's d

python3 - <8} {'varcope':>10} {'Cohen d':>9}")
print("-" * 62)

for name, cope, varcope in models:
    d = cope / math.sqrt(varcope)
    print(f"{name:<28} {cope:>8.4f} {varcope:>10.6f} {d:>9.4f}")
EOF

Answer key interpretation

Cohen's d can increase if nuisance regression lowers varcope more than it changes the cope. It can decrease if the original task estimate was inflated by task-correlated noise, such as motion that increased during the more difficult condition. A decrease is not necessarily bad. It may mean the nuisance model produced a less biased estimate of the true task effect.

Bonus. Residual power spectral density

This optional section extracts residual time series from stats/res4d.nii.gz and plots power spectral density for a small grey matter ROI.

Create or choose a small grey matter ROI

fslmeants -i M1.feat/stats/res4d.nii.gz -m gm_roi_3x3x3.nii.gz -o M1_residual_roi.txt
fslmeants -i M2.feat/stats/res4d.nii.gz -m gm_roi_3x3x3.nii.gz -o M2_residual_roi.txt
fslmeants -i M3.feat/stats/res4d.nii.gz -m gm_roi_3x3x3.nii.gz -o M3_residual_roi.txt

Plot the spectra in Python

python3 - <<'EOF'
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import welch

TR = 2.0  # Replace with your actual TR in seconds.
fs = 1 / TR

files = {
    "M1 baseline": "M1_residual_roi.txt",
    "M2 + motion": "M2_residual_roi.txt",
    "M3 + motion + aCompCor": "M3_residual_roi.txt",
}

plt.figure(figsize=(8, 5))

for label, path in files.items():
    y = np.loadtxt(path)
    y = np.asarray(y).squeeze()
    freqs, psd = welch(y, fs=fs, nperseg=min(128, len(y)))
    plt.plot(freqs, psd, label=label)

plt.xlabel("Frequency (Hz)")
plt.ylabel("Power spectral density")
plt.title("Residual power spectrum by nuisance model")
plt.legend()
plt.tight_layout()
plt.savefig("residual_psd_M1_M2_M3.png", dpi=300)
EOF

Answer key interpretation

Motion regressors often reduce transient or broadband residual power. aCompCor can reduce residual power at frequencies associated with respiration, cardiac pulsatility, and spatially shared physiological noise. The exact pattern depends on TR, sampling rate, preprocessing, and the physiological rhythms present during the scan.

Part 3. Reflection methods paragraph

15 pts Write a brief methods paragraph, approximately 150 words, describing the M3 confound strategy in third person and past tense.

Model answer

First-level analyses were conducted in FSL FEAT using a general linear model that included task regressors for the experimental conditions and nuisance regressors for non-neural sources of variance. Functional data were motion corrected using rigid-body volume registration, which estimated three translation and three rotation parameters for each volume. The final nuisance model included the six motion parameters and their six temporal derivatives, yielding 12 motion-related regressors. Slow scanner drift was addressed using FEAT's high-pass filtering procedure with a 1/128 Hz cutoff. Physiological noise was modeled using an anatomical CompCor approach. White matter and cerebrospinal fluid masks were used to identify tissue compartments unlikely to contain task-related neural signal, and principal components were extracted from those regions. The first five aCompCor components were included as additional nuisance regressors. All nuisance regressors were entered without HRF convolution, and task effects were estimated after accounting for variance shared with the confound regressors.