Building fMRI Regressors

Introduction: What is a Regressor?

In fMRI analysis, a regressor (or predictor) represents the expected BOLD signal timecourse associated with a specific experimental condition or event type. It’s typically created by convolving a series of event onsets (often represented as delta functions or “sticks”) with a hemodynamic response function (HRF).

fmrihrf provides the regressor() function to easily create these objects from event timings and an HRF. While these regressor objects are often constructed automatically by modeling functions in other packages, this vignette explores how to create and manipulate them directly, offering finer control over the model components.

Basic Regressor from Event Onsets

Suppose we have a simple event-related fMRI design with stimuli presented every 12 seconds. We want to model these events using the SPM canonical HRF (HRF_SPMG1). The events are brief, so we model them with a duration of 0 seconds (instantaneous).

# Define event onsets
onsets <- seq(0, 10 * 12, by = 12)

# Create the regressor object
# Uses HRF_SPMG1 by default if no hrf is specified
# Duration is 0 by default
reg1 <- regressor(onsets = onsets, hrf = HRF_SPMG1)

# Access components using helper functions
head(onsets(reg1))
#> [1]  0 12 24 36 48 60
nbasis(reg1)
#> [1] 1

Evaluating and Plotting a Regressor

A regressor object stores the event information but doesn’t automatically compute the timecourse. To get the predicted BOLD signal at specific time points (e.g., corresponding to scan acquisition times), we use the evaluate() function.

# Define a time grid corresponding to scan times (e.g., TR=2s)
TR <- 2
scan_times <- seq(0, 140, by = TR)

# Plot the regressor using the plot() method
# This automatically evaluates and plots with onset markers
plot(reg1, grid = scan_times)

Modeled SPM response to eleven events spaced 12 seconds apart. Dashed vertical lines mark event onsets; the response is sampled every 2 seconds.

Dashed lines mark event onsets; the curve is the modeled response.

Varying Event Durations

Sometimes events have different durations. The duration argument in regressor() can take a vector matching the length of onsets.

# Example onsets and durations
onsets_var_dur <- seq(0, 5 * 12, length.out = 6)
durations_var <- 1:length(onsets_var_dur) # Durations increase from 1s to 6s

# Create regressor with varying durations
reg_var_dur <- regressor(onsets_var_dur, HRF_SPMG1, duration = durations_var)

# Plot the regressor
scan_times_dur <- seq(0, max(onsets_var_dur) + 30, by = TR)
plot(reg_var_dur, grid = scan_times_dur)

Duration and Summation

By default (summate=TRUE), the predicted response accumulates if events overlap or have extended duration. Setting summate=FALSE preserves the same temporal profile but the peak amplitude does not grow with duration.

# Create regressor with varying durations, summate=FALSE
reg_var_dur_nosum <- regressor(onsets_var_dur, HRF_SPMG1,
                               duration = durations_var, summate = FALSE)

# Compare summating vs non-summating using plot_regressors()
plot_regressors(reg_var_dur, reg_var_dur_nosum,
                labels = c("summate=TRUE", "summate=FALSE"),
                grid = scan_times_dur,
                title = "Effect of Summation on Response",
                subtitle = "Same events with varying durations")

Varying Event Amplitudes (Parametric Modulation)

We can model variations in event intensity or some associated parameter by providing an amplitude vector. This creates a parametric regressor where the height of the HRF for each event is scaled by the corresponding amplitude value.

# Example onsets and amplitudes (e.g., representing task difficulty)
onsets_amp <- seq(0, 10 * 12, length.out = 11)
amplitudes_raw <- 1:length(onsets_amp)

# It's common practice to center the modulator
amplitudes_scaled <- scale(amplitudes_raw, center = TRUE, scale = FALSE)

# Create the parametric regressor
reg_amp <- regressor(onsets_amp, HRF_SPMG1, amplitude = amplitudes_scaled)

# Plot the parametric regressor
scan_times_amp <- seq(0, max(onsets_amp) + 30, by = TR)
plot(reg_amp, grid = scan_times_amp)

Continuous Features

