Hemodynamic Response Functions

Bradley R. Buchsbaum

2026-10-03

Introduction to Hemodynamic Response Functions (HRFs)

A hemodynamic response function (HRF) models the temporal evolution of the fMRI BOLD (Blood-Oxygen-Level-Dependent) signal in response to a brief neural event. Typically, the BOLD signal peaks 4-6 seconds after the event onset and then returns to baseline, often with a slight undershoot.

fmrihrf provides tools to define, manipulate, and visualize various HRFs commonly used in fMRI analysis.

Pre-defined HRF Objects

fmrihrf includes several pre-defined HRF objects, which are essentially functions with specific attributes defining their type, number of basis functions (nbasis), and effective duration (span).

Let’s look at two common examples: the SPM canonical HRF (HRF_SPMG1) and a Gaussian HRF (HRF_GAUSSIAN).

# SPM canonical HRF (based on difference of two gamma functions)
print(HRF_SPMG1)
#> -- HRF: SPMG1 --------------------------------------------- 
#>    Basis functions: 1 
#>    Span: 24 s
#>    Parameters: P1 = 5, P2 = 15, A1 = 0.008333

# Gaussian HRF
print(HRF_GAUSSIAN)
#> -- HRF: gaussian ------------------------------------------ 
#>    Basis functions: 1 
#>    Span: 24 s
#>    Parameters: mean = 6, sd = 2

These objects are functions themselves, so you can evaluate them at specific time points. The plot_hrfs() function draws one or more HRFs on shared axes; every figure in these vignettes uses its colours (see hrf_palette()):

time_points <- seq(0, 25, by = 0.1)

# normalize = TRUE scales each HRF to peak at 1.0
plot_hrfs(HRF_SPMG1, HRF_GAUSSIAN,
          labels = c("SPM canonical", "Gaussian"),
          normalize = TRUE,
          time = time_points,
          title = "SPM and Gaussian HRFs",
          subtitle = "Each curve normalized to peak at 1")

SPM and Gaussian HRFs normalized to a peak of one. The SPM response includes a negative undershoot after about 12 seconds; the Gaussian does not.

Note that the span attribute (e.g., 24 seconds) indicates the approximate time window over which the HRF is non-zero.

For a quick look at a single HRF, the base-graphics plot() method draws it and marks the peak: plot(HRF_SPMG1).

Choosing a Fixed HRF Scale

HRF scaling changes the units of fitted coefficients, so use an explicit convention when comparing designs across software. hrf_norm = "spm" matches the Nilearn/SPM reference-grid convention; "unit_peak" sets the canonical peak to one, and "unit_integral" sets its continuous integral to one.

spm_scaled <- gen_hrf(HRF_SPMG1, hrf_norm = "spm")
spm_grid <- seq(0, 32, length.out = 1600)
sum(spm_scaled(spm_grid))
#> [1] 1

Use normalize_hrf() to scale an existing HRF object. For SPM derivative bases, the modes above apply one canonical-derived factor to every column and therefore preserve their relative scale. The legacy normalise_hrf() function instead gives every basis column its own unit peak.

Modifying HRF Parameters with gen_hrf

The gen_hrf function is a flexible way to create new HRF functions, often by modifying the parameters of existing ones.

For example, the hrf_gaussian function takes mean and sd arguments. We can use gen_hrf to create Gaussian HRFs with different peak times (mean) and widths (sd).

# Create Gaussian HRFs with different parameters using gen_hrf
# Note: hrf_gaussian is the underlying function, not the HRF object HRF_GAUSSIAN
hrf_gauss_4_1 <- gen_hrf(hrf_gaussian, mean = 4, sd = 1, name = "Gaussian (Mean=4, SD=1)")
hrf_gauss_5_2 <- gen_hrf(hrf_gaussian, mean = 5, sd = 2, name = "Gaussian (Mean=5, SD=2)")
hrf_gauss_7_3 <- gen_hrf(hrf_gaussian, mean = 7, sd = 3, name = "Gaussian (Mean=7, SD=3)")
plot_hrfs(hrf_gauss_4_1, hrf_gauss_5_2, hrf_gauss_7_3,
          labels = c("mean 4, sd 1", "mean 5, sd 2", "mean 7, sd 3"),
          time = time_points, palette = "ordered",
          title = "Gaussian HRF parameters",
          subtitle = "Later mean, later peak; larger sd, broader")

