Overview and Workflow

Introduction

If you work with fMRI data, you already know the core analysis is multiple regression applied voxel-by-voxel — but specifying the design matrix, convolving events with hemodynamic response functions, and managing multi-run structures can be tedious and error-prone.

The fmrireg package gives you a concise, formula-based interface for these tasks, inspired by lm() and glm(). You describe your experimental design and HRF choices in a formula, and fmrireg handles convolution, baseline modeling, contrast computation, and model fitting.

This vignette starts with a complete, small fit and then unpacks the objects behind it. You begin with a time-by-voxel matrix and an event table; you finish with one coefficient estimate per condition and voxel.

A complete fit first

For a compact example, we create two voxels with different known responses to faces and scenes. In your own analysis, bold would be the matrix extracted from preprocessed images or regions of interest.

quick_events <- data.frame(
  onset = seq(10, 80, by = 10),
  condition = factor(rep(c("face", "scene"), 4)),
  run = 1L
)
quick_frame <- sampling_frame(blocklens = 80, TR = 2)
quick_event_model <- event_model(
  onset ~ hrf(condition), block = ~ run,
  data = quick_events, sampling_frame = quick_frame
)
quick_design <- as.matrix(design_matrix(quick_event_model))
quick_truth <- matrix(
  c(1.0, 0.6, 0.7, 1.1), nrow = ncol(quick_design),
  dimnames = list(colnames(quick_design), c("voxel1", "voxel2"))
)
# Set the noise level relative to the regressors' peak height, so the
# signal-to-noise ratio does not depend on the HRF's raw amplitude.
noise_sd <- 0.11 * max(abs(quick_design))
set.seed(2026)
bold <- quick_design %*% quick_truth +
  matrix(rnorm(80 * 2, sd = noise_sd), 80, 2)
quick_dataset <- matrix_frame(
  bold, TR = 2, run_length = 80, event_table = quick_events
)

face_vs_scene <- pair_contrast(
  ~ condition == "face", ~ condition == "scene",
  name = "face_vs_scene"
)
fit <- fmri_lm(
  onset ~ hrf(condition, contrasts = face_vs_scene), block = ~ run,
  dataset = quick_dataset
)

coefficient_table <- data.frame(
  condition = levels(quick_events$condition),
  truth_voxel1 = quick_truth[, 1],
  estimate_voxel1 = as.numeric(coef(fit)[1, ]),
  truth_voxel2 = quick_truth[, 2],
  estimate_voxel2 = as.numeric(coef(fit)[2, ])
)
knitr::kable(coefficient_table, digits = 2,
             row.names = FALSE,
             caption = "Known and recovered condition effects")

quick_contrast_numeric <- data.frame(
  voxel = colnames(quick_truth),
  estimate = as.numeric(coef(fit, type = "contrasts")$face_vs_scene),
  std_error = as.numeric(standard_error(fit, type = "contrasts")$face_vs_scene),
  statistic = as.numeric(stats(fit, type = "contrasts")$face_vs_scene),
  p_value = as.numeric(p_values(fit, type = "contrasts")$face_vs_scene)
)
quick_contrast <- transform(
  quick_contrast_numeric,
  p_value = format.pval(p_value, digits = 3, eps = 1e-4)
)
knitr::kable(quick_contrast, digits = 3,
             row.names = FALSE,
             caption = "Face-minus-scene contrast with uncertainty")

The first table separates the known effect sizes from their estimates. The second table asks a scientific question—face minus scene—and reports its estimate, standard error, test statistic, and p-value for each voxel. Here the first voxel favors faces while the second favors scenes, recovering the pattern used to generate the data. The next sections unpack the dataset, event, baseline, and contrast objects used for larger analyses.

Key Components and Workflow

A typical fMRI analysis workflow with fmrireg involves these main stages:

1. Data Representation: Defining how your fMRI data and experimental design are organized.
2. Model Building: Specifying the regression model using formulas, including:
   • Defining Hemodynamic Response Functions (HRFs).
   • Creating event-related regressors (event_model).
   • Modeling baseline components like drift and nuisance factors (baseline_model).
   • Defining contrasts of interest (contrast).
