| Title: | Hemodynamic Response Functions for fMRI Data Analysis |
|---|---|
| Description: | Creates, manipulates, and evaluates hemodynamic response functions and event-related regressors for functional magnetic resonance imaging data analysis. Supports multiple basis sets including Canonical, Gamma, Gaussian, B-spline, and Fourier bases. Features decorators for time-shifting and blocking, and efficient convolution algorithms for regressor construction. Methods are based on standard fMRI analysis techniques as described in Jezzard et al. (2001, ISBN:9780192630711). |
| Authors: | Bradley Buchsbaum [aut, cre] |
| Maintainer: | Bradley Buchsbaum <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 0.4.0 |
| Built: | 2026-10-02 23:16:08 UTC |
| Source: | https://github.com/bbuchsbaum/fmrihrf |
Calculate the onset time in seconds for each fMRI volume acquisition from the start of the experiment.
acquisition_onsets(x, ...) ## S3 method for class 'sampling_frame' acquisition_onsets(x, ...)acquisition_onsets(x, ...) ## S3 method for class 'sampling_frame' acquisition_onsets(x, ...)
x |
A sampling_frame object |
... |
Additional arguments (for extensibility) |
Returns the temporal onset of each brain volume acquisition, accounting
for TR, start_time, and run structure. This is essentially a convenience
wrapper around samples(x, global = TRUE) that provides clearer
semantic meaning for the common use case of getting acquisition times.
Note: The onset times include the start_time offset (default TR/2), so the first acquisition typically doesn't start at 0.
Numeric vector of acquisition onset times in seconds
samples for more flexible timing queries
# Single block with default start_time (TR/2 = 1) sf <- sampling_frame(blocklens = 100, TR = 2) onsets <- acquisition_onsets(sf) head(onsets) # Returns: 1, 3, 5, 7, 9, 11, ... # Multiple blocks with same TR sf2 <- sampling_frame(blocklens = c(100, 120), TR = 2) onsets2 <- acquisition_onsets(sf2) # First block: 1, 3, 5, ..., 199 # Second block: 201, 203, 205, ..., 439 # Variable TR per block sf3 <- sampling_frame(blocklens = c(100, 100), TR = c(2, 1.5)) onsets3 <- acquisition_onsets(sf3) # First block: 1, 3, 5, ..., 199 (TR=2) # Second block: 200.75, 202.25, 203.75, ... (TR=1.5, start_time=0.75) # Custom start times sf4 <- sampling_frame(blocklens = c(50, 50), TR = 2, start_time = 0) onsets4 <- acquisition_onsets(sf4) head(onsets4) # Returns: 0, 2, 4, 6, 8, 10, ...# Single block with default start_time (TR/2 = 1) sf <- sampling_frame(blocklens = 100, TR = 2) onsets <- acquisition_onsets(sf) head(onsets) # Returns: 1, 3, 5, 7, 9, 11, ... # Multiple blocks with same TR sf2 <- sampling_frame(blocklens = c(100, 120), TR = 2) onsets2 <- acquisition_onsets(sf2) # First block: 1, 3, 5, ..., 199 # Second block: 201, 203, 205, ..., 439 # Variable TR per block sf3 <- sampling_frame(blocklens = c(100, 100), TR = c(2, 1.5)) onsets3 <- acquisition_onsets(sf3) # First block: 1, 3, 5, ..., 199 (TR=2) # Second block: 200.75, 202.25, 203.75, ... (TR=1.5, start_time=0.75) # Custom start times sf4 <- sampling_frame(blocklens = c(50, 50), TR = 2, start_time = 0) onsets4 <- acquisition_onsets(sf4) head(onsets4) # Returns: 0, 2, 4, 6, 8, 10, ...
Generic accessor returning event amplitudes or scaling factors.
amplitudes(x, ...) ## S3 method for class 'Reg' amplitudes(x, ...)amplitudes(x, ...) ## S3 method for class 'Reg' amplitudes(x, ...)
x |
Object containing amplitude information |
... |
Additional arguments passed to methods |
Numeric vector of amplitudes
# Create a regressor with varying amplitudes reg <- regressor(onsets = c(1, 5, 10), hrf = HRF_SPMG1, amplitude = c(1, 0.5, 2), span = 20) amplitudes(reg)# Create a regressor with varying amplitudes reg <- regressor(onsets = c(1, 5, 10), hrf = HRF_SPMG1, amplitude = c(1, 0.5, 2), span = 20) amplitudes(reg)
Creates a new HRF object representing a response to a sustained (blocked) stimulus by convolving the input HRF with a boxcar function of a given width.
block_hrf( hrf, width, precision = 0.1, half_life = Inf, summate = TRUE, normalize = FALSE )block_hrf( hrf, width, precision = 0.1, half_life = Inf, summate = TRUE, normalize = FALSE )
hrf |
The HRF object (of class 'HRF') to block. |
width |
The width of the block in seconds. |
precision |
The sampling precision in seconds used for the internal convolution (default: 0.1). |
half_life |
The half-life of an optional exponential decay applied during the block (default: Inf, meaning no decay). |
summate |
Logical; if TRUE (default), responses within the block are integrated (summed). If FALSE, the integrated response is divided by the total block weight so amplitude does not grow with block width. |
normalize |
Logical; if TRUE, the resulting blocked HRF is scaled so that its peak value is 1 (default: FALSE). |
A new HRF object representing the blocked function.
Other HRF_decorator_functions:
lag_hrf(),
normalise_hrf(),
normalize_hrf()
blocked_spmg1 <- block_hrf(HRF_SPMG1, width = 5) t_vals <- seq(0, 30, by = 0.5) plot(t_vals, HRF_SPMG1(t_vals), type = 'l', col = "blue", ylab = "Response", xlab = "Time") lines(t_vals, blocked_spmg1(t_vals), col = "red") legend("topright", legend = c("Original", "Blocked (width=5)"), col = c("blue", "red"), lty = 1)blocked_spmg1 <- block_hrf(HRF_SPMG1, width = 5) t_vals <- seq(0, 30, by = 0.5) plot(t_vals, HRF_SPMG1(t_vals), type = 'l', col = "blue", ylab = "Response", xlab = "Time") lines(t_vals, blocked_spmg1(t_vals), col = "red") legend("topright", legend = c("Original", "Blocked (width=5)"), col = c("blue", "red"), lty = 1)
Generic accessor returning block indices for each sample or onset.
blockids(x, ...) ## S3 method for class 'sampling_frame' blockids(x, ...)blockids(x, ...) ## S3 method for class 'sampling_frame' blockids(x, ...)
x |
Object containing block structure |
... |
Additional arguments passed to methods |
Integer vector of block ids
# Get block identifiers from a sampling frame sframe <- sampling_frame(blocklens = c(100, 120, 80), TR = 2) blockids(sframe)# Get block identifiers from a sampling frame sframe <- sampling_frame(blocklens = c(100, 120, 80), TR = 2) blockids(sframe)
Generic accessor returning the number of scans in each block of a sampling frame or similar object.
blocklens(x, ...) ## S3 method for class 'sampling_frame' blocklens(x, ...)blocklens(x, ...) ## S3 method for class 'sampling_frame' blocklens(x, ...)
x |
Object containing block length information |
... |
Additional arguments passed to methods |
Numeric vector of block lengths
# Get block lengths from a sampling frame sframe <- sampling_frame(blocklens = c(100, 120, 80), TR = 2) blocklens(sframe)# Get block lengths from a sampling frame sframe <- sampling_frame(blocklens = c(100, 120, 80), TR = 2) blocklens(sframe)
Calculates the derivative of a Hemodynamic Response Function (HRF) at specified time points. This is useful for:
Understanding HRF dynamics and rate of change
Creating temporal derivative regressors for fMRI models
Analyzing HRF shape characteristics
Implementing advanced HRF basis sets
deriv(x, t, ...)deriv(x, t, ...)
x |
An HRF object |
t |
Numeric vector of time points at which to evaluate the derivative |
... |
Additional arguments passed to specific methods |
The derivative computation method depends on the HRF type:
Analytic derivatives are used when available (e.g., SPMG1, SPMG2, SPMG3)
Numeric finite-difference approximation is used as fallback
The default implementation uses numDeriv::grad for numerical
differentiation when analytic derivatives are not available.
Numeric vector or matrix of derivative values at the specified time points. For multi-basis HRFs, returns a matrix with one column per basis function.
[evaluate()], [HRF_objects], [numDeriv::grad()]
Other hrf:
HRF_objects,
penalty_matrix()
# Compute derivative of SPM canonical HRF t <- seq(0, 20, by = 0.1) hrf_deriv <- deriv(HRF_SPMG1, t) # Plot HRF and its derivative hrf_vals <- evaluate(HRF_SPMG1, t) plot(t, hrf_vals, type = "l", col = "black", ylab = "Response", xlab = "Time (s)") lines(t, hrf_deriv, col = "red", lty = 2) legend("topright", c("HRF", "Derivative"), col = c("black", "red"), lty = c(1, 2)) # For multi-basis HRFs, returns matrix deriv_matrix <- deriv(HRF_SPMG3, t) # Returns derivatives for all 3 basis functions# Compute derivative of SPM canonical HRF t <- seq(0, 20, by = 0.1) hrf_deriv <- deriv(HRF_SPMG1, t) # Plot HRF and its derivative hrf_vals <- evaluate(HRF_SPMG1, t) plot(t, hrf_vals, type = "l", col = "black", ylab = "Response", xlab = "Time (s)") lines(t, hrf_deriv, col = "red", lty = 2) legend("topright", c("HRF", "Derivative"), col = c("black", "red"), lty = c(1, 2)) # For multi-basis HRFs, returns matrix deriv_matrix <- deriv(HRF_SPMG3, t) # Returns derivatives for all 3 basis functions
Uses numerical differentiation via numDeriv::grad when analytic derivatives are not available for a specific HRF type.
Uses the analytic derivative formula for the SPM canonical HRF.
Returns derivatives for both the canonical HRF and its temporal derivative. The first column contains the derivative of the canonical HRF, and the second column contains the second derivative (derivative of the temporal derivative).
Returns derivatives for the canonical HRF and its two derivatives. Returns the first and second time derivatives of the canonical, followed by the time derivative of the dispersion basis.
## S3 method for class 'HRF' deriv(x, t, ...) ## S3 method for class 'SPMG1_HRF' deriv(x, t, ...) ## S3 method for class 'SPMG2_HRF' deriv(x, t, ...) ## S3 method for class 'SPMG3_HRF' deriv(x, t, ...)## S3 method for class 'HRF' deriv(x, t, ...) ## S3 method for class 'SPMG1_HRF' deriv(x, t, ...) ## S3 method for class 'SPMG2_HRF' deriv(x, t, ...) ## S3 method for class 'SPMG3_HRF' deriv(x, t, ...)
x |
An SPMG3_HRF object |
t |
Numeric vector of time points at which to evaluate the derivative |
... |
Additional arguments (currently unused) |
Numeric vector or matrix of derivative values
Numeric vector of derivative values
Matrix with 2 columns of derivative values
Matrix with 3 columns of derivative values
t <- seq(0, 30, by = 0.5) d <- deriv(HRF_SPMG1, t)t <- seq(0, 30, by = 0.5) d <- deriv(HRF_SPMG1, t)
Get durations of an object
durations(x, ...) ## S3 method for class 'Reg' durations(x, ...)durations(x, ...) ## S3 method for class 'Reg' durations(x, ...)
x |
The object to get durations from |
... |
Additional arguments passed to methods |
A numeric vector of durations
# Create a regressor with event durations reg <- regressor(onsets = c(1, 5, 10), hrf = HRF_SPMG1, duration = c(2, 3, 1), span = 20) durations(reg)# Create a regressor with event durations reg <- regressor(onsets = c(1, 5, 10), hrf = HRF_SPMG1, duration = c(2, 3, 1), span = 20) durations(reg)
Generic function to evaluate a regressor object over a specified time grid. Different types of regressors may have different evaluation methods.
evaluate(x, grid, ...) ## S3 method for class 'Reg' evaluate( x, grid, precision = 0.33, method = c("conv", "loop", "fft", "Rconv"), sparse = FALSE, normalize = FALSE, ... )evaluate(x, grid, ...) ## S3 method for class 'Reg' evaluate( x, grid, precision = 0.33, method = c("conv", "loop", "fft", "Rconv"), sparse = FALSE, normalize = FALSE, ... )
x |
A 'Reg' object (or an object inheriting from it, like 'regressor'). |
grid |
Numeric vector specifying the time points (seconds) for evaluation. |
... |
Additional arguments passed down (e.g., to 'evaluate.HRF' in the loop method). |
precision |
Numeric sampling precision for internal HRF evaluation and convolution (seconds). |
method |
The evaluation method:
|
sparse |
Logical indicating whether to return a sparse matrix (from the Matrix package). Default is FALSE. |
normalize |
Logical; if TRUE, scale evaluated regressor output to unit peak (maximum absolute value of 1). For multi-basis regressors, each basis column is normalized independently. |
A numeric vector or matrix containing the evaluated regressor values
[single_trial_regressor()], [regressor()]
# Create a regressor reg <- regressor(onsets = c(10, 30, 50), hrf = HRF_SPMG1) # Evaluate at specific time points times <- seq(0, 80, by = 0.1) response <- evaluate(reg, times) # Plot the response plot(times, response, type = "l", xlab = "Time (s)", ylab = "Response") # Create a regressor reg <- regressor(onsets = c(10, 30, 50), hrf = HRF_SPMG1) # Evaluate with default method (conv) times <- seq(0, 80, by = 0.5) response <- evaluate(reg, times) # Try different evaluation methods response_loop <- evaluate(reg, times, method = "loop") # With higher precision response_precise <- evaluate(reg, times, precision = 0.1)# Create a regressor reg <- regressor(onsets = c(10, 30, 50), hrf = HRF_SPMG1) # Evaluate at specific time points times <- seq(0, 80, by = 0.1) response <- evaluate(reg, times) # Plot the response plot(times, response, type = "l", xlab = "Time (s)", ylab = "Response") # Create a regressor reg <- regressor(onsets = c(10, 30, 50), hrf = HRF_SPMG1) # Evaluate with default method (conv) times <- seq(0, 80, by = 0.5) response <- evaluate(reg, times) # Try different evaluation methods response_loop <- evaluate(reg, times, method = "loop") # With higher precision response_precise <- evaluate(reg, times, precision = 0.1)
This function evaluates a hemodynamic response function (HRF) object for a given set of time points (grid) and other parameters. It handles both point evaluation (duration=0) and block evaluation (duration > 0).
## S3 method for class 'HRF' evaluate( x, grid, amplitude = 1, duration = 0, precision = 0.2, summate = TRUE, normalize = FALSE, ... )## S3 method for class 'HRF' evaluate( x, grid, amplitude = 1, duration = 0, precision = 0.2, summate = TRUE, normalize = FALSE, ... )
x |
The HRF object (inherits from 'HRF' and 'function'). |
grid |
A numeric vector of time points at which to evaluate the HRF. |
amplitude |
The scaling value for the event (default: 1). |
duration |
The duration of the event (seconds). A duration of 0 is an impulse of unit mass and returns the HRF itself; a duration > 0 is a unit-height boxcar and returns its integral against the HRF (default: 0). |
precision |
The quadrature step used to integrate over the block when duration > 0 (default: 0.2). It controls numerical accuracy only: the result converges as 'precision' decreases and does not otherwise depend on it. |
summate |
Logical; whether the HRF response should accumulate over the duration (default: TRUE), giving the integral over a unit-height block, so peak amplitude grows with duration. If FALSE, the result is divided by the block duration, giving the duration-averaged (unit-mass) response: the temporal profile is preserved, peak amplitude does not grow with duration, and the value approaches the impulse response as duration approaches 0. |
normalize |
Logical; scale output so that the peak absolute value is 1 (default: FALSE). Applied *after* amplitude scaling and duration processing. |
... |
Additional arguments (unused). |
A numeric vector or matrix of HRF values at the specified time points.
# Evaluate canonical HRF at specific times times <- seq(0, 20, by = 0.5) response <- evaluate(HRF_SPMG1, times) # Evaluate with amplitude scaling response_scaled <- evaluate(HRF_SPMG1, times, amplitude = 2) # Evaluate with duration (block design) response_block <- evaluate(HRF_SPMG1, times, duration = 5, summate = TRUE) # Multi-basis HRF evaluation response_multi <- evaluate(HRF_SPMG3, times) # Returns 3-column matrix# Evaluate canonical HRF at specific times times <- seq(0, 20, by = 0.5) response <- evaluate(HRF_SPMG1, times) # Evaluate with amplitude scaling response_scaled <- evaluate(HRF_SPMG1, times, amplitude = 2) # Evaluate with duration (block design) response_block <- evaluate(HRF_SPMG1, times, duration = 5, summate = TRUE) # Multi-basis HRF evaluation response_multi <- evaluate(HRF_SPMG3, times) # Returns 3-column matrix
Creates a regressor from a continuously sampled feature (for example RMS
energy of an acoustic stimulus) by treating each sample as a zero-order-hold
bin of width . The result is the Riemann-sum / ZOH approximation
of convolving the feature with an HRF, not a train of unit-mass impulses.
feature_regressor( values, hrf = HRF_SPMG1, times = NULL, dt = NULL, start = 0, center = TRUE, scale = c("none", "sd"), mask = NULL, span = NULL )feature_regressor( values, hrf = HRF_SPMG1, times = NULL, dt = NULL, start = 0, center = TRUE, scale = c("none", "sd"), mask = NULL, span = NULL )
values |
Numeric vector of feature samples. Matrix and array inputs are rejected; construct one feature regressor per column instead. |
hrf |
The hemodynamic response function to convolve with the feature.
Same types as |
times |
Numeric vector of sample times in seconds, same length as
|
dt |
Positive sampling interval in seconds. Mutually exclusive with
|
start |
Start time in seconds used only when |
center |
Logical; if |
scale |
Character; |
mask |
Optional logical vector the same length as |
span |
Temporal window in seconds for the HRF, passed to |
Amplitude modulation on this sampling grid is the same linear model as
convolving the sampled series: if is the (possibly centered) feature
and is the convolution operator, the regressor is . Zeros are
kept, so an all-TR / all-sample series stays an all-sample series. No
unmodulated companion regressor is added.
Centering and scaling are applied to the feature before convolution.
For a whole-run series (mask = NULL),
Away from the run edges, is nearly constant (overlapping
HRFs sum to a plateau), so with a GLM intercept the centered and raw
columns are affinely equivalent. They differ by the HRF-length ramp of
at the start and end of the run. Default center = TRUE
removes that boundary term; center = FALSE keeps it. scale = "sd"
only changes the feature's units (still pre-convolution). It is not
standardization of the final BOLD-space column after filtering.
Use mask when the feature should be centered only during an on-period
(stimulus present, task on, and so on). Center and scale then use only the
on-samples, and off-mask samples stay 0:
That is not an affine transform of the all-sample series. Pair it with a separate boxcar for the on-period if you want presence and intensity as two questions.
Each sample is a zero-order-hold bin of width , not a
unit-mass impulse. Evaluate with precision less than or equal to the
feature sampling interval. Compared with
regressor(times, amplitude = values, duration = 0), the predicted BOLD
is smaller by about .
An S3 object of class c("FeatureReg", "Reg", "list"). Evaluation
uses the same convolution path as regressor().
# 10 Hz envelope over 8 seconds, demeaned dt <- 0.1 t <- seq(0, 8, by = dt) rms <- abs(sin(2 * pi * t / 4)) feat <- feature_regressor(rms, dt = dt, hrf = HRF_SPMG1) grid <- seq(0, 12, by = 1) y <- evaluate(feat, grid, precision = dt) # Same ZOH encoding as a duration-dt event regressor (no centering) feat_raw <- feature_regressor(rms, dt = dt, center = FALSE, scale = "none") ev <- regressor(t, HRF_SPMG1, duration = dt, amplitude = rms) all.equal(evaluate(feat_raw, grid, precision = dt), evaluate(ev, grid, precision = dt))# 10 Hz envelope over 8 seconds, demeaned dt <- 0.1 t <- seq(0, 8, by = dt) rms <- abs(sin(2 * pi * t / 4)) feat <- feature_regressor(rms, dt = dt, hrf = HRF_SPMG1) grid <- seq(0, 12, by = 1) y <- evaluate(feat, grid, precision = dt) # Same ZOH encoding as a duration-dt event regressor (no centering) feat_raw <- feature_regressor(rms, dt = dt, center = FALSE, scale = "none") ev <- regressor(t, HRF_SPMG1, duration = dt, amplitude = rms) all.equal(evaluate(feat_raw, grid, precision = dt), evaluate(ev, grid, precision = dt))
Command line entrypoint for fmrihrf
fmrihrf_cli(args = commandArgs(trailingOnly = TRUE))fmrihrf_cli(args = commandArgs(trailingOnly = TRUE))
args |
Character vector of command line arguments. |
Integer exit status. '0' indicates success, '1' indicates a domain failure, and '2' indicates usage or runtime errors.
'gen_hrf' takes a base HRF function or object and applies optional lag, blocking, and normalization decorators based on arguments.
gen_hrf( hrf, lag = 0, width = 0, precision = 0.1, half_life = Inf, summate = TRUE, normalize = FALSE, name = NULL, span = NULL, hrf_norm = "none", ... )gen_hrf( hrf, lag = 0, width = 0, precision = 0.1, half_life = Inf, summate = TRUE, normalize = FALSE, name = NULL, span = NULL, hrf_norm = "none", ... )
hrf |
A function 'f(t)' or an existing 'HRF' object. |
lag |
Optional lag in seconds. If non-zero, applies 'lag_hrf'. |
width |
Optional block width in seconds. If non-zero, applies 'block_hrf'. |
precision |
Sampling precision for block convolution (passed to 'block_hrf'). Default is 0.1. |
half_life |
Half-life decay parameter for exponential decay in seconds (passed to 'block_hrf'). Default is Inf (no decay). |
summate |
Passed to 'block_hrf()' when 'width > 0'. If 'TRUE' (default), block responses are integrated; if 'FALSE', the integrated response is scaled by total block weight so amplitude does not grow with block width. |
normalize |
If TRUE, applies 'normalise_hrf' at the end. Default is FALSE. |
name |
Optional name for the *final* HRF object. If NULL (default), a name is generated based on the base HRF and applied decorators. |
span |
Optional span for the *final* HRF object. If NULL (default), the span is determined by the base HRF and decorators. |
hrf_norm |
Fixed-scale normalization mode passed to [normalize_hrf()]. The default '"none"' preserves the raw HRF. 'normalize = TRUE' remains a back-compatible shortcut for '"unit_peak_per_basis"'; do not supply both. |
... |
Extra arguments passed to the *base* HRF function if 'hrf' is a function. |
A final 'HRF' object, potentially modified by decorators.
# Lagged SPMG1 grf_lag <- gen_hrf(HRF_SPMG1, lag=3) # Blocked Gaussian grf_block <- gen_hrf(hrf_gaussian, width=5, precision=0.2) # Lagged and Blocked, then Normalized grf_both_norm <- gen_hrf(HRF_SPMG1, lag=2, width=4, normalize=TRUE)# Lagged SPMG1 grf_lag <- gen_hrf(HRF_SPMG1, lag=3) # Blocked Gaussian grf_block <- gen_hrf(hrf_gaussian, width=5, precision=0.2) # Lagged and Blocked, then Normalized grf_both_norm <- gen_hrf(HRF_SPMG1, lag=2, width=4, normalize=TRUE)
The 'gen_hrf_blocked' function creates a blocked HRF by convolving the input HRF with a boxcar function. This can be used to model block designs in fMRI analysis.
gen_hrf_blocked( hrf = hrf_gaussian, width = 5, precision = 0.1, half_life = Inf, summate = TRUE, normalize = FALSE, ... ) hrf_blocked( hrf = hrf_gaussian, width = 5, precision = 0.1, half_life = Inf, summate = TRUE, normalize = FALSE, ... )gen_hrf_blocked( hrf = hrf_gaussian, width = 5, precision = 0.1, half_life = Inf, summate = TRUE, normalize = FALSE, ... ) hrf_blocked( hrf = hrf_gaussian, width = 5, precision = 0.1, half_life = Inf, summate = TRUE, normalize = FALSE, ... )
hrf |
A function representing the hemodynamic response function. Default is 'hrf_gaussian'. |
width |
A numeric value specifying the width of the block in seconds. Default is 5. |
precision |
A numeric value specifying the sampling resolution in seconds. Default is 0.1. |
half_life |
A numeric value specifying the half-life of the exponential decay function, used to model response attenuation. Default is 'Inf', which means no decay. |
summate |
Logical; if TRUE (default), responses accumulate (peak grows with duration). If FALSE, the convolution is averaged so the temporal profile is preserved but peak amplitude does not grow with duration. |
normalize |
A logical value indicating whether to rescale the output so that the peak of the output is 1. Default is 'FALSE'. |
... |
Extra arguments passed to the HRF function. |
A function representing the blocked HRF.
A function representing the blocked HRF.
hrf_blocked(): alias for gen_hrf_blocked
Other gen_hrf:
gen_hrf_lagged()
# Deprecated: use gen_hrf(..., width = 10) or block_hrf(HRF, width = 10)# Deprecated: use gen_hrf(..., width = 10) or block_hrf(HRF, width = 10)
The 'gen_hrf_lagged' function takes an HRF function and applies a specified lag to it. This can be useful for modeling time-delayed hemodynamic responses.
gen_hrf_lagged(hrf, lag = 2, normalize = FALSE, ...) hrf_lagged(hrf, lag = 2, normalize = FALSE, ...)gen_hrf_lagged(hrf, lag = 2, normalize = FALSE, ...) hrf_lagged(hrf, lag = 2, normalize = FALSE, ...)
hrf |
A function representing the underlying HRF to be shifted. |
lag |
A numeric value specifying the lag or delay in seconds to apply to the HRF. This can also be a vector of lags, in which case the function returns an HRF set. |
normalize |
A logical value indicating whether to rescale the output so that the maximum absolute value is 1. Defaults to 'FALSE'. |
... |
Extra arguments supplied to the 'hrf' function. |
A function representing the lagged HRF. If 'lag' is a vector of lags, the function returns an HRF set.
an lagged hrf function
hrf_lagged(): alias for gen_hrf_lagged
Other gen_hrf:
gen_hrf_blocked()
Other gen_hrf:
gen_hrf_blocked()
hrf_lag5 <- gen_hrf_lagged(HRF_SPMG1, lag=5) hrf_lag5(0:20)hrf_lag5 <- gen_hrf_lagged(HRF_SPMG1, lag=5) hrf_lag5(0:20)
Retrieves an HRF by name from the registry and optionally applies decorators. This provides a unified interface for creating both pre-defined HRF objects and custom basis sets with specified parameters.
getHRF( name = "spmg1", nbasis = 5, span = 24, lag = 0, width = 0, summate = TRUE, normalize = FALSE, hrf_norm = "none", ... )getHRF( name = "spmg1", nbasis = 5, span = 24, lag = 0, width = 0, summate = TRUE, normalize = FALSE, hrf_norm = "none", ... )
name |
Character string specifying the HRF type. Options include:
|
nbasis |
Number of basis functions (for basis set types) |
span |
Temporal window in seconds (default: 24) |
lag |
Time lag in seconds to apply (default: 0) |
width |
Block width for block designs (default: 0) |
summate |
Whether to sum responses in block designs (default: TRUE) |
normalize |
Whether to normalize the HRF (default: FALSE) |
hrf_norm |
Fixed-scale normalization mode passed to [normalize_hrf()]. The default '"none"' leaves the HRF unchanged. 'normalize = TRUE' is the legacy per-basis unit-peak shortcut; do not supply both. |
... |
Additional arguments passed to generator functions (e.g., |
For single HRF types (spmg1, gamma, gaussian), the function returns pre-defined objects. For basis set types (fir, bspline, fourier, daguerre), it calls the appropriate generator function with the specified parameters.
An HRF object
# Get pre-defined canonical HRF canonical <- getHRF("spmg1") # Create custom FIR basis with 20 bins fir20 <- getHRF("fir", nbasis = 20, span = 30) # Create B-spline basis with lag bs_lag <- getHRF("bspline", nbasis = 8, lag = 2) # Create blocked Gaussian HRF block_gauss <- getHRF("gaussian", width = 5)# Get pre-defined canonical HRF canonical <- getHRF("spmg1") # Create custom FIR basis with 20 bins fir20 <- getHRF("fir", nbasis = 20, span = 30) # Create B-spline basis with lag bs_lag <- getHRF("bspline", nbasis = 8, lag = 2) # Create blocked Gaussian HRF block_gauss <- getHRF("gaussian", width = 5)
Generic accessor for converting block-wise onsets to global onsets.
global_onsets(x, ...) ## S3 method for class 'sampling_frame' global_onsets(x, onsets, blockids, ...)global_onsets(x, ...) ## S3 method for class 'sampling_frame' global_onsets(x, onsets, blockids, ...)
x |
Object describing the sampling frame |
... |
Additional arguments passed to methods |
onsets |
Numeric vector of onset times within blocks |
blockids |
Integer vector identifying the block for each onset. Values must be whole numbers with no NAs. |
Numeric vector of global onset times
# Convert block-relative onsets to global timing sframe <- sampling_frame(blocklens = c(100, 120), TR = 2) global_onsets(sframe, onsets = c(10, 20), blockids = c(1, 2))# Convert block-relative onsets to global timing sframe <- sampling_frame(blocklens = c(100, 120), TR = 2) global_onsets(sframe, onsets = c(10, 20), blockids = c(1, 2))
The 'HRF' function creates an object representing a hemodynamic response function (HRF). It is a class constructor for HRFs.
HRF(fun, name, nbasis = 1, span = 24, param_names = NULL)HRF(fun, name, nbasis = 1, span = 24, param_names = NULL)
fun |
A function representing the hemodynamic response, mapping from time to BOLD response. |
name |
A string specifying the name of the function. |
nbasis |
An integer representing the number of basis functions, e.g., the columnar dimension of the HRF. Default is 1. |
span |
A numeric value representing the span in seconds of the HRF. Default is 24. |
param_names |
A character vector containing the names of the parameters for the HRF function. |
The package provides several pre-defined HRF types that can be used in modeling fMRI responses:
**Canonical HRFs:** * ‘"spmg1"' or 'HRF_SPMG1': SPM’s canonical HRF (single basis function) * '"spmg2"' or 'HRF_SPMG2': SPM canonical + temporal derivative (2 basis functions) * '"spmg3"' or 'HRF_SPMG3': SPM canonical + temporal and dispersion derivatives (3 basis functions) * '"gaussian"' or 'HRF_GAUSSIAN': Gaussian-shaped HRF with peak around 5-6s * '"gamma"' or 'HRF_GAMMA': Gamma function-based HRF with longer tail
**Flexible basis sets:** * '"bspline"' or '"bs"' or 'HRF_BSPLINE': B-spline basis for flexible HRF modeling * '"tent"': Tent (triangular) basis functions for flexible HRF modeling * '"daguerre"' or 'HRF_DAGUERRE': Daguerre basis functions
To see a complete list of available HRF types with details, use the 'list_available_hrfs()' function.
An HRF object with the specified properties.
hrf <- HRF(hrf_gamma, "gamma", nbasis=1, param_names=c("shape", "rate")) resp <- evaluate(hrf, seq(0, 24, by=1)) # List all available HRF types list_available_hrfs(details = TRUE)hrf <- HRF(hrf_gamma, "gamma", nbasis=1, param_names=c("shape", "rate")) resp <- evaluate(hrf, seq(0, 24, by=1)) # List all available HRF types list_available_hrfs(details = TRUE)
Constructs the basis set for the Lag-Width-Undershoot (LWU) HRF model,
intended for Taylor expansion-based fitting. The basis consists of the
LWU HRF evaluated at a given expansion point theta0, and its
partial derivatives with respect to its parameters (tau, sigma, rho).
hrf_basis_lwu(theta0, t, normalize_primary = "none")hrf_basis_lwu(theta0, t, normalize_primary = "none")
theta0 |
A numeric vector of length 3 specifying the expansion point
|
t |
A numeric vector of time points (in seconds) at which to evaluate the basis. |
normalize_primary |
Character string, one of |
A numeric matrix of dimension length(t) x 4. The columns represent:
h0: LWU HRF evaluated at theta0
d_tau: Partial derivative with respect to tau at theta0
d_sigma: Partial derivative with respect to sigma at theta0
d_rho: Partial derivative with respect to rho at theta0
Other hrf_functions:
hrf_boxcar(),
hrf_bspline(),
hrf_gamma(),
hrf_gaussian(),
hrf_inv_logit(),
hrf_lwu(),
hrf_mexhat(),
hrf_sine(),
hrf_spmg1(),
hrf_time(),
hrf_weighted()
t_points <- seq(0, 30, by = 0.5) theta0_default <- c(tau = 6, sigma = 1, rho = 0.35) # Generate the basis set lwu_basis <- hrf_basis_lwu(theta0_default, t_points) dim(lwu_basis) # Should be length(t_points) x 4 head(lwu_basis) # Plot the basis functions matplot(t_points, lwu_basis, type = "l", lty = 1, main = "LWU HRF Basis Functions", ylab = "Value", xlab = "Time (s)") legend("topright", colnames(lwu_basis), col = 1:4, lty = 1, cex = 0.8) # Example with primary HRF normalization (not typical for Taylor fitting step) lwu_basis_norm_h0 <- hrf_basis_lwu(theta0_default, t_points, normalize_primary = "height") plot(t_points, lwu_basis_norm_h0[,1], type="l", main="Normalized h0 in Basis") max(abs(lwu_basis_norm_h0[,1])) # Should be 1t_points <- seq(0, 30, by = 0.5) theta0_default <- c(tau = 6, sigma = 1, rho = 0.35) # Generate the basis set lwu_basis <- hrf_basis_lwu(theta0_default, t_points) dim(lwu_basis) # Should be length(t_points) x 4 head(lwu_basis) # Plot the basis functions matplot(t_points, lwu_basis, type = "l", lty = 1, main = "LWU HRF Basis Functions", ylab = "Value", xlab = "Time (s)") legend("topright", colnames(lwu_basis), col = 1:4, lty = 1, cex = 0.8) # Example with primary HRF normalization (not typical for Taylor fitting step) lwu_basis_norm_h0 <- hrf_basis_lwu(theta0_default, t_points, normalize_primary = "height") plot(t_points, lwu_basis_norm_h0[,1], type="l", main="Normalized h0 in Basis") max(abs(lwu_basis_norm_h0[,1])) # Should be 1
Creates a simple boxcar (step function) HRF that is constant within a time window starting at t=0 and zero outside. Unlike traditional HRFs, this has no hemodynamic delay - it represents an instantaneous response.
hrf_boxcar(width, amplitude = 1, normalize = FALSE)hrf_boxcar(width, amplitude = 1, normalize = FALSE)
width |
Duration of the boxcar window in seconds. |
amplitude |
Height of the boxcar (default: 1). |
normalize |
Logical; if |
In a least-squares GLM, an isolated event's coefficient for a boxcar of
amplitude 1 estimates the mean signal in the window. With
normalize = TRUE the boxcar has unit area (amplitude
1/width), so the coefficient estimates the integrated signal over
the window (the mean multiplied by width).
For delayed windows (not starting at t=0), use lag_hrf to shift
the boxcar in time.
An HRF object that can be used with regressor() and other
fmrihrf functions.
The width is fixed when the HRF is created. The duration
parameter in regressor() does not modify the boxcar
width—it controls how long the neural input is sustained (which then gets
convolved with this HRF). For trial-varying boxcar widths, use a list of HRFs:
widths <- c(4, 6, 8) hrfs <- lapply(widths, function(w) hrf_boxcar(width = w, normalize = TRUE)) reg <- regressor(onsets = c(0, 20, 40), hrf = hrfs)
hrf_weighted for weighted/shaped boxcars,
lag_hrf to shift the window in time
Other hrf_functions:
hrf_basis_lwu(),
hrf_bspline(),
hrf_gamma(),
hrf_gaussian(),
hrf_inv_logit(),
hrf_lwu(),
hrf_mexhat(),
hrf_sine(),
hrf_spmg1(),
hrf_time(),
hrf_weighted()
# Simple boxcar of 5 seconds width hrf1 <- hrf_boxcar(width = 5) t <- seq(-1, 10, by = 0.1) plot(t, evaluate(hrf1, t), type = "s", main = "Simple Boxcar HRF") # Unit-area boxcar: beta estimates the integrated signal over 0-5 s # (5 x the mean); with amplitude 1, beta estimates the mean itself hrf2 <- hrf_boxcar(width = 5, normalize = TRUE) # Use in a regressor with trial-varying widths hrf_short <- hrf_boxcar(width = 4, normalize = TRUE) hrf_long <- hrf_boxcar(width = 8, normalize = TRUE) reg <- regressor(onsets = c(0, 20), hrf = list(hrf_short, hrf_long)) # For delayed windows, use lag_hrf decorator hrf_delayed <- lag_hrf(hrf_boxcar(width = 5), lag = 10) # Window from 10-15s# Simple boxcar of 5 seconds width hrf1 <- hrf_boxcar(width = 5) t <- seq(-1, 10, by = 0.1) plot(t, evaluate(hrf1, t), type = "s", main = "Simple Boxcar HRF") # Unit-area boxcar: beta estimates the integrated signal over 0-5 s # (5 x the mean); with amplitude 1, beta estimates the mean itself hrf2 <- hrf_boxcar(width = 5, normalize = TRUE) # Use in a regressor with trial-varying widths hrf_short <- hrf_boxcar(width = 4, normalize = TRUE) hrf_long <- hrf_boxcar(width = 8, normalize = TRUE) reg <- regressor(onsets = c(0, 20), hrf = list(hrf_short, hrf_long)) # For delayed windows, use lag_hrf decorator hrf_delayed <- lag_hrf(hrf_boxcar(width = 5), lag = 10) # Window from 10-15s
The 'hrf_bspline' function computes a B-spline basis for an HRF at time
points 't'. The 'N' basis functions are the interior functions of a clamped
B-spline basis on [0, span] with evenly spaced knots: the functions
anchored at the two boundaries are dropped, so every basis function (and
therefore every fitted HRF) is zero at onset and at the end of the span.
Outside [0, span] the basis is zero.
hrf_bspline(t, span = 24, N = 5, degree = 3, ...)hrf_bspline(t, span = 24, N = 5, degree = 3, ...)
t |
A vector of time points. |
span |
A numeric value representing the temporal window over which the basis set spans. Default value is 24. |
N |
An integer representing the number of basis functions. Must be at
least |
degree |
An integer representing the degree of the spline. Default value is 3. |
... |
Further arguments passed to |
A matrix representing the B-spline basis for the HRF at the given time points 't'.
Other hrf_functions:
hrf_basis_lwu(),
hrf_boxcar(),
hrf_gamma(),
hrf_gaussian(),
hrf_inv_logit(),
hrf_lwu(),
hrf_mexhat(),
hrf_sine(),
hrf_spmg1(),
hrf_time(),
hrf_weighted()
# Compute the B-spline HRF representation for time points from 0 to 20 with 0.5 increments hrfb <- hrf_bspline(seq(0, 20, by = .5), N = 4, degree = 2)# Compute the B-spline HRF representation for time points from 0 to 20 with 0.5 increments hrfb <- hrf_bspline(seq(0, 20, by = .5), N = 4, degree = 2)
Generates an HRF object using B-spline basis functions with custom parameters.
This is the generator function that creates HRF objects with variable numbers
of basis functions, unlike the pre-defined HRF_BSPLINE which has 5 functions.
hrf_bspline_generator(nbasis = 5, span = 24)hrf_bspline_generator(nbasis = 5, span = 24)
nbasis |
Number of basis functions (default: 5) |
span |
Temporal window in seconds (default: 24) |
An HRF object of class c("BSpline_HRF", "HRF", "function")
HRF_objects for pre-defined HRF objects,
getHRF for a unified interface to create HRFs
# Create B-spline basis with 10 functions custom_bs <- hrf_bspline_generator(nbasis = 10) t <- seq(0, 24, by = 0.1) response <- evaluate(custom_bs, t) matplot(t, response, type = "l", main = "B-spline HRF with 10 basis functions")# Create B-spline basis with 10 functions custom_bs <- hrf_bspline_generator(nbasis = 10) t <- seq(0, 24, by = 0.1) response <- evaluate(custom_bs, t) matplot(t, response, type = "l", main = "B-spline HRF with 10 basis functions")
Generates an HRF object using Daguerre spherical basis functions with custom parameters. These are orthogonal polynomials that naturally decay to zero.
hrf_daguerre_generator(nbasis = 3, scale = 4)hrf_daguerre_generator(nbasis = 3, scale = 4)
nbasis |
Number of basis functions (default: 3) |
scale |
Scale parameter for the time axis (default: 4) |
Daguerre basis functions are orthogonal polynomials on [0,Inf) with respect to the weight function w(x) = x^2 * exp(-x). They are particularly useful for modeling hemodynamic responses as they naturally decay to zero and can capture various response shapes with few parameters.
An HRF object of class c("Daguerre_HRF", "HRF", "function")
HRF_objects for pre-defined HRF objects,
getHRF for a unified interface to create HRFs
# Create Daguerre basis with 5 functions custom_dag <- hrf_daguerre_generator(nbasis = 5, scale = 3) t <- seq(0, 24, by = 0.1) response <- evaluate(custom_dag, t) matplot(t, response, type = "l", main = "Daguerre HRF with 5 basis functions")# Create Daguerre basis with 5 functions custom_dag <- hrf_daguerre_generator(nbasis = 5, scale = 3) t <- seq(0, 24, by = 0.1) response <- evaluate(custom_dag, t) matplot(t, response, type = "l", main = "Daguerre HRF with 5 basis functions")
Generates an HRF object using Finite Impulse Response (FIR) basis functions with custom parameters. Each basis function represents a time bin with a value of 1 in that bin and 0 elsewhere.
hrf_fir_generator(nbasis = 12, span = 24)hrf_fir_generator(nbasis = 12, span = 24)
nbasis |
Number of time bins (default: 12) |
span |
Temporal window in seconds (default: 24) |
The FIR basis divides the time window into nbasis equal bins.
Each basis function is an indicator function for its corresponding bin.
This provides maximum flexibility but requires more parameters than
smoother basis sets like B-splines.
An HRF object of class c("FIR_HRF", "HRF", "function")
HRF_objects for pre-defined HRF objects,
getHRF for a unified interface to create HRFs,
hrf_bspline_generator for a smoother alternative
# Create FIR basis with 20 bins over 30 seconds custom_fir <- hrf_fir_generator(nbasis = 20, span = 30) t <- seq(0, 30, by = 0.1) response <- evaluate(custom_fir, t) matplot(t, response, type = "l", main = "FIR HRF with 20 time bins") # Compare to default FIR with 12 bins default_fir <- HRF_FIR response_default <- evaluate(default_fir, t[1:241]) # 24 seconds matplot(t[1:241], response_default, type = "l", main = "Default FIR HRF (12 bins over 24s)")# Create FIR basis with 20 bins over 30 seconds custom_fir <- hrf_fir_generator(nbasis = 20, span = 30) t <- seq(0, 30, by = 0.1) response <- evaluate(custom_fir, t) matplot(t, response, type = "l", main = "FIR HRF with 20 time bins") # Compare to default FIR with 12 bins default_fir <- HRF_FIR response_default <- evaluate(default_fir, t[1:241]) # 24 seconds matplot(t[1:241], response_default, type = "l", main = "Default FIR HRF (12 bins over 24s)")
Generates a set of Fourier basis functions (sine and cosine pairs) over a given span.
hrf_fourier(t, span = 24, nbasis = 5)hrf_fourier(t, span = 24, nbasis = 5)
t |
A vector of time points. |
span |
The temporal window over which the basis functions span (default: 24). |
nbasis |
The number of basis functions (default: 5). Should be even for full sine-cosine pairs. |
A matrix of Fourier basis functions with nbasis columns.
# Create Fourier basis with 5 functions t <- seq(0, 24, by = 0.5) basis <- hrf_fourier(t, span = 24, nbasis = 5) matplot(t, basis, type = "l", main = "Fourier Basis Functions")# Create Fourier basis with 5 functions t <- seq(0, 24, by = 0.5) basis <- hrf_fourier(t, span = 24, nbasis = 5) matplot(t, basis, type = "l", main = "Fourier Basis Functions")
Generates an HRF object using Fourier basis functions (sine and cosine pairs) with custom parameters.
hrf_fourier_generator(nbasis = 5, span = 24)hrf_fourier_generator(nbasis = 5, span = 24)
nbasis |
Number of basis functions (default: 5). Should be even for complete sine-cosine pairs. |
span |
Temporal window in seconds (default: 24) |
The Fourier basis uses alternating sine and cosine functions with increasing frequencies. This provides a smooth, periodic basis set that can capture oscillatory components in the HRF.
An HRF object of class c("Fourier_HRF", "HRF", "function")
HRF_objects for pre-defined HRF objects,
getHRF for a unified interface to create HRFs
# Create Fourier basis with 8 functions custom_fourier <- hrf_fourier_generator(nbasis = 8) t <- seq(0, 24, by = 0.1) response <- evaluate(custom_fourier, t) matplot(t, response, type = "l", main = "Fourier HRF with 8 basis functions")# Create Fourier basis with 8 functions custom_fourier <- hrf_fourier_generator(nbasis = 8) t <- seq(0, 24, by = 0.1) response <- evaluate(custom_fourier, t) matplot(t, response, type = "l", main = "Fourier HRF with 8 basis functions")
Create a new HRF by linearly weighting the basis functions of an existing HRF. Useful when coefficients have been estimated for an FIR/bspline/SPMG3 basis and one wants a single functional HRF.
hrf_from_coefficients(hrf, h, ...) ## S3 method for class 'HRF' hrf_from_coefficients(hrf, h, name = NULL, ...)hrf_from_coefficients(hrf, h, ...) ## S3 method for class 'HRF' hrf_from_coefficients(hrf, h, name = NULL, ...)
hrf |
An object of class 'HRF'. |
h |
Numeric vector of length 'nbasis(hrf)' giving the weights. |
... |
Reserved for future extensions. |
name |
Optional name for the resulting HRF. |
A new 'HRF' object with 'nbasis = 1'.
# Create a custom HRF from SPMG3 basis coefficients coeffs <- c(1, 0.2, -0.1) # Main response + slight temporal shift - dispersion custom_hrf <- hrf_from_coefficients(HRF_SPMG3, coeffs) # Evaluate the custom HRF t <- seq(0, 20, by = 0.1) response <- evaluate(custom_hrf, t) # Create from FIR basis fir_coeffs <- c(0, 0.2, 0.5, 1, 0.8, 0.4, 0.1, 0, 0, 0, 0, 0) custom_fir <- hrf_from_coefficients(HRF_FIR, fir_coeffs)# Create a custom HRF from SPMG3 basis coefficients coeffs <- c(1, 0.2, -0.1) # Main response + slight temporal shift - dispersion custom_hrf <- hrf_from_coefficients(HRF_SPMG3, coeffs) # Evaluate the custom HRF t <- seq(0, 20, by = 0.1) response <- evaluate(custom_hrf, t) # Create from FIR basis fir_coeffs <- c(0, 0.2, 0.5, 1, 0.8, 0.4, 0.1, 0, 0, 0, 0, 0) custom_fir <- hrf_from_coefficients(HRF_FIR, fir_coeffs)
The 'hrf_gamma' function computes the gamma density-based HRF (hemodynamic response function) at given time points 't'.
hrf_gamma(t, shape = 6, rate = 1)hrf_gamma(t, shape = 6, rate = 1)
t |
A vector of time points. |
shape |
A numeric value representing the shape parameter for the gamma probability density function. Default value is 6. |
rate |
A numeric value representing the rate parameter for the gamma probability density function. Default value is 1. |
A numeric vector representing the gamma HRF at the given time points 't'.
Other hrf_functions:
hrf_basis_lwu(),
hrf_boxcar(),
hrf_bspline(),
hrf_gaussian(),
hrf_inv_logit(),
hrf_lwu(),
hrf_mexhat(),
hrf_sine(),
hrf_spmg1(),
hrf_time(),
hrf_weighted()
# Compute the gamma HRF representation for time points from 0 to 20 with 0.5 increments hrf_gamma_vals <- hrf_gamma(seq(0, 20, by = .5), shape = 6, rate = 1)# Compute the gamma HRF representation for time points from 0 to 20 with 0.5 increments hrf_gamma_vals <- hrf_gamma(seq(0, 20, by = .5), shape = 6, rate = 1)
The 'hrf_gaussian' function computes the Gaussian density-based HRF (hemodynamic response function) at given time points 't'.
hrf_gaussian(t, mean = 6, sd = 2)hrf_gaussian(t, mean = 6, sd = 2)
t |
A vector of time points. |
mean |
A numeric value representing the mean of the Gaussian probability density function. Default value is 6. |
sd |
A numeric value representing the standard deviation of the Gaussian probability density function. Default value is 2. |
A numeric vector representing the Gaussian HRF at the given time points 't'.
Other hrf_functions:
hrf_basis_lwu(),
hrf_boxcar(),
hrf_bspline(),
hrf_gamma(),
hrf_inv_logit(),
hrf_lwu(),
hrf_mexhat(),
hrf_sine(),
hrf_spmg1(),
hrf_time(),
hrf_weighted()
# Compute the Gaussian HRF representation for time points from 0 to 20 with 0.5 increments hrf_gaussian_vals <- hrf_gaussian(seq(0, 20, by = .5), mean = 6, sd = 2)# Compute the Gaussian HRF representation for time points from 0 to 20 with 0.5 increments hrf_gaussian_vals <- hrf_gaussian(seq(0, 20, by = .5), mean = 6, sd = 2)
Segments: 0->f1 (h1), f1->1 (h2), 1->f2 (h3), f2->0 (h4). Negative f1 gives an initial dip; negative f2 gives an undershoot. Peak is at t = h1 + h2 (amplitude 1 by construction).
hrf_half_cosine(t, h1 = 1, h2 = 5, h3 = 7, h4 = 7, f1 = 0, f2 = 0)hrf_half_cosine(t, h1 = 1, h2 = 5, h3 = 7, h4 = 7, f1 = 0, f2 = 0)
t |
Numeric vector of times (s) |
h1, h2, h3, h4
|
Segment durations (s). Must be > 0. |
f1 |
Initial dip level (default 0), typically in [-0.2, 0] |
f2 |
Undershoot level (default 0), typically in [-0.3, 0] |
Numeric vector same length as t
t <- seq(0, 30, by = 0.1) y <- hrf_half_cosine(t)t <- seq(0, 30, by = 0.1) y <- hrf_half_cosine(t)
A hemodynamic response function using the difference of two Inverse Logit functions.
hrf_inv_logit(t, mu1 = 6, s1 = 1, mu2 = 16, s2 = 1, lag = 0)hrf_inv_logit(t, mu1 = 6, s1 = 1, mu2 = 16, s2 = 1, lag = 0)
t |
A vector of times. |
mu1 |
The time-to-peak for the rising phase (mean of the first logistic function). |
s1 |
The width (slope) of the first logistic function. |
mu2 |
The time-to-peak for the falling phase (mean of the second logistic function). |
s2 |
The width (slope) of the second logistic function. |
lag |
The time delay (default: 0). |
A vector of the difference of two Inverse Logit HRF values.
Other hrf_functions:
hrf_basis_lwu(),
hrf_boxcar(),
hrf_bspline(),
hrf_gamma(),
hrf_gaussian(),
hrf_lwu(),
hrf_mexhat(),
hrf_sine(),
hrf_spmg1(),
hrf_time(),
hrf_weighted()
hrf_inv_logit_basis <- hrf_inv_logit(seq(0, 20, by = 0.5), mu1 = 6, s1 = 1, mu2 = 16, s2 = 1)hrf_inv_logit_basis <- hrf_inv_logit(seq(0, 20, by = 0.5), mu1 = 6, s1 = 1, mu2 = 16, s2 = 1)
Computes the Lag-Width-Undershoot (LWU) hemodynamic response function. This model uses two Gaussian components to model the main response and an optional undershoot.
hrf_lwu(t, tau = 6, sigma = 2.5, rho = 0.35, normalize = "none")hrf_lwu(t, tau = 6, sigma = 2.5, rho = 0.35, normalize = "none")
t |
A numeric vector of time points (in seconds). |
tau |
Lag of the main Gaussian component (time-to-peak of the positive lobe, in seconds). Default: 6. |
sigma |
Width (standard deviation) of the main Gaussian component (in seconds). Must be > 0.05. Default: 2.5. |
rho |
Amplitude of the undershoot Gaussian component, relative to the main component. Must be between 0 and 1.5. Default: 0.35. |
normalize |
Character string specifying normalization type. Either "none" for no normalization (default) or "height" to scale the HRF so its maximum absolute value is 1. |
The LWU model formula combines a positive Gaussian peak and a negative undershoot: h(t; tau, sigma, rho) = exp(-(t-tau)^2/(2*sigma^2)) - rho * exp(-(t-tau-2*sigma)^2/(2*(1.6*sigma)^2))
A numeric vector representing the LWU HRF values at the given time points 't'.
Other hrf_functions:
hrf_basis_lwu(),
hrf_boxcar(),
hrf_bspline(),
hrf_gamma(),
hrf_gaussian(),
hrf_inv_logit(),
hrf_mexhat(),
hrf_sine(),
hrf_spmg1(),
hrf_time(),
hrf_weighted()
t_points <- seq(0, 30, by = 0.1) # Default LWU HRF lwu_default <- hrf_lwu(t_points) plot(t_points, lwu_default, type = "l", main = "LWU HRF (Default Params)", ylab = "Amplitude") # LWU HRF with no undershoot lwu_no_undershoot <- hrf_lwu(t_points, rho = 0) lines(t_points, lwu_no_undershoot, col = "blue") # LWU HRF with a wider main peak and larger undershoot lwu_custom <- hrf_lwu(t_points, tau = 7, sigma = 1.5, rho = 0.5) lines(t_points, lwu_custom, col = "red") legend("topright", c("Default", "No Undershoot (rho=0)", "Custom (tau=7, sigma=1.5, rho=0.5)"), col = c("black", "blue", "red"), lty = 1, cex = 0.8) # Height-normalized HRF lwu_normalized <- hrf_lwu(t_points, tau = 6, sigma = 1, rho = 0.35, normalize = "height") plot(t_points, lwu_normalized, type = "l", main = "Height-Normalized LWU HRF", ylab = "Amplitude") abline(h = c(-1, 1), lty = 2, col = "grey") # Max absolute value should be 1t_points <- seq(0, 30, by = 0.1) # Default LWU HRF lwu_default <- hrf_lwu(t_points) plot(t_points, lwu_default, type = "l", main = "LWU HRF (Default Params)", ylab = "Amplitude") # LWU HRF with no undershoot lwu_no_undershoot <- hrf_lwu(t_points, rho = 0) lines(t_points, lwu_no_undershoot, col = "blue") # LWU HRF with a wider main peak and larger undershoot lwu_custom <- hrf_lwu(t_points, tau = 7, sigma = 1.5, rho = 0.5) lines(t_points, lwu_custom, col = "red") legend("topright", c("Default", "No Undershoot (rho=0)", "Custom (tau=7, sigma=1.5, rho=0.5)"), col = c("black", "blue", "red"), lty = 1, cex = 0.8) # Height-normalized HRF lwu_normalized <- hrf_lwu(t_points, tau = 6, sigma = 1, rho = 0.35, normalize = "height") plot(t_points, lwu_normalized, type = "l", main = "Height-Normalized LWU HRF", ylab = "Amplitude") abline(h = c(-1, 1), lty = 2, col = "grey") # Max absolute value should be 1
The 'hrf_mexhat' function computes the Mexican hat wavelet-based HRF (hemodynamic response function) at given time points 't'.
hrf_mexhat(t, mean = 6, sd = 2)hrf_mexhat(t, mean = 6, sd = 2)
t |
A vector of time points. |
mean |
A numeric value representing the mean of the Mexican hat wavelet. Default value is 6. |
sd |
A numeric value representing the standard deviation of the Mexican hat wavelet. Default value is 2. |
A numeric vector representing the Mexican hat wavelet-based HRF at the given time points 't'.
Other hrf_functions:
hrf_basis_lwu(),
hrf_boxcar(),
hrf_bspline(),
hrf_gamma(),
hrf_gaussian(),
hrf_inv_logit(),
hrf_lwu(),
hrf_sine(),
hrf_spmg1(),
hrf_time(),
hrf_weighted()
# Compute the Mexican hat HRF representation for time points from 0 to 20 with 0.5 increments hrf_mexhat_vals <- hrf_mexhat(seq(0, 20, by = .5), mean = 6, sd = 2)# Compute the Mexican hat HRF representation for time points from 0 to 20 with 0.5 increments hrf_mexhat_vals <- hrf_mexhat(seq(0, 20, by = .5), mean = 6, sd = 2)
A collection of pre-defined HRF objects for common fMRI analysis scenarios. These objects can be used directly in model specifications or as templates for creating custom HRFs.
HRF_GAMMA(t) HRF_GAUSSIAN(t) HRF_SPMG1(t) HRF_SPMG2(t) HRF_SPMG3(t) HRF_BSPLINE(t) HRF_FIR(t)HRF_GAMMA(t) HRF_GAUSSIAN(t) HRF_SPMG1(t) HRF_SPMG2(t) HRF_SPMG3(t) HRF_BSPLINE(t) HRF_FIR(t)
t |
Numeric vector of time points (in seconds) at which to evaluate the HRF |
SPMG bases are raw continuous kernels. SPMG2 adds the analytic
time derivative; SPMG3 also adds the response-dispersion difference
(h(d) - h(d + 0.01))/0.01, holding the response mean and mass fixed.
Columns are not implicitly orthogonalized or independently rescaled.
Use normalize_hrf() or the existing normalization flags to select scaling.
These conventions do not reproduce every sampled SPM/Nilearn design detail.
When called as functions, return numeric vectors or matrices of HRF values.
When used as objects, they are HRF objects with class c("HRF", "function").
HRF_SPMG1SPM canonical HRF (single basis function)
HRF_SPMG2SPM canonical HRF with temporal derivative (2 basis functions)
HRF_SPMG3SPM canonical HRF with temporal and dispersion derivatives (3 basis functions)
HRF_GAMMAGamma function-based HRF
HRF_GAUSSIANGaussian function-based HRF
HRF_BSPLINEB-spline basis HRF (5 basis functions)
HRF_FIRFinite Impulse Response (FIR) basis HRF (12 basis functions)
The pre-defined objects above have fixed numbers of basis functions. To create basis sets with custom parameters (e.g., different numbers of basis functions), use one of these approaches:
Using getHRF():
getHRF("fir", nbasis = 20) - FIR basis with 20 functions
getHRF("bspline", nbasis = 10, span = 30) - B-spline with 10 functions
getHRF("fourier", nbasis = 7) - Fourier basis with 7 functions
getHRF("daguerre", nbasis = 5, scale = 3) - Daguerre basis
Using generator functions directly:
hrf_fir_generator(nbasis = 20, span = 30)
hrf_bspline_generator(nbasis = 10, span = 30)
hrf_fourier_generator(nbasis = 7, span = 24)
hrf_daguerre_generator(nbasis = 5, scale = 3)
All HRF objects can be:
Called as functions with time argument: HRF_SPMG1(t)
Used in model specifications: hrf(condition, basis = HRF_SPMG1)
Evaluated with evaluate() method
Combined with decorators like lag_hrf() or block_hrf()
evaluate.HRF for evaluating HRF objects,
gen_hrf for creating HRFs with decorators,
list_available_hrfs for listing all HRF types,
getHRF for creating HRFs by name with custom parameters,
hrf_fir_generator, hrf_bspline_generator,
hrf_fourier_generator, hrf_daguerre_generator
for creating custom basis sets directly
Other hrf:
deriv(),
penalty_matrix()
# Evaluate HRFs at specific time points times <- seq(0, 20, by = 0.5) # Single basis canonical HRF canonical_response <- HRF_SPMG1(times) plot(times, canonical_response, type = "l", main = "SPM Canonical HRF") # Multi-basis HRF with derivatives multi_response <- HRF_SPMG3(times) # Returns 3-column matrix matplot(times, multi_response, type = "l", main = "SPM HRF with Derivatives") # Gamma and Gaussian HRFs gamma_response <- HRF_GAMMA(times) gaussian_response <- HRF_GAUSSIAN(times) # Compare different HRF shapes plot(times, canonical_response, type = "l", col = "blue", main = "HRF Comparison", ylab = "Response") lines(times, gamma_response, col = "red") lines(times, gaussian_response, col = "green") legend("topright", c("SPM Canonical", "Gamma", "Gaussian"), col = c("blue", "red", "green"), lty = 1) # Create custom FIR basis with 20 bins custom_fir <- getHRF("fir", nbasis = 20, span = 30) fir_response <- evaluate(custom_fir, times) matplot(times, fir_response, type = "l", main = "Custom FIR with 20 bins") # Create custom B-spline basis custom_bspline <- hrf_bspline_generator(nbasis = 8, span = 25) bspline_response <- evaluate(custom_bspline, times) matplot(times, bspline_response, type = "l", main = "Custom B-spline with 8 basis functions")# Evaluate HRFs at specific time points times <- seq(0, 20, by = 0.5) # Single basis canonical HRF canonical_response <- HRF_SPMG1(times) plot(times, canonical_response, type = "l", main = "SPM Canonical HRF") # Multi-basis HRF with derivatives multi_response <- HRF_SPMG3(times) # Returns 3-column matrix matplot(times, multi_response, type = "l", main = "SPM HRF with Derivatives") # Gamma and Gaussian HRFs gamma_response <- HRF_GAMMA(times) gaussian_response <- HRF_GAUSSIAN(times) # Compare different HRF shapes plot(times, canonical_response, type = "l", col = "blue", main = "HRF Comparison", ylab = "Response") lines(times, gamma_response, col = "red") lines(times, gaussian_response, col = "green") legend("topright", c("SPM Canonical", "Gamma", "Gaussian"), col = c("blue", "red", "green"), lty = 1) # Create custom FIR basis with 20 bins custom_fir <- getHRF("fir", nbasis = 20, span = 30) fir_response <- evaluate(custom_fir, times) matplot(times, fir_response, type = "l", main = "Custom FIR with 20 bins") # Create custom B-spline basis custom_bspline <- hrf_bspline_generator(nbasis = 8, span = 25) bspline_response <- evaluate(custom_bspline, times) matplot(times, bspline_response, type = "l", main = "Custom B-spline with 8 basis functions")
hrf_palette() returns the colours used by every plotting function in
fmrihrf, so figures made with plot_hrfs(), plot_regressors(), and the
plot() methods share one visual language.
hrf_palette(n = NULL, type = c("categorical", "ordered"))hrf_palette(n = NULL, type = c("categorical", "ordered"))
n |
Number of colours. If |
type |
Either |
Two palette types are provided:
"categorical": up to six distinct hues for unordered series (conditions,
HRF families). Every colour has at least 3:1 contrast against both white
and near-black backgrounds, and the first four stay distinguishable under
simulated deuteranopia, protanopia, and tritanopia. When more than six
colours are requested the ordered palette is used instead.
"ordered": a hue ramp (violet, blue, teal, green, ochre) for series
with a natural order, such as the functions of a basis set, lags, or
event durations. Lightness is held in the same mid range, so the ramp is
readable on light and dark backgrounds. Brick red, the first
categorical colour, is left out of the ramp so that it keeps one meaning.
A character vector of hex colours.
scale_colour_hrf() for the matching ggplot2 scales.
hrf_palette() hrf_palette(3) hrf_palette(8, type = "ordered")hrf_palette() hrf_palette(3) hrf_palette(8, type = "ordered")
A hemodynamic response function using the Sine Basis Set.
hrf_sine(t, span = 24, N = 5)hrf_sine(t, span = 24, N = 5)
t |
A vector of times. |
span |
The temporal window over which the basis sets span (default: 24). |
N |
The number of basis functions (default: 5). |
A matrix of sine basis functions.
Other hrf_functions:
hrf_basis_lwu(),
hrf_boxcar(),
hrf_bspline(),
hrf_gamma(),
hrf_gaussian(),
hrf_inv_logit(),
hrf_lwu(),
hrf_mexhat(),
hrf_spmg1(),
hrf_time(),
hrf_weighted()
hrf_sine_basis <- hrf_sine(seq(0, 20, by = 0.5), N = 4)hrf_sine_basis <- hrf_sine(seq(0, 20, by = 0.5), N = 4)
A hemodynamic response function based on the SPM canonical double gamma parameterization.
hrf_spmg1(t, P1 = 5, P2 = 15, A1 = 1/120)hrf_spmg1(t, P1 = 5, P2 = 15, A1 = 1/120)
t |
A vector of time points. |
P1 |
The first exponent parameter (default: 5). |
P2 |
The second exponent parameter (default: 15). |
A1 |
Amplitude scaling factor for the positive gamma function component; normally fixed at 1/120 |
This function models the hemodynamic response using the canonical double gamma parameterization in the SPM software. The HRF is defined by a linear combination of two gamma functions with different exponents (P1 and P2) and amplitudes (A1 and A2). It is commonly used in fMRI data analysis to estimate the BOLD (blood-oxygen-level-dependent) signal changes associated with neural activity.
A vector of HRF values at the given time points.
Other hrf_functions:
hrf_basis_lwu(),
hrf_boxcar(),
hrf_bspline(),
hrf_gamma(),
hrf_gaussian(),
hrf_inv_logit(),
hrf_lwu(),
hrf_mexhat(),
hrf_sine(),
hrf_time(),
hrf_weighted()
# Generate a time vector time_points <- seq(0, 30, by=0.1) # Compute the HRF values using the SPM canonical double gamma parameterization hrf_values <- hrf_spmg1(time_points) # Plot the HRF values plot(time_points, hrf_values, type='l', main='SPM Canonical Double Gamma HRF')# Generate a time vector time_points <- seq(0, 30, by=0.1) # Compute the HRF values using the SPM canonical double gamma parameterization hrf_values <- hrf_spmg1(time_points) # Plot the HRF values plot(time_points, hrf_values, type='l', main='SPM Canonical Double Gamma HRF')
Generates an HRF object using tent (piecewise linear) basis functions with
custom parameters. This generator mirrors HRF_TENT but allows callers
to control the number of basis elements and temporal span.
hrf_tent_generator(nbasis = 5, span = 24)hrf_tent_generator(nbasis = 5, span = 24)
nbasis |
Number of tent basis functions (default: 5) |
span |
Temporal window in seconds (default: 24) |
An HRF object of class c("Tent_HRF", "HRF", "function")
HRF_objects for pre-defined HRF objects,
getHRF for a unified interface to create HRFs,
hrf_bspline_generator for a smoother alternative
# Create a tent basis with 6 functions over a 20 second window custom_tent <- hrf_tent_generator(nbasis = 6, span = 20) t <- seq(0, 20, by = 0.1) response <- evaluate(custom_tent, t) matplot(t, response, type = "l", main = "Tent HRF with 6 basis functions")# Create a tent basis with 6 functions over a 20 second window custom_tent <- hrf_tent_generator(nbasis = 6, span = 20) t <- seq(0, 20, by = 0.1) response <- evaluate(custom_tent, t) matplot(t, response, type = "l", main = "Tent HRF with 6 basis functions")
The 'hrf_time' function computes the value of an HRF, which is a simple linear function of time 't', when 't' is greater than 0 and less than 'maxt'.
hrf_time(t, maxt = 22)hrf_time(t, maxt = 22)
t |
A numeric value representing time in seconds. |
maxt |
A numeric value representing the maximum time point in the domain. Default value is 22. |
A numeric value representing the value of the HRF at the given time 't'.
Other hrf_functions:
hrf_basis_lwu(),
hrf_boxcar(),
hrf_bspline(),
hrf_gamma(),
hrf_gaussian(),
hrf_inv_logit(),
hrf_lwu(),
hrf_mexhat(),
hrf_sine(),
hrf_spmg1(),
hrf_weighted()
# Compute the HRF value for t = 5 seconds with the default maximum time hrf_val <- hrf_time(5) # Compute the HRF value for t = 5 seconds with a custom maximum time of 30 seconds hrf_val_custom_maxt <- hrf_time(5, maxt = 30)# Compute the HRF value for t = 5 seconds with the default maximum time hrf_val <- hrf_time(5) # Compute the HRF value for t = 5 seconds with a custom maximum time of 30 seconds hrf_val_custom_maxt <- hrf_time(5, maxt = 30)
Create a Toeplitz matrix for hemodynamic response function (HRF) convolution.
hrf_toeplitz(hrf, time, len, sparse = FALSE)hrf_toeplitz(hrf, time, len, sparse = FALSE)
hrf |
The hemodynamic response function. |
time |
A numeric vector representing the time points. |
len |
The length of the output Toeplitz matrix. |
sparse |
Logical, if TRUE, the output Toeplitz matrix is returned as a sparse matrix (default: FALSE). |
A Toeplitz matrix for HRF convolution.
# Create HRF and time points hrf_fun <- function(t) hrf_spmg1(t) times <- seq(0, 30, by = 1) # Create Toeplitz matrix H <- hrf_toeplitz(hrf_fun, times, len = 50) # Create sparse version H_sparse <- hrf_toeplitz(hrf_fun, times, len = 50, sparse = TRUE)# Create HRF and time points hrf_fun <- function(t) hrf_spmg1(t) times <- seq(0, 30, by = 1) # Create Toeplitz matrix H <- hrf_toeplitz(hrf_fun, times, len = 50) # Create sparse version H_sparse <- hrf_toeplitz(hrf_fun, times, len = 50, sparse = TRUE)
Creates a flexible weighted HRF starting at t=0 with user-specified weights. Unlike traditional HRFs, this has no built-in hemodynamic delay - it directly maps weights to time points, allowing for arbitrary temporal response shapes.
hrf_weighted( weights, width = NULL, times = NULL, method = c("constant", "linear"), normalize = FALSE )hrf_weighted( weights, width = NULL, times = NULL, method = c("constant", "linear"), normalize = FALSE )
weights |
Numeric vector of weights. Required. |
width |
Total duration of the window in seconds. If provided without
|
times |
Numeric vector of time points (in seconds, relative to t=0) where
weights are specified. Must be strictly increasing and start at 0 for
consistency with other HRFs. If provided, |
method |
Interpolation method between time points:
|
normalize |
Logical; if |
This is useful for summarising the signal in chosen post-stimulus windows
with a chosen temporal profile. In a least-squares GLM, an isolated event's
coefficient is the amplitude of the weight profile that best
matches the data, . This equals the
mean signal in the window only when all non-zero weights are 1 (a boxcar);
in general it is not a weighted mean. Normalizing rescales the coefficient
but does not change this.
There are two ways to specify the temporal structure:
width + weights: the window [0, width) is divided
into length(weights) equal bins (constant method) or the weights
are placed at evenly spaced points from 0 to width (linear method)
times + weights: explicit time points for each weight
(relative to t=0)
For delayed windows (not starting at t=0), use lag_hrf to shift
the weighted HRF in time.
An HRF object that can be used with regressor() and other
fmrihrf functions.
The temporal structure (width or times) is fixed when the HRF
is created. The duration parameter in regressor() does
not modify the weighted HRF's structure—it controls how long the
neural input is sustained (which then gets convolved with this HRF). For
trial-varying weighted HRFs, use a list of HRFs:
hrf_early <- hrf_weighted(width = 6, weights = c(1, 1, 0, 0), normalize = TRUE) hrf_late <- hrf_weighted(width = 6, weights = c(0, 0, 1, 1), normalize = TRUE) reg <- regressor(onsets = c(0, 20), hrf = list(hrf_early, hrf_late))
hrf_boxcar for simple uniform boxcars,
lag_hrf to shift the window in time,
empirical_hrf for HRFs from measured data
Other hrf_functions:
hrf_basis_lwu(),
hrf_boxcar(),
hrf_bspline(),
hrf_gamma(),
hrf_gaussian(),
hrf_inv_logit(),
hrf_lwu(),
hrf_mexhat(),
hrf_sine(),
hrf_spmg1(),
hrf_time()
# Simple: 6 s window split into 4 bins of 1.5 s (starting at 0, 1.5, 3, 4.5 s) hrf1 <- hrf_weighted(width = 6, weights = c(0.2, 0.5, 0.8, 0.3)) t <- seq(-1, 10, by = 0.1) plot(t, evaluate(hrf1, t), type = "s", main = "Weighted HRF (width + weights)") # Explicit times for precise control hrf2 <- hrf_weighted( times = c(0, 1, 3, 5, 6), weights = c(0.1, 0.5, 0.8, 0.5, 0.1), method = "linear" ) plot(t, evaluate(hrf2, t), type = "l", main = "Smooth Weighted HRF") # Normalized weights: all four bin weights sum to 1 hrf3 <- hrf_weighted( width = 8, weights = c(1, 2, 2, 1), normalize = TRUE ) # Trial-varying weighted HRFs hrf_early <- hrf_weighted(width = 6, weights = c(1, 1, 0, 0), normalize = TRUE) hrf_late <- hrf_weighted(width = 6, weights = c(0, 0, 1, 1), normalize = TRUE) reg <- regressor(onsets = c(0, 20), hrf = list(hrf_early, hrf_late)) # For delayed windows, use lag_hrf hrf_delayed <- lag_hrf(hrf_weighted(width = 5, weights = c(1, 2, 1)), lag = 10)# Simple: 6 s window split into 4 bins of 1.5 s (starting at 0, 1.5, 3, 4.5 s) hrf1 <- hrf_weighted(width = 6, weights = c(0.2, 0.5, 0.8, 0.3)) t <- seq(-1, 10, by = 0.1) plot(t, evaluate(hrf1, t), type = "s", main = "Weighted HRF (width + weights)") # Explicit times for precise control hrf2 <- hrf_weighted( times = c(0, 1, 3, 5, 6), weights = c(0.1, 0.5, 0.8, 0.5, 0.1), method = "linear" ) plot(t, evaluate(hrf2, t), type = "l", main = "Smooth Weighted HRF") # Normalized weights: all four bin weights sum to 1 hrf3 <- hrf_weighted( width = 8, weights = c(1, 2, 2, 1), normalize = TRUE ) # Trial-varying weighted HRFs hrf_early <- hrf_weighted(width = 6, weights = c(1, 1, 0, 0), normalize = TRUE) hrf_late <- hrf_weighted(width = 6, weights = c(0, 0, 1, 1), normalize = TRUE) reg <- regressor(onsets = c(0, 20), hrf = list(hrf_early, hrf_late)) # For delayed windows, use lag_hrf hrf_delayed <- lag_hrf(hrf_weighted(width = 5, weights = c(1, 2, 1)), lag = 10)
Install fmrihrf command line wrappers
install_cli(dest_dir = "~/.local/bin", overwrite = FALSE, commands = NULL)install_cli(dest_dir = "~/.local/bin", overwrite = FALSE, commands = NULL)
dest_dir |
Directory where wrapper commands should be copied. |
overwrite |
Logical; overwrite an existing command only when 'TRUE'. |
commands |
Optional character vector of command names to install. |
Invisibly, a named character vector of installed command paths.
Creates a new HRF object by applying a temporal lag to an existing HRF object.
lag_hrf(hrf, lag)lag_hrf(hrf, lag)
hrf |
The HRF object (of class 'HRF') to lag. |
lag |
The time lag in seconds to apply. Positive values shift the response later in time. |
A new HRF object representing the lagged function.
Other HRF_decorator_functions:
block_hrf(),
normalise_hrf(),
normalize_hrf()
lagged_spmg1 <- lag_hrf(HRF_SPMG1, 5) # Evaluate at time 10; equivalent to HRF_SPMG1(10 - 5) lagged_spmg1(10) HRF_SPMG1(5)lagged_spmg1 <- lag_hrf(HRF_SPMG1, 5) # Evaluate at time 10; equivalent to HRF_SPMG1(10 - 5) lagged_spmg1(10) HRF_SPMG1(5)
Reads the internal HRF registry to list available HRF types.
list_available_hrfs(details = FALSE)list_available_hrfs(details = FALSE)
details |
Logical; if TRUE, attempt to add descriptions (basic for now). |
A data frame with columns: name, type (object/generator), nbasis_default.
# List all available HRFs hrfs <- list_available_hrfs() print(hrfs) # List with details hrfs_detailed <- list_available_hrfs(details = TRUE) print(hrfs_detailed)# List all available HRFs hrfs <- list_available_hrfs() print(hrfs) # List with details hrfs_detailed <- list_available_hrfs(details = TRUE) print(hrfs_detailed)
'make_hrf' resolves a basis specification to an 'HRF' object and applies an optional temporal lag. The basis may be given as the name of a built-in HRF, as a generating function, or as an existing 'HRF' object.
make_hrf(basis, lag, nbasis = 1)make_hrf(basis, lag, nbasis = 1)
basis |
Character name of a built-in HRF, a function that generates HRF values, or an object of class 'HRF'. |
lag |
Numeric scalar giving the shift in seconds applied to the HRF. |
nbasis |
Integer specifying the number of basis functions when 'basis' is provided as a name. |
An object of class 'HRF' representing the lagged basis.
# Canonical SPM HRF delayed by 2 seconds h <- make_hrf("spmg1", lag = 2) h(0:5)# Canonical SPM HRF delayed by 2 seconds h <- make_hrf("spmg1", lag = 2) h(0:5)
Return the number of basis functions represented by an object.
nbasis(x, ...) ## S3 method for class 'HRF' nbasis(x, ...) ## S3 method for class 'Reg' nbasis(x, ...)nbasis(x, ...) ## S3 method for class 'HRF' nbasis(x, ...) ## S3 method for class 'Reg' nbasis(x, ...)
x |
Object containing HRF or regressor information. |
... |
Additional arguments passed to methods. |
This information is typically used when constructing penalty matrices or understanding the complexity of an HRF model or regressor.
Integer scalar giving the number of basis functions.
# Number of basis functions for different HRF types nbasis(HRF_SPMG1) # 1 basis function nbasis(HRF_SPMG3) # 3 basis functions (canonical + 2 derivatives) nbasis(HRF_BSPLINE) # 5 basis functions (default) # For a regressor reg <- regressor(onsets = c(10, 30, 50), hrf = HRF_SPMG3) nbasis(reg) # 3 (inherits from the HRF)# Number of basis functions for different HRF types nbasis(HRF_SPMG1) # 1 basis function nbasis(HRF_SPMG3) # 3 basis functions (canonical + 2 derivatives) nbasis(HRF_BSPLINE) # 5 basis functions (default) # For a regressor reg <- regressor(onsets = c(10, 30, 50), hrf = HRF_SPMG3) nbasis(reg) # 3 (inherits from the HRF)
Converts event timing information into a neural input function representing the underlying neural activity before HRF convolution. This function is useful for:
neural_input(x, ...) ## S3 method for class 'Reg' neural_input(x, start = 0, end = NULL, resolution = 0.33, ...)neural_input(x, ...) ## S3 method for class 'Reg' neural_input(x, start = 0, end = NULL, resolution = 0.33, ...)
x |
A regressor object containing event timing information |
... |
Additional arguments passed to methods |
start |
Numeric; start time of the input function |
end |
Numeric; end time of the input function |
resolution |
Numeric; temporal resolution in seconds (default: 0.33) |
Creating stimulus functions for fMRI analysis
Modeling sustained vs. transient neural activity
Generating inputs for HRF convolution
Visualizing the temporal structure of experimental designs
A list containing:
Numeric vector of time points
Numeric vector of input amplitudes at each time point
regressor, evaluate.Reg, HRF_SPMG1
# Create a regressor with multiple events reg <- regressor( onsets = c(10, 30, 50), duration = c(2, 2, 2), amplitude = c(1, 1.5, 0.8), hrf = HRF_SPMG1 ) # Generate neural input function input <- neural_input(reg, start = 0, end = 60, resolution = 0.5) # Plot the neural input function plot(input$time, input$neural_input, type = "l", xlab = "Time (s)", ylab = "Neural Input", main = "Neural Input Function") # Create regressor with varying durations reg_sustained <- regressor( onsets = c(10, 30), duration = c(5, 10), # sustained activity amplitude = c(1, 1), hrf = HRF_SPMG1 ) # Generate and compare neural inputs input_sustained <- neural_input( reg_sustained, start = 0, end = 60, resolution = 0.5 )# Create a regressor with multiple events reg <- regressor( onsets = c(10, 30, 50), duration = c(2, 2, 2), amplitude = c(1, 1.5, 0.8), hrf = HRF_SPMG1 ) # Generate neural input function input <- neural_input(reg, start = 0, end = 60, resolution = 0.5) # Plot the neural input function plot(input$time, input$neural_input, type = "l", xlab = "Time (s)", ylab = "Neural Input", main = "Neural Input Function") # Create regressor with varying durations reg_sustained <- regressor( onsets = c(10, 30), duration = c(5, 10), # sustained activity amplitude = c(1, 1), hrf = HRF_SPMG1 ) # Generate and compare neural inputs input_sustained <- neural_input( reg_sustained, start = 0, end = 60, resolution = 0.5 )
Back-compatible spelling and behavior for independently peak-normalizing every basis column. New code can use 'normalize_hrf(hrf, "unit_peak_per_basis")' explicitly.
normalise_hrf(hrf)normalise_hrf(hrf)
hrf |
An object of class 'HRF'. |
A unit-peak 'HRF' object.
Other HRF_decorator_functions:
block_hrf(),
lag_hrf(),
normalize_hrf()
gauss_unnorm <- as_hrf(function(t) 5 * dnorm(t, 6, 2), name = "unnorm_gauss") gauss_norm <- normalise_hrf(gauss_unnorm) max(gauss_norm(seq(0, 20, by = 0.1)))gauss_unnorm <- as_hrf(function(t) 5 * dnorm(t, 6, 2), name = "unnorm_gauss") gauss_norm <- normalise_hrf(gauss_unnorm) max(gauss_norm(seq(0, 20, by = 0.1)))
Creates an HRF whose evaluations are divided by constants computed once on a fixed reference grid. The resulting scale therefore does not depend on the time points supplied in later calls.
normalize_hrf(hrf, mode)normalize_hrf(hrf, mode)
hrf |
An object of class 'HRF'. |
mode |
Normalization convention: '"spm"' divides by the sum of the canonical basis on the 1,600-point 0–32 second Nilearn reference grid; '"unit_peak"' divides every basis by the canonical basis peak; '"unit_integral"' divides every basis by the canonical basis trapezoidal integral; '"unit_peak_per_basis"' scales each basis independently; and '"none"' returns 'hrf' unchanged. |
For multi-basis HRFs, '"spm"', '"unit_peak"', and '"unit_integral"' use one scalar computed from the first (canonical) column and apply it uniformly. This preserves the relative scale of derivative bases. Only '"unit_peak_per_basis"' rescales columns independently.
An 'HRF' object with fixed normalization.
Other HRF_decorator_functions:
block_hrf(),
lag_hrf(),
normalise_hrf()
spm_scaled <- normalize_hrf(HRF_SPMG1, "spm") reference_grid <- seq(0, 32, length.out = 1600) sum(spm_scaled(reference_grid)) peak_scaled <- normalize_hrf(HRF_SPMG2, "unit_peak") max(abs(peak_scaled(seq(0, 24, by = 0.01))[, 1]))spm_scaled <- normalize_hrf(HRF_SPMG1, "spm") reference_grid <- seq(0, 32, length.out = 1600) sum(spm_scaled(reference_grid)) peak_scaled <- normalize_hrf(HRF_SPMG2, "unit_peak") max(abs(peak_scaled(seq(0, 24, by = 0.01))[, 1]))
Generic accessor returning event onset times in seconds.
onsets(x, ...) ## S3 method for class 'Reg' onsets(x, ...)onsets(x, ...) ## S3 method for class 'Reg' onsets(x, ...)
x |
Object containing onset information |
... |
Additional arguments passed to methods |
Numeric vector of onsets
# Create a regressor with event onsets reg <- regressor(onsets = c(1, 5, 10, 15), hrf = HRF_SPMG1, span = 20) onsets(reg)# Create a regressor with event onsets reg <- regressor(onsets = c(1, 5, 10, 15), hrf = HRF_SPMG1, span = 20) onsets(reg)
Generate a penalty matrix for regularizing HRF basis coefficients. The penalty matrix encodes shape priors that discourage implausible or overly wiggly HRF estimates. Different HRF types use different penalty structures:
FIR/B-spline/Tent bases: Roughness penalties based on discrete derivatives
SPM canonical + derivatives: Differential shrinkage of derivative terms
Fourier bases: Penalties on high-frequency components
Daguerre bases: Increasing weights on higher-order terms
Default: Identity matrix (ridge penalty)
penalty_matrix(x, ...) ## S3 method for class 'HRF' penalty_matrix(x, order = 2, ...) ## S3 method for class 'BSpline_HRF' penalty_matrix(x, order = 2, ...) ## S3 method for class 'Tent_HRF' penalty_matrix(x, order = 2, ...) ## S3 method for class 'FIR_HRF' penalty_matrix(x, order = 2, ...) ## S3 method for class 'SPMG2_HRF' penalty_matrix(x, order = 2, shrink_deriv = 2, ...) ## S3 method for class 'SPMG3_HRF' penalty_matrix(x, order = 2, shrink_deriv = 2, ...) ## S3 method for class 'Fourier_HRF' penalty_matrix(x, order = 2, ...) ## S3 method for class 'Daguerre_HRF' penalty_matrix(x, order = 2, ...)penalty_matrix(x, ...) ## S3 method for class 'HRF' penalty_matrix(x, order = 2, ...) ## S3 method for class 'BSpline_HRF' penalty_matrix(x, order = 2, ...) ## S3 method for class 'Tent_HRF' penalty_matrix(x, order = 2, ...) ## S3 method for class 'FIR_HRF' penalty_matrix(x, order = 2, ...) ## S3 method for class 'SPMG2_HRF' penalty_matrix(x, order = 2, shrink_deriv = 2, ...) ## S3 method for class 'SPMG3_HRF' penalty_matrix(x, order = 2, shrink_deriv = 2, ...) ## S3 method for class 'Fourier_HRF' penalty_matrix(x, order = 2, ...) ## S3 method for class 'Daguerre_HRF' penalty_matrix(x, order = 2, ...)
x |
The HRF object or basis specification |
... |
Additional arguments passed to specific methods |
order |
Integer specifying the order of the penalty (default: 2) |
shrink_deriv |
Numeric; penalty weight for derivative terms in SPMG2/SPMG3 bases (default: 2) |
The penalty matrix R is used in regularized estimation as lambda * h^T R h, where h are the basis coefficients and lambda is the regularization parameter. Well-designed penalty matrices can significantly improve HRF estimation by encoding smoothness or other shape constraints.
A symmetric positive definite penalty matrix of dimension nbasis(x) x nbasis(x)
[nbasis()], [HRF_objects]
Other hrf:
HRF_objects,
deriv()
# FIR basis with smoothness penalty fir_hrf <- HRF_FIR R_fir <- penalty_matrix(fir_hrf) # B-spline basis with second-order smoothness bspline_hrf <- HRF_BSPLINE R_bspline <- penalty_matrix(bspline_hrf, order = 2) # SPM canonical with derivative shrinkage spmg3_hrf <- HRF_SPMG3 R_spmg3 <- penalty_matrix(spmg3_hrf, shrink_deriv = 4)# FIR basis with smoothness penalty fir_hrf <- HRF_FIR R_fir <- penalty_matrix(fir_hrf) # B-spline basis with second-order smoothness bspline_hrf <- HRF_BSPLINE R_bspline <- penalty_matrix(bspline_hrf, order = 2) # SPM canonical with derivative shrinkage spmg3_hrf <- HRF_SPMG3 R_spmg3 <- penalty_matrix(spmg3_hrf, shrink_deriv = 4)
Plots one or more HRF objects on shared axes. Multi-basis HRFs are
expanded into one curve per basis function, so plot_hrfs(HRF_SPMG3)
shows the canonical response and both derivatives. Uses ggplot2 when
available, otherwise base graphics. Colours come from hrf_palette().
plot_hrfs( ..., time = NULL, normalize = FALSE, labels = NULL, title = NULL, subtitle = NULL, use_ggplot = TRUE, draw = TRUE, basis = c("all", "first"), layout = c("overlay", "stack"), palette = c("auto", "categorical", "ordered"), reference = NULL, reference_label = NULL, scales = c("free_y", "fixed") )plot_hrfs( ..., time = NULL, normalize = FALSE, labels = NULL, title = NULL, subtitle = NULL, use_ggplot = TRUE, draw = TRUE, basis = c("all", "first"), layout = c("overlay", "stack"), palette = c("auto", "categorical", "ordered"), reference = NULL, reference_label = NULL, scales = c("free_y", "fixed") )
... |
HRF objects to compare. Can be passed as individual arguments or as a named list. |
time |
Numeric vector of time points. If NULL (default), uses seq(0, max_span, by = 0.1) where max_span is the maximum span across all HRFs. |
normalize |
Logical; if TRUE, normalize all HRFs to peak at 1. Useful for comparing shapes regardless of amplitude. Default is FALSE. |
labels |
Character vector of labels, either one per HRF or one per
plotted curve (after basis expansion). If NULL (default), uses the 'name'
attribute of each HRF; basis functions of a single basis set are labelled
|
title |
Character string for the plot title. If NULL (default), uses the HRF name for a single HRF and "HRF comparison" otherwise. |
subtitle |
Character string for the plot subtitle. If NULL (default), no subtitle is shown. |
use_ggplot |
Logical; if TRUE and ggplot2 is available, use ggplot2 for plotting. If FALSE, use base R graphics. Default is TRUE. |
draw |
Logical; draw the plot (default TRUE). With |
basis |
Either |
layout |
Either |
palette |
Colour palette: |
reference |
Optional HRF drawn as a dashed grey reference curve (for example the canonical HRF), outside the colour legend. It is evaluated on the same time grid and normalized like the other HRFs. |
reference_label |
Text naming the reference curve, shown as a caption. Defaults to the reference HRF's name. |
scales |
For |
Inside a knitr document the result is returned visibly and printed by knitr, which lets document themes (for example dark-mode figure twins) handle the ggplot. At the console the plot is drawn immediately and the data are returned invisibly.
A data frame in long format with columns 'time', 'HRF', and
'response'. With ggplot2, the plot is stored in the "plot" attribute.
The value is invisible except when drawing inside knitr, where it prints
like a ggplot object: it is drawn when it is the visible value of a chunk
(not when assigned), and subsetting or modifying it, or calling
as.data.frame(), returns a plain data frame. Call as.data.frame()
before combining results with rbind() or dplyr verbs.
# Compare canonical HRFs plot_hrfs(HRF_SPMG1, HRF_GAMMA, HRF_GAUSSIAN) # A basis set: one curve per basis function plot_hrfs(HRF_SPMG3, labels = c("Canonical", "Temporal derivative", "Dispersion derivative")) # HRFs ordered by a parameter use the ordered palette plot_hrfs(block_hrf(HRF_SPMG1, width = 1), block_hrf(HRF_SPMG1, width = 3), block_hrf(HRF_SPMG1, width = 5), labels = c("1 s", "3 s", "5 s"), palette = "ordered", title = "Effect of event duration") # Normalize for shape comparison plot_hrfs(HRF_SPMG1, HRF_GAMMA, HRF_GAUSSIAN, normalize = TRUE, subtitle = "All HRFs normalized to peak at 1") # Use base R graphics instead of ggplot2 plot_hrfs(HRF_SPMG1, HRF_GAMMA, use_ggplot = FALSE)# Compare canonical HRFs plot_hrfs(HRF_SPMG1, HRF_GAMMA, HRF_GAUSSIAN) # A basis set: one curve per basis function plot_hrfs(HRF_SPMG3, labels = c("Canonical", "Temporal derivative", "Dispersion derivative")) # HRFs ordered by a parameter use the ordered palette plot_hrfs(block_hrf(HRF_SPMG1, width = 1), block_hrf(HRF_SPMG1, width = 3), block_hrf(HRF_SPMG1, width = 5), labels = c("1 s", "3 s", "5 s"), palette = "ordered", title = "Effect of event duration") # Normalize for shape comparison plot_hrfs(HRF_SPMG1, HRF_GAMMA, HRF_GAUSSIAN, normalize = TRUE, subtitle = "All HRFs normalized to peak at 1") # Use base R graphics instead of ggplot2 plot_hrfs(HRF_SPMG1, HRF_GAMMA, use_ggplot = FALSE)
Plots one or more regressors on a shared time axis. Regressors built from a
basis set are expanded into one curve per basis function, and a
regressor_set() is expanded into one curve per condition. Event onsets
(and durations) are drawn as bars under the curves; zero-duration events get
a minimum bar width of 0.4% of the time range so that they stay visible. Uses ggplot2 when
available, otherwise base graphics. Colours come from hrf_palette().
plot_regressors( ..., grid = NULL, labels = NULL, title = NULL, subtitle = NULL, show_onsets = NULL, onset_alpha = 0.8, precision = NULL, use_ggplot = TRUE, draw = TRUE, basis = c("all", "first"), layout = c("overlay", "stack"), samples = NULL, palette = c("auto", "categorical", "ordered"), scales = c("free_y", "fixed") )plot_regressors( ..., grid = NULL, labels = NULL, title = NULL, subtitle = NULL, show_onsets = NULL, onset_alpha = 0.8, precision = NULL, use_ggplot = TRUE, draw = TRUE, basis = c("all", "first"), layout = c("overlay", "stack"), samples = NULL, palette = c("auto", "categorical", "ordered"), scales = c("free_y", "fixed") )
... |
Regressor objects ( |
grid |
Numeric vector of time points for evaluation. If NULL (default),
a 0.25 s grid covering all regressors is used. Use a fine grid for the
curve and |
labels |
Character vector of labels, one per regressor or one per plotted curve. If NULL (default), uses list names, condition levels, or "Regressor_1", "Regressor_2", etc. |
title |
Character string for the plot title. If NULL (default), uses "Regressor comparison". |
subtitle |
Character string for the plot subtitle. If NULL, no subtitle. |
show_onsets |
Logical or character. If TRUE, mark events for every
regressor in its own colour. If "first", mark only the first regressor's
events, in grey. If FALSE, hide event marks. If NULL (default), uses
"first" for |
onset_alpha |
Alpha transparency for event marks. Default is 0.8. |
precision |
Numeric sampling precision for HRF evaluation. If NULL (default), the grid spacing capped at 0.33 s. |
use_ggplot |
Logical; if TRUE and ggplot2 is available, use ggplot2 for plotting. If FALSE, use base R graphics. Default is TRUE. |
draw |
Logical; draw the plot (default TRUE). With |
basis |
Either |
layout |
Either |
samples |
Optional numeric vector of sample times (for example scan acquisition times). The regressors are evaluated there and drawn as points on the curves. |
palette |
Colour palette: |
scales |
For |
Inside a knitr document the result is returned visibly and printed by knitr, which lets document themes (for example dark-mode figure twins) handle the ggplot. At the console the plot is drawn immediately and the data are returned invisibly.
A data frame in long format with columns 'time', 'Regressor', and
'response'. With ggplot2, the plot is stored in the "plot" attribute.
The value is invisible except when drawing inside knitr, where it prints
like a ggplot object (see plot_hrfs()).
# Create regressors with different HRFs onsets <- c(10, 30, 50) reg1 <- regressor(onsets, HRF_SPMG1) reg2 <- regressor(onsets, HRF_GAMMA) reg3 <- regressor(onsets, HRF_GAUSSIAN) # Compare regressors plot_regressors(reg1, reg2, reg3, labels = c("SPM Canonical", "Gamma", "Gaussian")) # Show the scan-time samples of a regressor (TR = 2 s) plot_regressors(reg1, samples = seq(0, 80, by = 2), labels = "SPMG1") # One panel per basis function of a basis-set regressor plot_regressors(regressor(c(10, 40), HRF_SPMG3), layout = "stack") # Compare original vs shifted regressor reg_shifted <- shift(reg1, 5) plot_regressors(reg1, reg_shifted, labels = c("Original", "Shifted +5s"), show_onsets = TRUE)# Create regressors with different HRFs onsets <- c(10, 30, 50) reg1 <- regressor(onsets, HRF_SPMG1) reg2 <- regressor(onsets, HRF_GAMMA) reg3 <- regressor(onsets, HRF_GAUSSIAN) # Compare regressors plot_regressors(reg1, reg2, reg3, labels = c("SPM Canonical", "Gamma", "Gaussian")) # Show the scan-time samples of a regressor (TR = 2 s) plot_regressors(reg1, samples = seq(0, 80, by = 2), labels = "SPMG1") # One panel per basis function of a basis-set regressor plot_regressors(regressor(c(10, 40), HRF_SPMG3), layout = "stack") # Compare original vs shifted regressor reg_shifted <- shift(reg1, 5) plot_regressors(reg1, reg_shifted, labels = c("Original", "Shifted +5s"), show_onsets = TRUE)
Plot a Feature Regressor
## S3 method for class 'FeatureReg' plot( x, grid = NULL, show_onsets = FALSE, onset_color = NULL, onset_alpha = 0.5, precision = NULL, layout = c("stack", "overlay"), ... )## S3 method for class 'FeatureReg' plot( x, grid = NULL, show_onsets = FALSE, onset_color = NULL, onset_alpha = 0.5, precision = NULL, layout = c("stack", "overlay"), ... )
x |
A |
grid |
Numeric vector of time points for evaluation. If |
show_onsets |
Logical; if |
onset_color |
Colour for sample-time ticks. If NULL (default), a neutral grey. |
onset_alpha |
Alpha transparency for sample-time ticks. Default is 0.5. |
precision |
Numeric sampling precision for HRF evaluation. If NULL (default), the grid spacing capped at 0.33 s. |
layout |
For multi-basis HRFs, |
... |
Additional arguments passed to the underlying plot functions. |
Invisibly returns a data frame with the time and response values.
feat <- feature_regressor(abs(sin(seq(0, 8, by = 0.1))), dt = 0.1) plot(feat, grid = seq(0, 12, by = 0.5))feat <- feature_regressor(abs(sin(seq(0, 8, by = 0.1))), dt = 0.1) plot(feat, grid = seq(0, 12, by = 0.5))
Draws an HRF with base graphics. Single-basis HRFs show the response curve
with its peak annotated. Multi-basis HRFs (e.g., HRF_SPMG3) show every
basis function, coloured along the ordered hrf_palette().
## S3 method for class 'HRF' plot(x, time = NULL, normalize = FALSE, show_peak = TRUE, ...)## S3 method for class 'HRF' plot(x, time = NULL, normalize = FALSE, show_peak = TRUE, ...)
x |
An HRF object |
time |
Numeric vector of time points. If NULL (default), uses seq(0, span, by = 0.1) where span is the HRF's span attribute. |
normalize |
Logical; if TRUE, normalize responses to peak at 1. Default is FALSE. |
show_peak |
Logical; if TRUE (default for single-basis HRFs), annotate the peak time on the plot. |
... |
Additional arguments passed to |
Invisibly returns a data frame with the time and response values (useful for further customization).
plot_hrfs() for ggplot2 comparisons of several HRFs.
# Plot single-basis HRF plot(HRF_SPMG1) # Plot multi-basis HRF plot(HRF_SPMG3) # Plot with normalization plot(HRF_GAMMA, normalize = TRUE) # Custom time range plot(HRF_SPMG1, time = seq(0, 30, by = 0.5))# Plot single-basis HRF plot(HRF_SPMG1) # Plot multi-basis HRF plot(HRF_SPMG3) # Plot with normalization plot(HRF_GAMMA, normalize = TRUE) # Custom time range plot(HRF_SPMG1, time = seq(0, 30, by = 0.5))
Draws the predicted BOLD time course of a regressor with base graphics. Event onsets are marked with ticks along the time axis. Regressors built from a basis set are drawn one basis function per panel by default.
## S3 method for class 'Reg' plot( x, grid = NULL, show_onsets = TRUE, onset_color = NULL, onset_alpha = 0.5, precision = NULL, layout = c("stack", "overlay"), ... )## S3 method for class 'Reg' plot( x, grid = NULL, show_onsets = TRUE, onset_color = NULL, onset_alpha = 0.5, precision = NULL, layout = c("stack", "overlay"), ... )
x |
A |
grid |
Numeric vector of time points for evaluation. If NULL (default), uses a grid from 0 to max(onsets) + span with step 0.25 s. |
show_onsets |
Logical; if TRUE (default), mark event onsets with ticks on the time axis. |
onset_color |
Colour for onset ticks. If NULL (default), a neutral grey. |
onset_alpha |
Alpha transparency for onset ticks. Default is 0.5. |
precision |
Numeric sampling precision for HRF evaluation. If NULL (default), the grid spacing capped at 0.33 s, so sharp HRF edges are drawn where they occur. |
layout |
For multi-basis regressors, |
... |
Additional arguments passed to |
Invisibly returns a data frame with the time and response values.
plot_regressors() for ggplot2 comparisons of several regressors.
# Create and plot a simple regressor reg <- regressor(onsets = c(10, 30, 50), hrf = HRF_SPMG1) plot(reg) # Plot with custom time grid plot(reg, grid = seq(0, 80, by = 1)) # Plot without onset markers plot(reg, show_onsets = FALSE) # A basis-set regressor: one panel per basis function plot(regressor(c(10, 40), HRF_SPMG3))# Create and plot a simple regressor reg <- regressor(onsets = c(10, 30, 50), hrf = HRF_SPMG1) plot(reg) # Plot with custom time grid plot(reg, grid = seq(0, 80, by = 1)) # Plot without onset markers plot(reg, show_onsets = FALSE) # A basis-set regressor: one panel per basis function plot(regressor(c(10, 40), HRF_SPMG3))
Displays a concise summary of an HRF object including its name, number of basis functions, temporal span, and parameters (if any).
## S3 method for class 'HRF' print(x, ...)## S3 method for class 'HRF' print(x, ...)
x |
An HRF object |
... |
Additional arguments (unused) |
Invisibly returns the HRF object
# Print canonical HRF print(HRF_SPMG1) # Print multi-basis HRF print(HRF_SPMG3) # Print Gaussian HRF print(HRF_GAUSSIAN)# Print canonical HRF print(HRF_SPMG1) # Print multi-basis HRF print(HRF_SPMG3) # Print Gaussian HRF print(HRF_GAUSSIAN)
Provides a concise summary of the regressor object using the cli package.
## S3 method for class 'Reg' print(x, ...) ## S3 method for class 'sampling_frame' print(x, ...)## S3 method for class 'Reg' print(x, ...) ## S3 method for class 'sampling_frame' print(x, ...)
x |
A 'Reg' object. |
... |
Not used. |
No return value, called for side effects (prints to console)
r <- regressor(onsets = c(1, 10, 20), hrf = HRF_SPMG1, duration = 0, amplitude = 1, span = 40) print(r)r <- regressor(onsets = c(1, 10, 20), hrf = HRF_SPMG1, duration = 0, amplitude = 1, span = 40) print(r)
Create a new HRF by linearly weighting the basis functions of an existing HRF. This is useful for turning estimated basis coefficients into a single functional HRF.
S3 method for 'HRF' objects that returns a matrix mapping basis coefficients to sampled HRF values at the provided time grid. For single-basis HRFs, this returns a one-column matrix. For multi-basis HRFs (e.g., SPMG2/SPMG3, FIR, B-spline), this returns a matrix with one column per basis function.
reconstruction_matrix(hrf, sframe, ...) ## S3 method for class 'HRF' reconstruction_matrix(hrf, sframe, ...)reconstruction_matrix(hrf, sframe, ...) ## S3 method for class 'HRF' reconstruction_matrix(hrf, sframe, ...)
hrf |
An object of class 'HRF'. |
sframe |
A numeric vector of times, or a 'sampling_frame' object from which times are extracted via 'samples()'. |
... |
Additional arguments passed to 'samples()' when 'sframe' is a 'sampling_frame', and to 'evaluate()' for HRF evaluation. |
Reconstruction matrix for an HRF basis
Returns a matrix that converts basis coefficients into a
sampled HRF shape.
A numeric matrix with one column per basis function.
A numeric matrix of dimension 'length(times) x nbasis(hrf)'.
# Create reconstruction matrix for basis functions hrf <- HRF_SPMG2 # 2-basis HRF times <- seq(0, 20, by = 0.5) rmat <- reconstruction_matrix(hrf, times) dim(rmat) # Shows dimensions# Create reconstruction matrix for basis functions hrf <- HRF_SPMG2 # 2-basis HRF times <- seq(0, 20, by = 0.5) rmat <- reconstruction_matrix(hrf, times) dim(rmat) # Shows dimensions
Creates an object representing event-related regressors for fMRI modeling. This function defines event onsets and associates them with a hemodynamic response function (HRF) to generate predicted time courses.
regressor( onsets, hrf = HRF_SPMG1, duration = 0, amplitude = 1, span = 40, summate = TRUE, drop_zero_amplitude = TRUE )regressor( onsets, hrf = HRF_SPMG1, duration = 0, amplitude = 1, span = 40, summate = TRUE, drop_zero_amplitude = TRUE )
onsets |
A numeric vector of event onset times in seconds. |
hrf |
The hemodynamic response function (HRF) to convolve with the events. This can be:
Defaults to 'HRF_SPMG1'. |
duration |
A numeric scalar or vector specifying the duration of each event in seconds. If scalar, it's applied to all events. Defaults to 0 (impulse events). |
amplitude |
A numeric scalar or vector specifying the amplitude (scaling factor) for each event. If scalar, it's applied to all events. Defaults to 1. |
span |
The temporal window (in seconds) over which the HRF is defined or evaluated. This influences the length of the convolution. If not provided, it may be inferred from the 'hrf' object or default to 40s. For list HRFs, the maximum span across all HRFs is used. **Note:** Unlike some previous versions, the 'span' is not automatically adjusted based on 'duration'; ensure the provided or inferred 'span' is sufficient for your longest event duration. |
summate |
Logical scalar; if 'TRUE' (default), the HRF response amplitude scales with the duration of sustained events (via weighted integration). If 'FALSE', weighted integration is normalized by total block weight so amplitude does not grow with duration. |
drop_zero_amplitude |
Logical scalar; if 'TRUE' (default), remove zero-amplitude events. Set to 'FALSE' to retain all supplied onsets, including zero-valued samples in a continuous feature series. |
This function serves as the main public interface for creating regressor objects. Internally, it utilizes the 'Reg()' constructor which performs validation and efficient storage. The resulting object can be evaluated at specific time points using the 'evaluate()' function.
Events with an amplitude of 0 are filtered out when 'drop_zero_amplitude = TRUE'.
## Trial-Varying HRFs
When 'hrf' is a list of HRF objects, each event can have its own HRF. This is useful for trial-wise analyses where different events may have different temporal characteristics (e.g., different boxcar windows, different weights). The list must have either length 1 (recycled to all events) or length equal to the number of onsets.
An S3 object of class 'Reg' and 'list' containing processed event information and the HRF specification. The object includes a 'filtered_all' attribute indicating whether all events were removed due to zero or 'NA' amplitudes, and an 'hrf_is_list' attribute indicating whether trial-varying HRFs are used.
# Create a simple regressor with 3 events reg <- regressor(onsets = c(10, 30, 50), hrf = HRF_SPMG1) # Regressor with durations and amplitudes reg2 <- regressor( onsets = c(10, 30, 50), duration = c(2, 2, 2), amplitude = c(1, 1.5, 0.8), hrf = HRF_SPMG1 ) # Using different HRF types reg_gamma <- regressor(onsets = c(10, 30), hrf = "gamma") # Evaluate regressor at specific time points times <- seq(0, 60, by = 0.1) response <- evaluate(reg, times) # Trial-varying HRFs: different boxcar windows for each event hrf1 <- hrf_boxcar(width = 4, normalize = TRUE) hrf2 <- hrf_boxcar(width = 6, normalize = TRUE) reg_varying <- regressor( onsets = c(10, 30), hrf = list(hrf1, hrf2) )# Create a simple regressor with 3 events reg <- regressor(onsets = c(10, 30, 50), hrf = HRF_SPMG1) # Regressor with durations and amplitudes reg2 <- regressor( onsets = c(10, 30, 50), duration = c(2, 2, 2), amplitude = c(1, 1.5, 0.8), hrf = HRF_SPMG1 ) # Using different HRF types reg_gamma <- regressor(onsets = c(10, 30), hrf = "gamma") # Evaluate regressor at specific time points times <- seq(0, 60, by = 0.1) response <- evaluate(reg, times) # Trial-varying HRFs: different boxcar windows for each event hrf1 <- hrf_boxcar(width = 4, normalize = TRUE) hrf2 <- hrf_boxcar(width = 6, normalize = TRUE) reg_varying <- regressor( onsets = c(10, 30), hrf = list(hrf1, hrf2) )
'regressor_design' extends [regressor_set()] by allowing onsets to be specified relative to individual blocks and by directly returning the evaluated design matrix.
regressor_design( onsets, fac, block, sframe, hrf = HRF_SPMG1, duration = 0, amplitude = 1, span = 40, precision = 0.33, method = c("conv", "loop", "fft", "Rconv"), sparse = FALSE, summate = TRUE )regressor_design( onsets, fac, block, sframe, hrf = HRF_SPMG1, duration = 0, amplitude = 1, span = 40, precision = 0.33, method = c("conv", "loop", "fft", "Rconv"), sparse = FALSE, summate = TRUE )
onsets |
Numeric vector of event onset times, expressed relative to the start of their corresponding block. |
fac |
A factor (or object coercible to a factor) indicating the condition for each onset. |
block |
Integer vector identifying the block for each onset. Values must be valid block indices for 'sframe'. |
sframe |
A [sampling_frame] describing the temporal structure of the experiment. |
hrf |
Hemodynamic response function shared by all conditions. |
duration |
Numeric scalar or vector of event durations. |
amplitude |
Numeric scalar or vector of event amplitudes. |
span |
Numeric scalar giving the HRF span in seconds. |
precision |
Numeric precision used during convolution. |
method |
Evaluation method passed to [evaluate()]. |
sparse |
Logical; if 'TRUE' a sparse design matrix is returned. |
summate |
Logical; passed to [regressor()]. |
A numeric matrix (or sparse matrix) with one column per factor level and one row per sample defined by 'sframe'.
# Create a sampling frame for 2 blocks, 100 scans each, TR=2 sframe <- sampling_frame(blocklens = c(100, 100), TR = 2) # Events in block-relative time onsets <- c(10, 30, 50, 20, 40, 60) conditions <- factor(c("A", "B", "A", "B", "A", "B")) blocks <- c(1, 1, 1, 2, 2, 2) # Build design matrix design <- regressor_design( onsets = onsets, fac = conditions, block = blocks, sframe = sframe, hrf = HRF_SPMG1 ) # Design matrix has 200 rows (total scans) and 2 columns (conditions) dim(design)# Create a sampling frame for 2 blocks, 100 scans each, TR=2 sframe <- sampling_frame(blocklens = c(100, 100), TR = 2) # Events in block-relative time onsets <- c(10, 30, 50, 20, 40, 60) conditions <- factor(c("A", "B", "A", "B", "A", "B")) blocks <- c(1, 1, 1, 2, 2, 2) # Build design matrix design <- regressor_design( onsets = onsets, fac = conditions, block = blocks, sframe = sframe, hrf = HRF_SPMG1 ) # Design matrix has 200 rows (total scans) and 2 columns (conditions) dim(design)
Creates a set of regressors, one for each level of a factor. Each condition shares the same HRF and other parameters but has distinct onsets, durations and amplitudes.
regressor_set( onsets, fac, hrf = HRF_SPMG1, duration = 0, amplitude = 1, span = 40, summate = TRUE ) ## S3 method for class 'RegSet' evaluate( x, grid, precision = 0.33, method = c("conv", "loop", "fft", "Rconv"), sparse = FALSE, ... )regressor_set( onsets, fac, hrf = HRF_SPMG1, duration = 0, amplitude = 1, span = 40, summate = TRUE ) ## S3 method for class 'RegSet' evaluate( x, grid, precision = 0.33, method = c("conv", "loop", "fft", "Rconv"), sparse = FALSE, ... )
onsets |
Numeric vector of event onset times. |
fac |
A factor (or object coercible to a factor) indicating the condition for each onset. |
hrf |
Hemodynamic response function used for all conditions. |
duration |
Numeric scalar or vector of event durations. |
amplitude |
Numeric scalar or vector of event amplitudes. |
span |
Numeric scalar giving the HRF span in seconds. |
summate |
Logical; passed to [regressor()]. |
x |
A RegSet object |
grid |
Numeric vector of time points at which to evaluate |
precision |
Numeric precision for evaluation |
method |
Evaluation method |
sparse |
Logical whether to return sparse matrix |
... |
Additional arguments passed to evaluate |
An object of class 'RegSet' containing one 'Reg' per factor level.
# Create events for 3 conditions onsets <- c(10, 20, 30, 40, 50, 60) conditions <- factor(c("A", "B", "C", "A", "B", "C")) # Create regressor set rset <- regressor_set(onsets, conditions, hrf = HRF_SPMG1) # With durations and amplitudes rset2 <- regressor_set( onsets = onsets, fac = conditions, duration = 2, amplitude = c(1, 1.5, 0.8, 1, 1.5, 0.8), hrf = HRF_SPMG1 ) # Evaluate the regressor set times <- seq(0, 80, by = 0.1) design_matrix <- evaluate(rset, times)# Create events for 3 conditions onsets <- c(10, 20, 30, 40, 50, 60) conditions <- factor(c("A", "B", "C", "A", "B", "C")) # Create regressor set rset <- regressor_set(onsets, conditions, hrf = HRF_SPMG1) # With durations and amplitudes rset2 <- regressor_set( onsets = onsets, fac = conditions, duration = 2, amplitude = c(1, 1.5, 0.8, 1, 1.5, 0.8), hrf = HRF_SPMG1 ) # Evaluate the regressor set times <- seq(0, 80, by = 0.1) design_matrix <- evaluate(rset, times)
Generic function retrieving sampling times from a sampling frame or related object.
samples(x, ...) ## S3 method for class 'sampling_frame' samples(x, blockids = NULL, global = FALSE, ...)samples(x, ...) ## S3 method for class 'sampling_frame' samples(x, blockids = NULL, global = FALSE, ...)
x |
Object describing the sampling grid |
... |
Additional arguments passed to methods |
blockids |
Integer vector of block identifiers to include (default: all blocks) |
global |
Logical indicating whether to return global times (default: FALSE) |
Numeric vector of sample times
# Get sample times from a sampling frame sframe <- sampling_frame(blocklens = c(100, 120), TR = 2) samples(sframe, blockids = 1) # First block only samples(sframe, global = TRUE) # All blocks, global timing# Get sample times from a sampling frame sframe <- sampling_frame(blocklens = c(100, 120), TR = 2) samples(sframe, blockids = 1) # First block only samples(sframe, global = TRUE) # All blocks, global timing
sampling_frame describes the block structure and temporal sampling of an fMRI paradigm.A sampling_frame describes the block structure and temporal sampling of an fMRI paradigm.
sampling_frame(blocklens, TR, start_time = TR/2, precision = 0.1)sampling_frame(blocklens, TR, start_time = TR/2, precision = 0.1)
blocklens |
A numeric vector representing the number of scans in each block. |
TR |
A numeric value or vector representing the repetition time in seconds (i.e., the spacing between consecutive image acquisitions). When a vector is provided, its length must be 1 or equal to the number of blocks. |
start_time |
A numeric value or vector representing the offset of the first
scan of each block (default is |
precision |
A numeric value representing the discrete sampling interval used for convolution with the hemodynamic response function (default is 0.1). |
A list with class "sampling_frame" describing the block structure and temporal sampling of an fMRI paradigm.
frame <- sampling_frame(blocklens = c(100, 100, 100), TR = 2, precision = 0.5) # The relative time (with respect to the last block) in seconds of each sample/acquisition sam <- samples(frame) # The global time (with respect to the first block) of each sample/acquisition gsam <- samples(frame, global = TRUE) # Block identifiers for each acquisition can be retrieved using # blockids(frame)frame <- sampling_frame(blocklens = c(100, 100, 100), TR = 2, precision = 0.5) # The relative time (with respect to the last block) in seconds of each sample/acquisition sam <- samples(frame) # The global time (with respect to the first block) of each sample/acquisition gsam <- samples(frame, global = TRUE) # Block identifiers for each acquisition can be retrieved using # blockids(frame)
Discrete colour and fill scales built on hrf_palette(). Use them in your
own ggplot2 figures to match the output of plot_hrfs() and
plot_regressors().
scale_colour_hrf( type = c("categorical", "ordered"), ..., aesthetics = "colour" ) scale_color_hrf(type = c("categorical", "ordered"), ..., aesthetics = "colour") scale_fill_hrf(type = c("categorical", "ordered"), ..., aesthetics = "fill")scale_colour_hrf( type = c("categorical", "ordered"), ..., aesthetics = "colour" ) scale_color_hrf(type = c("categorical", "ordered"), ..., aesthetics = "colour") scale_fill_hrf(type = c("categorical", "ordered"), ..., aesthetics = "fill")
type |
Either |
... |
Further arguments passed to |
aesthetics |
The aesthetics the scale applies to. |
A ggplot2 scale object.
if (requireNamespace("ggplot2", quietly = TRUE)) { t <- seq(0, 24, by = 0.2) df <- data.frame( time = rep(t, 3), response = c(HRF_SPMG1(t), HRF_GAMMA(t), HRF_GAUSSIAN(t)), hrf = rep(c("SPMG1", "Gamma", "Gaussian"), each = length(t)) ) ggplot2::ggplot(df, ggplot2::aes(time, response, colour = hrf)) + ggplot2::geom_line() + scale_colour_hrf() }if (requireNamespace("ggplot2", quietly = TRUE)) { t <- seq(0, 24, by = 0.2) df <- data.frame( time = rep(t, 3), response = c(HRF_SPMG1(t), HRF_GAMMA(t), HRF_GAUSSIAN(t)), hrf = rep(c("SPMG1", "Gamma", "Gaussian"), each = length(t)) ) ggplot2::ggplot(df, ggplot2::aes(time, response, colour = hrf)) + ggplot2::geom_line() + scale_colour_hrf() }
Apply a temporal shift to a time series object. This function shifts the values in time while preserving the structure of the object. Common uses include:
Aligning regressors with different temporal offsets
Applying temporal derivatives to time series
Correcting for timing differences between signals
shift(x, ...) ## S3 method for class 'FeatureReg' shift(x, shift_amount, ...) ## S3 method for class 'Reg' shift(x, shift_amount, ...)shift(x, ...) ## S3 method for class 'FeatureReg' shift(x, shift_amount, ...) ## S3 method for class 'Reg' shift(x, shift_amount, ...)
x |
An object representing a time series or a time-based data structure |
... |
Additional arguments passed to methods |
shift_amount |
Numeric; amount to shift by (positive = forward, negative = backward) |
An object of the same class as the input, with values shifted in time:
Values are moved by the specified offset
Object structure and dimensions are preserved
Empty regions are filled with padding value
[regressor()], [evaluate()]
# Create a simple time series with events event_data <- data.frame( onsets = c(1, 10, 20, 30), run = c(1, 1, 1, 1) ) # Create regressor from events reg <- regressor( onsets = event_data$onsets, hrf = HRF_SPMG1, duration = 0, amplitude = 1 ) # Shift regressor forward by 2 seconds reg_forward <- shift(reg, shift_amount = 2) # Shift regressor backward by 1 second reg_backward <- shift(reg, shift_amount = -1) # Evaluate original and shifted regressors times <- seq(0, 50, by = 2) orig_values <- evaluate(reg, times) shifted_values <- evaluate(reg_forward, times)# Create a simple time series with events event_data <- data.frame( onsets = c(1, 10, 20, 30), run = c(1, 1, 1, 1) ) # Create regressor from events reg <- regressor( onsets = event_data$onsets, hrf = HRF_SPMG1, duration = 0, amplitude = 1 ) # Shift regressor forward by 2 seconds reg_forward <- shift(reg, shift_amount = 2) # Shift regressor backward by 1 second reg_backward <- shift(reg, shift_amount = -1) # Evaluate original and shifted regressors times <- seq(0, 50, by = 2) orig_values <- evaluate(reg, times) shifted_values <- evaluate(reg_forward, times)
Creates a regressor object for modeling a single trial event in an fMRI experiment. This is particularly useful for trial-wise analyses where each trial needs to be modeled separately. The regressor represents the predicted BOLD response for a single event using a specified hemodynamic response function (HRF).
single_trial_regressor( onsets, hrf = HRF_SPMG1, duration = 0, amplitude = 1, span = 24 )single_trial_regressor( onsets, hrf = HRF_SPMG1, duration = 0, amplitude = 1, span = 24 )
onsets |
the event onset in seconds, must be of length 1. |
hrf |
a hemodynamic response function, e.g. |
duration |
duration of the event (default is 0), must be length 1. |
amplitude |
scaling vector (default is 1), must be length 1. |
span |
the temporal window of the impulse response function (default is 24). |
This is a convenience wrapper around 'regressor' that ensures inputs have length 1.
A 'Reg' object (inheriting from 'regressor' and 'list').
# Create single trial regressor at 10 seconds str1 <- single_trial_regressor(onsets = 10, hrf = HRF_SPMG1) # Single trial with duration and custom amplitude str2 <- single_trial_regressor( onsets = 15, duration = 3, amplitude = 2, hrf = HRF_SPMG1 ) # Evaluate the response times <- seq(0, 40, by = 0.1) response <- evaluate(str1, times)# Create single trial regressor at 10 seconds str1 <- single_trial_regressor(onsets = 10, hrf = HRF_SPMG1) # Single trial with duration and custom amplitude str2 <- single_trial_regressor( onsets = 15, duration = 3, amplitude = 2, hrf = HRF_SPMG1 ) # Evaluate the response times <- seq(0, 40, by = 0.1) response <- evaluate(str1, times)