Three Gaussian HRFs. Later means peak later, and larger standard deviations give lower, broader curves.

gen_hrf can also directly incorporate lags and durations (see later sections).

Modeling Event Duration with block_hrf

fMRI events often have a duration (e.g., a stimulus presented for several seconds). The block_hrf function (or gen_hrf with a width argument) modifies an HRF to model the response to a sustained event of a specific width (duration). Internally, it convolves the original HRF with a boxcar function of the specified width.

The precision argument controls the sampling resolution used for this convolution.

# Create blocked HRFs using the SPM canonical HRF with different durations
hrf_spm_w1 <- block_hrf(HRF_SPMG1, width = 1)
hrf_spm_w2 <- block_hrf(HRF_SPMG1, width = 2)
hrf_spm_w4 <- block_hrf(HRF_SPMG1, width = 4)

The three durations below and in the next sections use the ordered palette: violet is the shortest event and ochre the longest.

plot_hrfs(hrf_spm_w1, hrf_spm_w2, hrf_spm_w4,
          labels = c("1 s event", "2 s event", "4 s event"),
          time = time_points, palette = "ordered",
          title = "Longer events: larger responses",
          subtitle = "block_hrf(HRF_SPMG1, width = 1, 2, 4)")

SPM canonical responses to 1, 2 and 4 second events. Longer events produce higher and later peaks.

Normalization

By default, longer durations lead to higher peak responses (assuming summation, see next section). Setting normalize=TRUE in block_hrf (or gen_hrf) rescales the response so the peak amplitude is 1, regardless of duration. What remains is the change in shape: longer events rise later and stay up longer.

# Create normalized blocked HRFs
hrf_spm_w1_norm <- block_hrf(HRF_SPMG1, width = 1, normalize = TRUE)
hrf_spm_w2_norm <- block_hrf(HRF_SPMG1, width = 2, normalize = TRUE)
hrf_spm_w4_norm <- block_hrf(HRF_SPMG1, width = 4, normalize = TRUE)
plot_hrfs(hrf_spm_w1_norm, hrf_spm_w2_norm, hrf_spm_w4_norm,
          labels = c("1 s event", "2 s event", "4 s event"),
          time = time_points, palette = "ordered",
          title = "Normalized: only shape changes",
          subtitle = "block_hrf(..., normalize = TRUE)")

Normalized SPM responses to 1, 2 and 4 second events. All peak at 1; longer events peak later and are broader.

Averaging instead of summing: summate

The summate argument in block_hrf controls how the response to a sustained event is built. With summate = TRUE (the default) the responses to each moment of the event add up, so longer events give larger peaks. With summate = FALSE they are averaged instead: the result is the summated response divided by the event width. The peak no longer grows with duration; for long events it falls, because the same response is spread over a longer time.

# Create non-summating blocked HRFs
hrf_spm_w2_nosum <- block_hrf(HRF_SPMG1, width = 2, summate = FALSE)
hrf_spm_w4_nosum <- block_hrf(HRF_SPMG1, width = 4, summate = FALSE)
hrf_spm_w8_nosum <- block_hrf(HRF_SPMG1, width = 8, summate = FALSE)
plot_hrfs(hrf_spm_w2_nosum, hrf_spm_w4_nosum, hrf_spm_w8_nosum,
          labels = c("2 s event", "4 s event", "8 s event"),
          time = time_points, palette = "ordered",
          title = "Averaging: peaks do not grow",
          subtitle = "block_hrf(..., summate = FALSE)")

Averaged (summate = FALSE) SPM responses to 2, 4 and 8 second events. Peaks do not grow with duration; the 8 second response is lower and broader.

summate = FALSE and normalize = TRUE can be combined; with a unit peak, only the shape differences shown in the normalized figure above remain.

Modeling Temporal Shifts with lag_hrf

Sometimes, the hemodynamic response might be delayed or advanced relative to the event onset. The lag_hrf function (or gen_hrf_lagged) shifts an existing HRF in time by a specified lag (in seconds). A positive lag delays the response, while a negative lag advances it.