A sampled feature such as RMS energy is a time series, not a list of trials. feature_regressor() encodes each sample as a zero-order-hold bin of width dt and convolves that signal with the HRF. That is the same linear model as amplitude modulation on this sampling grid: the predicted BOLD is $Hx$ .

By default the series is demeaned before convolution (center = TRUE) and left in native units. For a whole-run series that is $H(x-\mu 1)=Hx-\mu H1$ . In the interior of the run, overlapping HRFs make $H1$ nearly constant, so with a GLM intercept the centered and raw columns test the same effect. They differ by the HRF-length ramp of $H1$ at the run boundaries; centering removes that onset/offset transient. scale = "sd" z-scores the feature, not the final filtered design column.

dt <- 0.1
feat_times <- seq(0, 20, by = dt)
# Simulated acoustic envelope
rms <- abs(sin(2 * pi * feat_times / 8)) * (0.5 + 0.5 * sin(2 * pi * feat_times / 20))

feat <- feature_regressor(rms, dt = dt, hrf = HRF_SPMG1)
scan_times_feat <- seq(0, max(feat_times) + 30, by = TR)
plot(feat, grid = scan_times_feat, precision = dt)

Using regressor(times, amplitude = rms, duration = 0) instead would treat each sample as a unit-mass impulse and scale the predicted BOLD by about 1/dt. Pass duration = dt (and skip centering) if you need the same ZOH encoding from regressor().

If you want intensity conditional on an on-period, pass a mask for those samples (off-period stays 0 after centering) and a separate boxcar for presence. That is not an affine transform of the all-sample series. Do not drop off-period samples before centering: that silently turns the all-sample model into the sparse one.

Combining Duration and Amplitude Modulation

You can provide both duration and amplitude vectors to model events that vary in both aspects.

set.seed(123)
onsets_comb <- seq(0, 10 * 12, length.out = 11)
amps_comb <- scale(1:length(onsets_comb), center = TRUE, scale = FALSE)
durs_comb <- sample(1:5, length(onsets_comb), replace = TRUE)

reg_comb <- regressor(onsets_comb, HRF_SPMG1, 
                      amplitude = amps_comb, duration = durs_comb)

# Plot the combined regressor
scan_times_comb <- seq(0, max(onsets_comb) + 30, by = TR)
plot(reg_comb, grid = scan_times_comb)

Regressors with HRF Basis Sets

If you use an HRF object with multiple basis functions (e.g., HRF_SPMG3, HRF_BSPLINE), the regressor object will represent multiple timecourses, one for each basis function. evaluate() will return a matrix.

# Use a B-spline basis set
onsets_basis <- seq(0, 10 * 12, length.out = 11)
hrf_basis <- HRF_BSPLINE # Uses N=5 basis functions by default

reg_basis <- regressor(onsets_basis, hrf_basis)
nbasis(reg_basis) # Should be 5
#> [1] 5

# Evaluate - this returns a matrix
scan_times_basis <- seq(0, max(onsets_basis) + 30, by = TR)
pred_basis_matrix <- evaluate(reg_basis, scan_times_basis)
dim(pred_basis_matrix) # rows = time points, cols = basis functions
#> [1] 76  5

# Plot using the plot() method - automatically handles multi-basis
plot(reg_basis, grid = scan_times_basis)

Shifting Regressors

You can temporally shift all onsets within a regressor using the shift() method.

# Original regressor
reg_orig <- regressor(onsets = c(10, 30, 50), hrf = HRF_SPMG1)

# Shifted regressor (delay by 5 seconds)
reg_shifted <- shift(reg_orig, shift_amount = 5)

onsets(reg_orig)
#> [1] 10 30 50
onsets(reg_shifted) # Onsets are now 15, 35, 55
#> [1] 15 35 55

# Compare original and shifted using plot_regressors()
scan_times_shift <- seq(0, 80, by = TR)
plot_regressors(reg_orig, reg_shifted,
                labels = c("Original", "Shifted +5s"),
                grid = scan_times_shift,
                show_onsets = TRUE,  # Show onsets for both
                title = "Shifting a Regressor")

Bradley R. Buchsbaum

2026-09-28