3. Model Estimation: Fitting the specified model to the data (e.g., using fmri_lm).
4. Advanced Topics: Exploring single-trial estimation, simulation, group analysis, or accelerated engines.

# Example: one run, 120 scans, TR=2s
sframe <- sampling_frame(blocklens = 120, TR = 2)
sframe
#> Sampling Frame
#> ==============
#> 
#> Structure:
#>   1 block
#>   Total scans: 120
#> 
#> Timing:
#>   TR: 2 s
#>   Precision: 0.1 s
#> 
#> Duration:
#>   Total time: 240.0 s

# Plot the SPM canonical HRF
time <- seq(0, 24, by = 0.2)
plot(time, HRF_SPMG1(time), type = 'l', xlab = "Time (s)", ylab = "BOLD Response",
     main = "SPM Canonical HRF (HRF_SPMG1)")

# A simple design table with two stimulus types
design_table <- data.frame(
  onset = c(10, 40, 70, 100, 130, 160, 190, 220),
  stim_type = factor(rep(c("face", "scene"), 4)),
  run_id = rep(1, 8)
)

# Model the effect of stim_type, convolved with the SPM canonical HRF
event_mod <- event_model(onset ~ hrf(stim_type),
                         data = design_table,
                         block = ~ run_id,
                         sampling_frame = sframe)
event_mod

# Example: Cubic B-spline drift model + run intercepts
# (Requires 'sframe' defined earlier)
base_mod <- baseline_model(basis = "bs", degree = 3, sframe = sframe)
print(base_mod)
#> ================================================
#>            Baseline Model                       
#> ================================================
#>   Drift Components                           
#>     * Name: baseline_bs_3                
#>     * Basis type: bs                     
#>     * Degree: 3                          
#>     * Drift columns: 3                   
#> 
#>   Additional Components                    
#>     * Constant columns: 1                
#>     * Nuisance columns: 0                
#> 
#>   Model Summary                            
#>     * Total columns: 4                   
#> 
#>   Design Matrix Preview                    
#>      0.000  0.000  0.000  1.000          
#>      0.025  0.000  0.000  1.000          
#>      0.049  0.001  0.000  1.000          
#> ================================================

# Example: Contrast comparing 'face' vs 'scene' levels of 'stim_type'
face_vs_scene <- pair_contrast(~ stim_type == "face", ~ stim_type == "scene", name = "FvS")
print(face_vs_scene)
#> contrast: FvS 
#>  A:  ~stim_type == "face" 
#>  B:  ~stim_type == "scene"

# Create three example voxel series over the same one-run design
set.seed(42)
overview_data <- matrix_frame(
  matrix(rnorm(120 * 3), 120, 3),
  TR = 2, run_length = 120, event_table = design_table
)

# Combine the design with its dataset, then fit it
full_model <- fmri_model(event_mod, base_mod, dataset = overview_data)
overview_fit <- fmri_lm(full_model)
overview_coef_matrix <- coef(overview_fit) # voxels x conditions
overview_coef_table <- data.frame(
  condition = levels(design_table$stim_type),
  voxel_1 = overview_coef_matrix[1, ],
  voxel_2 = overview_coef_matrix[2, ],
  voxel_3 = overview_coef_matrix[3, ],
  check.names = FALSE
)
knitr::kable(
  overview_coef_table, digits = 2, row.names = FALSE,
  caption = "Condition coefficients from the object-building route"
)

Getting Started

This overview provides a map of the fmrireg package. For detailed explanations and examples, please refer to the specific vignettes linked above.

Next Steps

• vignette("a_10_dataset", package = "fmrireg") — Dataset management
• vignette("a_08_simulation", package = "fmrireg") — Build controlled examples
• vignette("a_09_linear_model", package = "fmrireg") — Fit and diagnose a GLM
• vignette("a_05_contrasts", package = "fmrireg") — Test planned hypotheses