# Create lagged versions of the Gaussian HRF
hrf_gauss_lag_neg2 <- lag_hrf(HRF_GAUSSIAN, lag = -2)
hrf_gauss_lag_0 <- HRF_GAUSSIAN # Original (lag=0)
hrf_gauss_lag_pos3 <- lag_hrf(HRF_GAUSSIAN, lag = 3)
plot_hrfs(hrf_gauss_lag_neg2, hrf_gauss_lag_0, hrf_gauss_lag_pos3,
          labels = c("lag -2 s", "lag 0 s", "lag +3 s"),
          time = time_points, palette = "ordered",
          title = "Lag moves the response in time",
          subtitle = "lag_hrf(HRF_GAUSSIAN, lag = -2, 0, 3)")

Gaussian HRF shifted by -2, 0 and +3 seconds. The shape is unchanged; only the peak time moves.

Combining Lag and Duration

We can combine lag_hrf and block_hrf using the pipe operator (%>%) from dplyr (or magrittr).

# Create HRFs that are both lagged and blocked
hrf_lb_1 <- HRF_GAUSSIAN %>% lag_hrf(1) %>% block_hrf(width = 1, normalize = TRUE)
hrf_lb_3 <- HRF_GAUSSIAN %>% lag_hrf(3) %>% block_hrf(width = 3, normalize = TRUE)
hrf_lb_5 <- HRF_GAUSSIAN %>% lag_hrf(5) %>% block_hrf(width = 5, normalize = TRUE)
plot_hrfs(hrf_lb_1, hrf_lb_3, hrf_lb_5,
          labels = c("lag 1 s, width 1 s", "lag 3 s, width 3 s", "lag 5 s, width 5 s"),
          time = time_points, palette = "ordered",
          title = "Lag and duration combined",
          subtitle = "lag_hrf() %>% block_hrf(normalize = TRUE)")

Gaussian HRFs with lag and width both equal to 1, 3 and 5 seconds, normalized to peak at 1. Larger values peak later and are broader.

Alternatively, gen_hrf can apply lag and width directly, and gives the same function:

hrf_lb_gen_3 <- gen_hrf(hrf_gaussian, lag = 3, width = 3, normalize = TRUE)

# Largest difference from the piped version over 0-25 s
max(abs(hrf_lb_gen_3(time_points) - hrf_lb_3(time_points)))
#> [1] 0.0004053635

Multivariate HRFs: Basis Sets

Instead of assuming a fixed HRF shape, we can model the response using a linear combination of multiple basis functions. This allows for more flexibility in capturing variations in HRF shape across brain regions or individuals. The resulting HRF function returns a matrix where each column corresponds to a basis function, and plot_hrfs() draws one curve per column.

SPM Basis Sets

fmrihrf provides pre-defined HRF objects for the SPM canonical HRF plus its temporal derivative (HRF_SPMG2), and additionally its dispersion derivative (HRF_SPMG3).

# SPM + Temporal Derivative (2 basis functions)
print(HRF_SPMG2)
#> -- HRF: SPMG2 --------------------------------------------- 
#>    Basis functions: 2 
#>    Span: 24 s

# SPM + Temporal + Dispersion Derivatives (3 basis functions)
print(HRF_SPMG3)
#> -- HRF: SPMG3 --------------------------------------------- 
#>    Basis functions: 3 
#>    Span: 24 s

HRF_SPMG2 is the first two columns of HRF_SPMG3, so one figure shows both:

plot_hrfs(HRF_SPMG3,
          labels = c("Canonical", "Temporal derivative", "Dispersion derivative"),
          time = time_points,
          title = "SPM canonical and derivatives")

Canonical SPM response and its temporal and dispersion derivatives, with all positive and negative values retained.

Why the temporal derivative? Adding a small multiple of it to the canonical HRF shifts the peak in time. This is how a GLM with HRF_SPMG2 absorbs small latency differences:

spm_early <- gen_empirical_hrf(time_points,
  drop(HRF_SPMG2(time_points) %*% c(1, 0.5)))
spm_late <- gen_empirical_hrf(time_points,
  drop(HRF_SPMG2(time_points) %*% c(1, -0.5)))

