| 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.3.1 |
| Built: | 2026-05-19 10:03:22 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()
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. Since SPMG3 already includes first and second derivatives as basis functions, this method returns their derivatives (1st, 2nd, and 3rd derivatives of the original HRF).
## 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", "fft", "Rconv", "loop"), sparse = FALSE, normalize = FALSE, ... )evaluate(x, grid, ...) ## S3 method for class 'Reg' evaluate( x, grid, precision = 0.33, method = c("conv", "fft", "Rconv", "loop"), 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). If > 0, the HRF is evaluated over this duration (default: 0). |
precision |
The temporal resolution for evaluating responses when duration > 0 (default: 0.2). |
summate |
Logical; whether the HRF response should accumulate over the duration (default: TRUE). If FALSE, the convolution is averaged so the temporal profile is preserved but peak amplitude does not grow with duration. |
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
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, ... )gen_hrf( hrf, lag = 0, width = 0, precision = 0.1, half_life = Inf, summate = TRUE, normalize = FALSE, name = NULL, span = NULL, ... )
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. |
... |
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, ... )getHRF( name = "spmg1", nbasis = 5, span = 24, lag = 0, width = 0, summate = TRUE, normalize = FALSE, ... )
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) |
... |
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 |
When used in a GLM, the estimated coefficient represents a (weighted) average
of the data within the specified time window. If normalize = TRUE, the
coefficient directly estimates the mean signal in that window.
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") # Normalized boxcar - coefficient will estimate mean signal in window hrf2 <- hrf_boxcar(width = 5, normalize = TRUE) # integral is now 1, so beta estimates mean(Y[0:5]) # 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") # Normalized boxcar - coefficient will estimate mean signal in window hrf2 <- hrf_boxcar(width = 5, normalize = TRUE) # integral is now 1, so beta estimates mean(Y[0:5]) # 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 the B-spline representation of an HRF (hemodynamic response function) at given time points 't'.
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 20. |
N |
An integer representing the number of basis functions. Default value is 5. |
degree |
An integer representing the degree of the spline. Default value is 3. |
... |
Additional arguments passed to 'splines::bs'. |
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 |
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")
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 = 0.0833)hrf_spmg1(t, P1 = 5, P2 = 15, A1 = 0.0833)
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 .0833 |
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 extracting weighted averages of data at specific time points.
When normalize = TRUE and the HRF is used in a GLM, the estimated
coefficient represents a weighted mean of the data at the specified times.
There are two ways to specify the temporal structure:
width + weights: Weights are evenly spaced from 0 to width
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: 6s window with 4 evenly-spaced weights (at 0, 2, 4, 6s) 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 - coefficient estimates weighted mean of signal 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: 6s window with 4 evenly-spaced weights (at 0, 2, 4, 6s) 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 - coefficient estimates weighted mean of signal 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()
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 )
Creates a new HRF object whose output is scaled such that the maximum absolute value of the response is 1.
normalise_hrf(hrf)normalise_hrf(hrf)
hrf |
The HRF object (of class 'HRF') to normalise. |
For multi-basis HRFs, each basis function (column) is normalised independently.
A new HRF object representing the normalised function.
Other HRF_decorator_functions:
block_hrf(),
lag_hrf()
# Create a gaussian HRF with a peak value != 1 gauss_unnorm <- as_hrf(function(t) 5 * dnorm(t, 6, 2), name="unnorm_gauss") # Normalise it gauss_norm <- normalise_hrf(gauss_unnorm) t_vals <- seq(0, 20, by = 0.1) max(gauss_unnorm(t_vals)) # Peak is > 1 max(gauss_norm(t_vals)) # Peak is 1# Create a gaussian HRF with a peak value != 1 gauss_unnorm <- as_hrf(function(t) 5 * dnorm(t, 6, 2), name="unnorm_gauss") # Normalise it gauss_norm <- normalise_hrf(gauss_unnorm) t_vals <- seq(0, 20, by = 0.1) max(gauss_unnorm(t_vals)) # Peak is > 1 max(gauss_norm(t_vals)) # Peak is 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)
Creates a comparison plot of multiple HRF objects. This function provides a convenient way to visualize different HRFs on the same plot, with options for normalization and customization. Uses ggplot2 if available for publication-quality figures, otherwise falls back to base R graphics.
plot_hrfs( ..., time = NULL, normalize = FALSE, labels = NULL, title = NULL, subtitle = NULL, use_ggplot = TRUE )plot_hrfs( ..., time = NULL, normalize = FALSE, labels = NULL, title = NULL, subtitle = NULL, use_ggplot = TRUE )
... |
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 for each HRF. If NULL (default), uses the 'name' attribute of each HRF, or "HRF_1", "HRF_2", etc. |
title |
Character string for the plot title. If NULL (default), uses "HRF Comparison". |
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. |
Invisibly returns a data frame in long format with columns 'time', 'HRF', and 'response'. If use_ggplot is TRUE and ggplot2 is available, also returns a ggplot object as an attribute 'plot'.
# Compare canonical HRFs plot_hrfs(HRF_SPMG1, HRF_GAMMA, HRF_GAUSSIAN) # Compare with custom labels plot_hrfs(HRF_SPMG1, HRF_GAMMA, labels = c("SPM Canonical", "Gamma")) # Normalize for shape comparison plot_hrfs(HRF_SPMG1, HRF_GAMMA, HRF_GAUSSIAN, normalize = TRUE, title = "HRF Shape Comparison", subtitle = "All HRFs normalized to peak at 1") # Compare blocked HRFs with different durations hrf_1s <- block_hrf(HRF_SPMG1, width = 1) hrf_3s <- block_hrf(HRF_SPMG1, width = 3) hrf_5s <- block_hrf(HRF_SPMG1, width = 5) plot_hrfs(hrf_1s, hrf_3s, hrf_5s, labels = c("1s duration", "3s duration", "5s duration"), title = "Effect of Event Duration on HRF") # 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) # Compare with custom labels plot_hrfs(HRF_SPMG1, HRF_GAMMA, labels = c("SPM Canonical", "Gamma")) # Normalize for shape comparison plot_hrfs(HRF_SPMG1, HRF_GAMMA, HRF_GAUSSIAN, normalize = TRUE, title = "HRF Shape Comparison", subtitle = "All HRFs normalized to peak at 1") # Compare blocked HRFs with different durations hrf_1s <- block_hrf(HRF_SPMG1, width = 1) hrf_3s <- block_hrf(HRF_SPMG1, width = 3) hrf_5s <- block_hrf(HRF_SPMG1, width = 5) plot_hrfs(hrf_1s, hrf_3s, hrf_5s, labels = c("1s duration", "3s duration", "5s duration"), title = "Effect of Event Duration on HRF") # Use base R graphics instead of ggplot2 plot_hrfs(HRF_SPMG1, HRF_GAMMA, use_ggplot = FALSE)
Creates a comparison plot of multiple regressor objects. This function provides a convenient way to visualize different regressors on the same plot, with options for showing event onsets and customization. Uses ggplot2 if available for publication-quality figures, otherwise falls back to base R graphics.
plot_regressors( ..., grid = NULL, labels = NULL, title = NULL, subtitle = NULL, show_onsets = "first", onset_alpha = 0.3, precision = 0.33, use_ggplot = TRUE )plot_regressors( ..., grid = NULL, labels = NULL, title = NULL, subtitle = NULL, show_onsets = "first", onset_alpha = 0.3, precision = 0.33, use_ggplot = TRUE )
... |
Regressor objects to compare. Can be passed as individual arguments or as a named list. |
grid |
Numeric vector of time points for evaluation. If NULL (default), automatically generates a grid covering all regressors. |
labels |
Character vector of labels for each regressor. If NULL (default), uses "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, show onset lines for all regressors. If "first", show only for the first regressor. If FALSE, hide onsets. Default is "first". |
onset_alpha |
Alpha transparency for onset lines. Default is 0.3. |
precision |
Numeric sampling precision for HRF evaluation. Default is 0.33. |
use_ggplot |
Logical; if TRUE and ggplot2 is available, use ggplot2 for plotting. If FALSE, use base R graphics. Default is TRUE. |
Invisibly returns a data frame in long format with columns 'time', 'Regressor', and 'response'.
# 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")) # Compare regressors with different event timings reg_fast <- regressor(seq(0, 60, by = 10), HRF_SPMG1) reg_slow <- regressor(seq(0, 60, by = 20), HRF_SPMG1) plot_regressors(reg_fast, reg_slow, labels = c("Fast (10s ISI)", "Slow (20s ISI)"), title = "Effect of Inter-Stimulus Interval") # Compare original vs shifted regressor reg_orig <- regressor(c(10, 30, 50), HRF_SPMG1) reg_shifted <- shift(reg_orig, 5) plot_regressors(reg_orig, reg_shifted, labels = c("Original", "Shifted +5s"))# 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")) # Compare regressors with different event timings reg_fast <- regressor(seq(0, 60, by = 10), HRF_SPMG1) reg_slow <- regressor(seq(0, 60, by = 20), HRF_SPMG1) plot_regressors(reg_fast, reg_slow, labels = c("Fast (10s ISI)", "Slow (20s ISI)"), title = "Effect of Inter-Stimulus Interval") # Compare original vs shifted regressor reg_orig <- regressor(c(10, 30, 50), HRF_SPMG1) reg_shifted <- shift(reg_orig, 5) plot_regressors(reg_orig, reg_shifted, labels = c("Original", "Shifted +5s"))
Creates a visualization of an HRF object. For single-basis HRFs, shows the response curve with peak annotation. For multi-basis HRFs (e.g., HRF_SPMG3), shows all basis functions on the same plot.
## 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 and amplitude on the plot. |
... |
Additional arguments passed to underlying plot functions. |
Invisibly returns a data frame with the time and response values (useful for further customization).
# 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))
Creates a visualization of a regressor object showing the predicted BOLD response over time. Optionally displays event onsets as vertical lines.
## S3 method for class 'Reg' plot( x, grid = NULL, show_onsets = TRUE, onset_color = "red", onset_alpha = 0.5, precision = 0.33, ... )## S3 method for class 'Reg' plot( x, grid = NULL, show_onsets = TRUE, onset_color = "red", onset_alpha = 0.5, precision = 0.33, ... )
x |
A 'Reg' object created by 'regressor()'. |
grid |
Numeric vector of time points for evaluation. If NULL (default), automatically generates a grid from 0 to max(onsets) + span with step 0.5s. |
show_onsets |
Logical; if TRUE (default), show vertical dashed lines at event onset times. |
onset_color |
Color for onset lines. Default is "red". |
onset_alpha |
Alpha transparency for onset lines. Default is 0.5. |
precision |
Numeric sampling precision for HRF evaluation. Default is 0.33. |
... |
Additional arguments passed to underlying plot functions. |
Invisibly returns a data frame with the time and response values.
# 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 )# 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 )
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 )regressor( onsets, hrf = HRF_SPMG1, duration = 0, amplitude = 1, span = 40, summate = 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. |
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 automatically filtered out.
## 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", "fft", "Rconv", "loop"), 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", "fft", "Rconv", "loop"), 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", "fft", "Rconv", "loop"), 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", "fft", "Rconv", "loop"), 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)
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 'Reg' shift(x, shift_amount, ...)shift(x, ...) ## 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)