--- title: Hemodynamic Response Functions author: Bradley R. Buchsbaum date: '`r Sys.Date()`' output: albersdown::albers_vignette: family: red preset: interaction toc: yes toc_depth: 2.0 vignette: | %\VignetteIndexEntry{Hemodynamic Response Functions} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", # Sharper website figures; keep the CRAN vignette archive compact. fig.retina = if (identical(Sys.getenv("IN_PKGDOWN"), "true")) 2 else 1, fig.width = 7, fig.height = 4, message = FALSE, warning = FALSE ) # CRAN builds: skip dark-mode figure twins to keep the source package small; # the pkgdown site (IN_PKGDOWN = "true") keeps them. if (!identical(Sys.getenv("IN_PKGDOWN"), "true")) { options(albersdown.dark_figures = FALSE) } library(fmrihrf) library(dplyr) # For pipe operator %>% library(ggplot2) # For plotting library(tidyr) # For data manipulation ``` ## 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`). ```{r print_hrfs} # SPM canonical HRF (based on difference of two gamma functions) print(HRF_SPMG1) # Gaussian HRF print(HRF_GAUSSIAN) ``` 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()`): ```{r evaluate_basic_hrfs, fig.alt="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."} 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") ``` 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. ```{r fixed-hrf-scale} spm_scaled <- gen_hrf(HRF_SPMG1, hrf_norm = "spm") spm_grid <- seq(0, 32, length.out = 1600) sum(spm_scaled(spm_grid)) ``` ```{r check-fixed-hrf-scale, include = FALSE} stopifnot(abs(sum(spm_scaled(spm_grid)) - 1) < 1e-12) ``` 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`). ```{r modify_gaussian_params} # 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)") ``` ```{r modify_gaussian_params_plot, fig.alt="Three Gaussian HRFs. Later means peak later, and larger standard deviations give lower, broader curves."} 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") ``` `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. ```{r blocked_hrfs} # 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. ```{r blocked_hrfs_plot, fig.alt="SPM canonical responses to 1, 2 and 4 second events. Longer events produce higher and later peaks."} 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)") ``` ### 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. ```{r blocked_normalized} # 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) ``` ```{r blocked_normalized_plot, fig.alt="Normalized SPM responses to 1, 2 and 4 second events. All peak at 1; longer events peak later and are broader."} 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)") ``` ### 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. ```{r blocked_summate_false} # 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) ``` ```{r blocked_summate_false_plot, fig.alt="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."} 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)") ``` `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. ```{r lagged_hrfs} # 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) ``` ```{r lagged_hrfs_plot, fig.alt="Gaussian HRF shifted by -2, 0 and +3 seconds. The shape is unchanged; only the peak time moves."} 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)") ``` ## Combining Lag and Duration We can combine `lag_hrf` and `block_hrf` using the pipe operator (`%>%`) from `dplyr` (or `magrittr`). ```{r lagged_blocked_hrfs} # 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) ``` ```{r lagged_blocked_hrfs_plot, fig.alt="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."} 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)") ``` Alternatively, `gen_hrf` can apply lag and width directly, and gives the same function: ```{r gen_hrf_lag_width} 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))) ``` ## 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`). ```{r spm_basis_sets} # SPM + Temporal Derivative (2 basis functions) print(HRF_SPMG2) # SPM + Temporal + Dispersion Derivatives (3 basis functions) print(HRF_SPMG3) ``` `HRF_SPMG2` is the first two columns of `HRF_SPMG3`, so one figure shows both: ```{r spm_basis_sets_plot, fig.height = 4.4, fig.alt="Canonical SPM response and its temporal and dispersion derivatives, with all positive and negative values retained."} plot_hrfs(HRF_SPMG3, labels = c("Canonical", "Temporal derivative", "Dispersion derivative"), time = time_points, title = "SPM canonical and derivatives") ``` 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: ```{r spm_derivative_shift, fig.height = 4.4, fig.alt="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."} 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") ``` ### 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`. ```{r bspline_basis} # 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) # 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) ``` 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. ```{r bspline_basis_plot, fig.alt="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."} bspline_times <- seq(0, 24, by = 0.1) plot_hrfs(hrf_bs_5_3, time = bspline_times, title = "Cubic B-spline basis (N = 5)") ``` ```{r tent_basis_plot, fig.alt="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."} plot_hrfs(hrf_bs_10_1, time = bspline_times, title = "Tent basis: linear B-splines (N = 11)") ``` ### 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: ```{r sine_basis} hrf_sin_5 <- gen_hrf(hrf_sine, N = 5, name = "Sine Basis (N=5)") print(hrf_sin_5) ``` ```{r sine_basis_plot, fig.height = 5, fig.alt="Five sine basis functions in separate panels, with one to five full cycles over the 24 second window."} plot_hrfs(hrf_sin_5, time = bspline_times, layout = "stack", title = "Sine basis (N = 5)") ``` ### 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: ```{r half_cosine, fig.alt="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."} 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)") ``` ## Other HRF Shapes `fmrihrf` has several other parametric shapes. Each is created with `gen_hrf()`: ```{r other_shapes} # 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: ```{r other_shapes_plot, fig.height = 5, fig.alt="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."} 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") ``` ## 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: - Extracting raw signal averages from specific post-stimulus windows - Trial-wise analyses where you want the mean (or weighted mean) of the BOLD signal - Comparing signal in different temporal windows directly - Creating custom temporal weighting schemes ### 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`. ```{r boxcar_basic} # Create a boxcar of width 5 seconds (from 0 to 5 seconds) hrf_box <- hrf_boxcar(width = 5) print(hrf_box) ``` To create a boxcar that starts at a later time point—useful for capturing signal in a specific post-stimulus window—use `lag_hrf()`: ```{r boxcar_delayed} # 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) ``` ```{r boxcar_delayed_plot, fig.alt="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."} 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") ``` ### 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. ```{r boxcar_normalized} # 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") cat("Expected (1/width):", 1/4, "\n") # Verify integral ≈ 1 integral <- sum(resp_norm) * 0.01 cat("Integral of normalized boxcar:", round(integral, 3), "\n") ``` A quick least-squares check with a signal of 5 inside the 4-8 s window: ```{r boxcar_beta_check} 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 ``` 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: - Use `width` + `weights`: evenly space the weights across the specified width - Use `times` + `weights`: explicitly specify the time points for each weight 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. ```{r weighted_width} # 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" ) ``` ```{r weighted_plot, fig.height = 5.2, fig.alt="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."} 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") ``` #### Sub-second Precision The `hrf_weighted` function supports sub-second time intervals, which is useful for fine-grained temporal weighting: ```{r weighted_subsecond, fig.alt="A Gaussian-shaped weighting function centred at 7 seconds, built from weights every 0.25 seconds between 4 and 10 seconds."} # 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") ``` ### 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. ```{r weighted_normalized} 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)) ``` ### 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: ```{r early_late_comparison} # 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: ```{r early_late_comparison_plot, fig.alt="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."} 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") ``` 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: ```{r boxcar_regressor} # 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): ```{r boxcar_regressor_plot, fig.height = 4.4, fig.alt="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."} 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") ``` ## 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: ```{r custom_basis_lagged} # 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 ``` ```{r custom_basis_lagged_plot, fig.alt="Six Gaussian basis functions, labelled B1 to B6, peaking every 2 seconds from 6 to 16 seconds."} plot_hrfs(hrf_custom_set, time = time_points, title = "Lagged-Gaussian basis set", subtitle = "HRF_GAUSSIAN lagged by 0, 2, ..., 10 s") ``` ## 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. ```{r empirical_hrf_single} # 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) ``` The points are the measured profile; the line is the interpolating HRF: ```{r empirical_hrf_single_plot, echo = FALSE, fig.alt="Empirical HRF drawn as a line through 25 measured points, one per second, peaking near 6 seconds with an undershoot after 12 seconds."} emp_df <- data.frame(time = seq(0, 24, by = 0.1)) emp_df$response <- emp_hrf(emp_df$time) ggplot(emp_df, aes(time, response)) + geom_hline(yintercept = 0, colour = "grey70", linewidth = 0.3) + geom_line(colour = hrf_palette(1), linewidth = 0.9) + geom_point(data = data.frame(time = sim_times, response = sim_profile_norm), colour = hrf_palette(1), size = 1.8) + labs(title = "Empirical HRF from a profile", x = "Time (s)", y = "Response / peak") + theme(plot.title.position = "plot") ``` ### 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. ```{r empirical_hrf_pca} # 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: ```{r empirical_hrf_pca_plot1, echo = FALSE, fig.alt="Fifty simulated HRFs drawn as thin grey lines, with their mean as a thick line. Peak times range from about 3 to 8 seconds."} sim_df <- data.frame(time = rep(sim_times, n_sim), response = as.vector(sim_mat), id = rep(seq_len(n_sim), each = length(sim_times))) ggplot(sim_df, aes(time, response)) + geom_hline(yintercept = 0, colour = "grey70", linewidth = 0.3) + geom_line(aes(group = id), colour = "grey60", linewidth = 0.3, alpha = 0.6) + geom_line(data = data.frame(time = sim_times, response = rowMeans(sim_mat)), colour = hrf_palette(1), linewidth = 1.1) + labs(title = "50 simulated HRFs and their mean", subtitle = "Random lag and duration", x = "Time (s)", y = "Response") + theme(plot.title.position = "plot") ``` ```{r empirical_hrf_pca2} # 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") # 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) ``` ```{r empirical_hrf_pca_plot2, fig.alt="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."} 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") ``` This empirical basis set can then be used in regression models just like any other pre-defined or custom basis set.