plot_hrfs(spm_early, spm_late,
          labels = c("canonical + 0.5 x derivative", "canonical - 0.5 x derivative"),
          time = time_points, reference = HRF_SPMG1, reference_label = "canonical",
          title = "The derivative shifts the peak")

Canonical SPM HRF (dashed grey) and the canonical plus or minus 0.5 times its temporal derivative. Adding the derivative moves the peak earlier; subtracting it moves the peak later.

B-Spline Basis Set

The hrf_bspline function generates a B-spline basis set. We typically use it within gen_hrf to create an HRF object. Key parameters are N (number of basis functions) and degree.

# B-spline basis with N=5 basis functions, degree=3 (cubic)
hrf_bs_5_3 <- gen_hrf(hrf_bspline, N = 5, degree = 3, name = "B-spline (N=5, deg=3)")
print(hrf_bs_5_3)
#> -- HRF: B-spline (N=5, deg=3) ----------------------------- 
#>    Basis functions: 5 
#>    Span: 24 s

# B-spline basis with N=11 basis functions, degree=1 (linear -> tent functions):
# tents peak every 2 s, at 2, 4, ..., 22 s
hrf_bs_10_1 <- gen_hrf(hrf_bspline, N = 11, degree = 1, name = "Tent Set (N=11)")
print(hrf_bs_10_1)
#> -- HRF: Tent Set (N=11) ----------------------------------- 
#>    Basis functions: 11 
#>    Span: 24 s

Each basis function covers a different stretch of the 24-second window; their weighted sum can take almost any smooth shape. Every function is zero at 0 s and at 24 s, so any weighted combination starts and ends at baseline.

bspline_times <- seq(0, 24, by = 0.1)
plot_hrfs(hrf_bs_5_3, time = bspline_times,
          title = "Cubic B-spline basis (N = 5)")

Five cubic B-spline basis functions on 0 to 24 seconds, labelled B1 to B5 at their peaks, tiling the window from early to late.

plot_hrfs(hrf_bs_10_1, time = bspline_times,
          title = "Tent basis: linear B-splines (N = 11)")

Eleven piecewise-linear tent functions peaking every 2 seconds from 2 to 22 seconds, each overlapping its neighbours and all zero at 0 and 24 seconds.

Sine Basis Set

The hrf_sine function creates a basis set using sine waves of different frequencies. Overlaid, five sine waves are hard to read, so each gets its own panel:

hrf_sin_5 <- gen_hrf(hrf_sine, N = 5, name = "Sine Basis (N=5)")
print(hrf_sin_5)
#> -- HRF: Sine Basis (N=5) ---------------------------------- 
#>    Basis functions: 5 
#>    Span: 24 s
plot_hrfs(hrf_sin_5, time = bspline_times, layout = "stack",
          title = "Sine basis (N = 5)")

Five sine basis functions in separate panels, with one to five full cycles over the 24 second window.

Half-Cosine Basis Set (FLOBS-like)

The hrf_half_cosine function implements the half-cosine HRF described by Woolrich et al. (2004), the model behind FSL’s FLOBS (FMRIB’s Linear Optimal Basis Sets). Four half-cosine segments of durations h1–h4 model the initial dip, rise, fall, and recovery; f1 and f2 set the height of the initial dip and the undershoot. The defaults set both to zero, giving a single bump; negative values add the dip and undershoot:

hrf_hc_default <- gen_hrf(hrf_half_cosine, name = "Half-cosine (default)")
hrf_hc_dip <- gen_hrf(hrf_half_cosine, f1 = -0.1, f2 = -0.2,
                      name = "Half-cosine (f1 = -0.1, f2 = -0.2)")
plot_hrfs(hrf_hc_default, hrf_hc_dip,
          labels = c("defaults (f1 = f2 = 0)", "f1 = -0.1, f2 = -0.2"),
          time = time_points,
          title = "Half-cosine HRF",
          subtitle = "Woolrich et al. (2004)")

Two half-cosine HRFs. The default has no dip or undershoot; with f1 = -0.1 and f2 = -0.2 an initial dip to -0.1 and an undershoot to -0.2 at about 13 seconds appear.

Other HRF Shapes

fmrihrf has several other parametric shapes. Each is created with gen_hrf():

# Gamma probability density
hrf_gam <- gen_hrf(hrf_gamma, shape = 6, rate = 1, name = "Gamma (shape=6, rate=1)")

# Mexican hat wavelet (second derivative of a Gaussian)
hrf_mh <- gen_hrf(hrf_mexhat, mean = 6, sd = 1.5, name = "Mexican Hat (mean=6, sd=1.5)")

# Difference of two inverse-logit (sigmoid) functions: separate rise and fall
hrf_il <- gen_hrf(hrf_inv_logit, mu1 = 5, s1 = 1, mu2 = 15, s2 = 1.5, name = "Inv. Logit Diff.")

The gamma peaks at (shape - 1) / rate = 5 s; the Mexican hat has negative lobes on both sides of its peak; the inverse-logit difference rises around mu1 and falls around mu2, giving a plateau:

plot_hrfs(hrf_gam, hrf_mh, hrf_il,
          labels = c("Gamma (shape 6, rate 1)", "Mexican hat (mean 6, sd 1.5)",
                     "Inverse-logit (rise 5 s, fall 15 s)"),
          time = time_points, layout = "stack",
          title = "Other HRF shapes")

Three HRF shapes in separate panels: a gamma peaking at 5 seconds, a Mexican hat with negative side lobes around a peak at 6 seconds, and an inverse-logit difference that rises around 5 seconds and falls around 15 seconds.

Boxcar and Weighted HRFs (No Hemodynamic Delay)

Traditional HRFs model the hemodynamic delay—the sluggish blood flow response that peaks several seconds after neural activity. However, sometimes you want to extract signal from specific time windows without assuming any hemodynamic transformation. This is useful for:

Simple Boxcar HRF (hrf_boxcar)

The hrf_boxcar function creates a simple step function that is constant within a time window and zero outside. Unlike traditional HRFs, there is no built-in hemodynamic delay—the HRF starts at time 0 (event onset) and extends for the specified width.

# Create a boxcar of width 5 seconds (from 0 to 5 seconds)
hrf_box <- hrf_boxcar(width = 5)
print(hrf_box)
#> -- HRF: boxcar[5] ----------------------------------------- 
#>    Basis functions: 1 
#>    Span: 5 s
#>    Parameters: width = 5, amplitude = 1, normalize = FALSE

To create a boxcar that starts at a later time point—useful for capturing signal in a specific post-stimulus window—use lag_hrf():

# Boxcar from 4-8 seconds post-stimulus (capturing the expected BOLD peak)
# Use lag_hrf() to delay a 4-second boxcar by 4 seconds
hrf_delayed <- hrf_boxcar(width = 4) %>% lag_hrf(lag = 4)
plot_hrfs(hrf_box, hrf_delayed,
          labels = c("Boxcar, 0-5 s", "Boxcar, 4-8 s"),
          normalize = TRUE, time = seq(0, 25, by = 0.05),
          reference = HRF_SPMG1, reference_label = "SPM canonical (peak = 1)",
          title = "Boxcars select a time window",
          subtitle = "No hemodynamic delay: a boxcar weights only its window")

A 0 to 5 second boxcar, a 4 to 8 second boxcar, and the normalized SPM canonical HRF as a dashed grey reference. The delayed boxcar covers the rising edge and peak of the canonical response.

What β Estimates: Boxcar Height

The height of a boxcar sets the units of its regression coefficient. With the default amplitude of 1, an isolated event’s β estimates the mean signal in the window. With normalize = TRUE, the boxcar is scaled to unit area (height 1/width), so β estimates the integrated signal over the window: the mean multiplied by the window width.

# Unit-area boxcar: a 4-second window lagged by 4 seconds (4-8 s)
hrf_norm <- hrf_boxcar(width = 4, normalize = TRUE) %>% lag_hrf(lag = 4)

# Check: amplitude should be 1/4 = 0.25
t_fine <- seq(0, 12, by = 0.01)
resp_norm <- evaluate(hrf_norm, t_fine)
cat("Amplitude of normalized boxcar:", max(resp_norm), "\n")
#> Amplitude of normalized boxcar: 0.25
cat("Expected (1/width):", 1/4, "\n")
#> Expected (1/width): 0.25

# Verify integral ≈ 1
integral <- sum(resp_norm) * 0.01
cat("Integral of normalized boxcar:", round(integral, 3), "\n")
#> Integral of normalized boxcar: 1

A quick least-squares check with a signal of 5 inside the 4-8 s window:

scan_t <- seq(0, 40, by = 0.5)
y <- ifelse(scan_t >= 4 & scan_t < 8, 5, 0)
beta <- function(hrf) {
  x <- evaluate(regressor(0, hrf), scan_t)
  unname(coef(lm(y ~ x))["x"])
}
c(height_1 = beta(hrf_boxcar(width = 4) %>% lag_hrf(lag = 4)),  # about 5: the mean
  unit_area = beta(hrf_norm))                                   # about 20: 5 x 4 s
#>  height_1 unit_area 
#>  4.861713 19.446853

Use the default height when you want β in signal units (the window mean); use normalize = TRUE when you want the window’s integrated response.

Weighted HRF (hrf_weighted)

The hrf_weighted function provides more flexibility by allowing you to specify different weights at different time points. You can either:

This creates either a step function (method = "constant") or a smoothly interpolated function (method = "linear"). With method = "constant", each weight fills one time bin: with width, the window is split into equal bins; with times, weight i runs from times[i] to times[i + 1], and the last bin is as wide as the one before it.

# 6 weights over a 10-second window: six equal bins of 10/6 s
hrf_wt_width <- hrf_weighted(
  weights = c(0.1, 0.3, 1.0, 1.0, 0.3, 0.1),
  width = 10,
  method = "constant"
)

# The same weights in 2-second bins starting at 2, 4, ..., 12 s (window 2-14 s)
hrf_wt <- hrf_weighted(
  weights = c(0.1, 0.3, 1.0, 1.0, 0.3, 0.1),
  times = c(2, 4, 6, 8, 10, 12),
  method = "constant"
)

# Smooth weights using linear interpolation between the same time points
hrf_smooth <- hrf_weighted(
  weights = c(0.1, 0.3, 1.0, 1.0, 0.3, 0.1),
  times = c(2, 4, 6, 8, 10, 12),
  method = "linear"
)
plot_hrfs(hrf_wt_width, hrf_wt, hrf_smooth,
          labels = c("width = 10: six bins, 0-10 s",
                     "times = 2, ..., 12: bins, 2-14 s",
                     "method = linear: 2-12 s"),
          time = seq(0, 16, by = 0.02), layout = "stack",
          title = "Weighted HRFs", subtitle = "Weights 0.1, 0.3, 1, 1, 0.3, 0.1")

Three weighted HRFs in separate panels: steps from 0 to 10 seconds, the same steps from 2 to 12 seconds, and a linear interpolation through the same weights from 2 to 12 seconds.

Sub-second Precision

The hrf_weighted function supports sub-second time intervals, which is useful for fine-grained temporal weighting:

# Sub-second intervals: create a Gaussian-shaped weight function
times_fine <- seq(4, 10, by = 0.25)
weights_gaussian <- dnorm(times_fine, mean = 7, sd = 1)

hrf_gauss_wt <- hrf_weighted(weights_gaussian, times = times_fine, method = "linear")

plot_hrfs(hrf_gauss_wt, labels = "Gaussian weights, 0.25 s spacing",
          time = seq(0, 14, by = 0.02),
          title = "Weights at sub-second spacing")

A Gaussian-shaped weighting function centred at 7 seconds, built from weights every 0.25 seconds between 4 and 10 seconds.

Normalized Weighted HRF

When normalize = TRUE, the weights are scaled to sum to 1 (method = "constant") or the curve to integrate to 1 (method = "linear"). This fixes the scale of the weighting profile, which makes coefficients comparable across profiles. It does not turn β into a weighted mean: least squares fits the amplitude of the whole profile, β = Σ w(t) y(t) / Σ w(t)², which equals a plain mean only for a boxcar. To summarise a window as a weighted mean, compute that mean directly from the data.

hrf_wt_norm <- hrf_weighted(
  weights = c(1, 2, 2, 1),  # Will be normalized
  times = c(4, 6, 8, 10),
  method = "constant",
  normalize = TRUE
)

# Four 2-second bins (4-6, 6-8, 8-10, 10-12 s); weights 1, 2, 2, 1 are
# rescaled to 1/6, 2/6, 2/6, 1/6
evaluate(hrf_wt_norm, c(5, 7, 9, 11, 13))
#> [1] 0.1666667 0.3333333 0.3333333 0.1666667 0.0000000

Practical Example: Comparing Early vs. Late Response Windows

A common analysis compares BOLD signal in early vs. late portions of a trial. Here’s how to set up HRFs for this:

# Early window: 2-6 seconds (4-second boxcar lagged by 2 seconds)
hrf_early <- hrf_boxcar(width = 4) %>% lag_hrf(lag = 2)

# Late window: 8-12 seconds (4-second boxcar lagged by 8 seconds)
hrf_late <- hrf_boxcar(width = 4) %>% lag_hrf(lag = 8)

These boxcars have height 1, so each β is the mean signal in its window. Drawn against the canonical response (scaled to the same height), the early window covers the rise and peak, and the late window the return to baseline:

spm_ref <- gen_empirical_hrf(time_points,
  HRF_SPMG1(time_points) / max(HRF_SPMG1(time_points)))
plot_hrfs(hrf_early, hrf_late,
          labels = c("Early window, 2-6 s", "Late window, 8-12 s"),
          time = seq(0, 25, by = 0.05),
          reference = spm_ref, reference_label = "SPM canonical, scaled to peak 1",
          title = "Early and late windows",
          subtitle = "Each beta is the mean signal in its window")

Early (2 to 6 seconds) and late (8 to 12 seconds) boxcar windows of height 1, drawn with the SPM canonical HRF scaled to the same height. The early window covers the peak; the late window covers the decline.

Using these HRFs in separate regressors allows you to estimate and compare the mean BOLD signal in each window.

Using Boxcar/Weighted HRFs with Regressors

These HRFs integrate seamlessly with the regressor() function:

# Create a regressor with boxcar HRF (4-second window starting 4s after onset)
reg_boxcar <- regressor(
  onsets = c(0, 20, 40),
  hrf = hrf_boxcar(width = 4, normalize = TRUE) %>% lag_hrf(lag = 4)
)

# Compare with traditional SPM HRF
reg_spm <- regressor(onsets = c(0, 20, 40), hrf = HRF_SPMG1)

The two regressors have different units, so each gets its own panel (and y axis). Grey bars under each panel mark the event onsets. The small notches at 24 and 44 s are where the previous event’s canonical response is cut off at the end of its 24-second span (see the regressor vignette):

plot_regressors(reg_spm, reg_boxcar,
                labels = c("SPM canonical HRF", "Boxcar HRF, 4-8 s window"),
                grid = seq(0, 60, by = 0.05), layout = "stack",
                title = "Boxcar vs. canonical regressor",
                subtitle = "Events at t = 0, 20, 40 s")

Two regressors for events at 0, 20 and 40 seconds, in separate panels: the smooth SPM canonical regressor, and a boxcar regressor with 4 second plateaus from 4 to 8 seconds after each event.

Creating Custom Basis Sets with hrf_set

The hrf_set function allows you to combine any set of HRF functions into a single multivariate HRF object (a basis set).

For example, we can create a basis set from a series of lagged Gaussian HRFs:

# Create a list of lagged Gaussian HRFs
lag_times <- seq(0, 10, by = 2)
list_of_hrfs <- lapply(lag_times, function(lag) {
  lag_hrf(HRF_GAUSSIAN, lag = lag)
})

# Combine them into a single HRF basis set object
hrf_custom_set <- do.call(hrf_set, list_of_hrfs)
print(hrf_custom_set) # Note: name is default 'hrf_set', nbasis is 6
#> -- HRF: hrf_set ------------------------------------------- 
#>    Basis functions: 6 
#>    Span: 34 s
plot_hrfs(hrf_custom_set, time = time_points,
          title = "Lagged-Gaussian basis set",
          subtitle = "HRF_GAUSSIAN lagged by 0, 2, ..., 10 s")

Six Gaussian basis functions, labelled B1 to B6, peaking every 2 seconds from 6 to 16 seconds.

Creating Empirical HRFs

From a Single Measured Response (gen_empirical_hrf)

If you have a measured or estimated hemodynamic response profile (e.g., from deconvolution), you can turn it into an HRF function using gen_empirical_hrf. It uses linear interpolation between the provided points.

# Simulate an average measured response profile
sim_times <- 0:24
set.seed(42) # For reproducibility
sim_profile <- rowMeans(replicate(20, {
  h <- HRF_SPMG1 %>% lag_hrf(lag = runif(n = 1, min = -1, max = 1)) %>%
                    block_hrf(width = runif(n = 1, min = 0, max = 2))
  h(sim_times)
}))

# Normalize profile to max = 1 for better visualization
sim_profile_norm <- sim_profile / max(sim_profile)

# Create the empirical HRF function from the normalized profile
emp_hrf <- gen_empirical_hrf(sim_times, sim_profile_norm)
print(emp_hrf)
#> -- HRF: empirical_hrf ------------------------------------- 
#>    Basis functions: 1 
#>    Span: 24 s

The points are the measured profile; the line is the interpolating HRF:

Empirical HRF drawn as a line through 25 measured points, one per second, peaking near 6 seconds with an undershoot after 12 seconds.

Empirical Basis Set via PCA

You can create an empirical basis set by applying dimensionality reduction (like PCA) to a collection of observed or simulated HRFs.

# 1. Simulate a matrix of diverse HRFs
set.seed(123) # for reproducibility
n_sim <- 50
sim_mat <- replicate(n_sim, {
  hrf_func <- HRF_SPMG1 %>%
              lag_hrf(lag = runif(1, -2, 2)) %>%
              block_hrf(width = runif(1, 0, 3))
  hrf_func(sim_times)
})

The 50 simulated responses vary in latency and width; PCA finds the few shapes that capture most of this variation:

Fifty simulated HRFs drawn as thin grey lines, with their mean as a thick line. Peak times range from about 3 to 8 seconds.

# 2. Perform PCA on the transpose (each column = one HRF, each row = one time point)
pca_res <- prcomp(t(sim_mat), center = TRUE, scale. = FALSE)
n_components <- 3

# Print variance explained by top components
variance_explained <- summary(pca_res)$importance[2, 1:n_components]
cat("Variance explained by top", n_components, "components:",
    paste0(round(variance_explained * 100, 1), "%"), "\n")
#> Variance explained by top 3 components: 67.3% 30% 2.5%

# Extract the top principal components
pc_vectors <- pca_res$rotation[, 1:n_components]

# 3. Convert principal components into HRF functions
list_pc_hrfs <- list()

for (i in 1:n_components) {
  pc_vec <- pc_vectors[, i]
  # The sign of a principal component is arbitrary; make its largest
  # deviation positive so PC1 looks like an HRF rather than its mirror image
  pc_vec <- pc_vec * sign(pc_vec[which.max(abs(pc_vec - pc_vec[1]))] - pc_vec[1])
  pc_vec_zeroed <- pc_vec - pc_vec[1]
  max_abs <- max(abs(pc_vec_zeroed))
  pc_vec_norm <- pc_vec_zeroed / max_abs
  list_pc_hrfs[[i]] <- gen_empirical_hrf(sim_times, pc_vec_norm)
}

# 4. Combine PC HRFs into a basis set using hrf_set
emp_pca_basis <- do.call(hrf_set, list_pc_hrfs)
print(emp_pca_basis)
#> -- HRF: hrf_set ------------------------------------------- 
#>    Basis functions: 3 
#>    Span: 24 s
plot_hrfs(emp_pca_basis,
          labels = paste0("PC", 1:n_components, " (",
                          round(variance_explained * 100), "% of variance)"),
          time = seq(0, 24, by = 0.1),
          title = "Empirical basis from PCA",
          subtitle = "Each component scaled to a maximum absolute value of 1")

The first three principal components of the simulated HRFs, each scaled to a maximum absolute value of one. PC1 resembles the average HRF; PC2 and PC3 have alternating positive and negative lobes that shift and reshape the response.

This empirical basis set can then be used in regression models just like any other pre-defined or custom basis set.