---
title: Building fMRI Regressors
author: Bradley R. Buchsbaum
date: '`r Sys.Date()`'
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Building fMRI Regressors}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
params:
  family: red
  preset: homage
css: albers.css
resource_files:
- albers.css
- albers.js
includes:
  in_header: |-
    <script src="albers.js"></script>
    <script>document.addEventListener('DOMContentLoaded',function(){document.body.classList.add('palette-red');});</script>

---

```{r setup, include = FALSE}
if (requireNamespace("ggplot2", quietly = TRUE)) ggplot2::theme_set(ggplot2::theme_minimal())
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7, 
  fig.height = 5,
  message = FALSE,
  warning = FALSE
)
library(fmrihrf)
library(dplyr)
library(ggplot2)
library(tidyr)
```

## Introduction: What is a Regressor?

In fMRI analysis, a **regressor** (or predictor) represents the expected BOLD signal timecourse associated with a specific experimental condition or event type. It's typically created by convolving a series of event onsets (often represented as delta functions or "sticks") with a hemodynamic response function (HRF).

`fmrihrf` provides the `regressor()` function to easily create these objects from event timings and an HRF. While these regressor objects are often constructed automatically by modeling functions in other packages, this vignette explores how to create and manipulate them directly, offering finer control over the model components.

## Basic Regressor from Event Onsets

Suppose we have a simple event-related fMRI design with stimuli presented every 12 seconds. We want to model these events using the SPM canonical HRF (`HRF_SPMG1`). The events are brief, so we model them with a duration of 0 seconds (instantaneous).

```{r basic_regressor}
# Define event onsets
onsets <- seq(0, 10 * 12, by = 12)

# Create the regressor object
# Uses HRF_SPMG1 by default if no hrf is specified
# Duration is 0 by default
reg1 <- regressor(onsets = onsets, hrf = HRF_SPMG1)

# Print the regressor object to see its properties (uses new print.Reg method)
print(reg1)

# Access components using helper functions
head(onsets(reg1))
nbasis(reg1)
```

## Evaluating and Plotting a Regressor

A `regressor` object stores the event information but doesn't automatically compute the timecourse. To get the predicted BOLD signal at specific time points (e.g., corresponding to scan acquisition times), we use the `evaluate()` function.

```{r evaluate_plot_basic}
# Define a time grid corresponding to scan times (e.g., TR=2s)
TR <- 2
scan_times <- seq(0, 140, by = TR)

# Plot the regressor using the plot() method
# This automatically evaluates and plots with onset markers
plot(reg1, grid = scan_times)
```

## Varying Event Durations

Sometimes events have different durations. The `duration` argument in `regressor()` can take a vector matching the length of `onsets`.

```{r varying_duration}
# Example onsets and durations
onsets_var_dur <- seq(0, 5 * 12, length.out = 6)
durations_var <- 1:length(onsets_var_dur) # Durations increase from 1s to 6s

# Create regressor with varying durations
reg_var_dur <- regressor(onsets_var_dur, HRF_SPMG1, duration = durations_var)
print(reg_var_dur)

# Plot the regressor
scan_times_dur <- seq(0, max(onsets_var_dur) + 30, by = TR)
plot(reg_var_dur, grid = scan_times_dur)
```

### Duration and Summation

By default (`summate=TRUE`), the predicted response accumulates if events overlap or have extended duration. Setting `summate=FALSE` preserves the same temporal profile but the peak amplitude does not grow with duration.

```{r duration_no_summate}
# Create regressor with varying durations, summate=FALSE
reg_var_dur_nosum <- regressor(onsets_var_dur, HRF_SPMG1,
                               duration = durations_var, summate = FALSE)

# Compare summating vs non-summating using plot_regressors()
plot_regressors(reg_var_dur, reg_var_dur_nosum,
                labels = c("summate=TRUE", "summate=FALSE"),
                grid = scan_times_dur,
                title = "Effect of Summation on Response",
                subtitle = "Same events with varying durations")
```

## Varying Event Amplitudes (Parametric Modulation)

We can model variations in event intensity or some associated parameter by providing an `amplitude` vector. This creates a *parametric regressor* where the height of the HRF for each event is scaled by the corresponding amplitude value.

```{r parametric_modulation}
# Example onsets and amplitudes (e.g., representing task difficulty)
onsets_amp <- seq(0, 10 * 12, length.out = 11)
amplitudes_raw <- 1:length(onsets_amp)

# It's common practice to center the modulator
amplitudes_scaled <- scale(amplitudes_raw, center = TRUE, scale = FALSE)

# Create the parametric regressor
reg_amp <- regressor(onsets_amp, HRF_SPMG1, amplitude = amplitudes_scaled)
print(reg_amp)

# Plot the parametric regressor
scan_times_amp <- seq(0, max(onsets_amp) + 30, by = TR)
plot(reg_amp, grid = scan_times_amp)
```

## Combining Duration and Amplitude Modulation

You can provide both `duration` and `amplitude` vectors to model events that vary in both aspects.

```{r duration_amplitude}
set.seed(123)
onsets_comb <- seq(0, 10 * 12, length.out = 11)
amps_comb <- scale(1:length(onsets_comb), center = TRUE, scale = FALSE)
durs_comb <- sample(1:5, length(onsets_comb), replace = TRUE)

reg_comb <- regressor(onsets_comb, HRF_SPMG1, 
                      amplitude = amps_comb, duration = durs_comb)
print(reg_comb)

# Plot the combined regressor
scan_times_comb <- seq(0, max(onsets_comb) + 30, by = TR)
plot(reg_comb, grid = scan_times_comb)
```

## Regressors with HRF Basis Sets

If you use an HRF object with multiple basis functions (e.g., `HRF_SPMG3`, `HRF_BSPLINE`), the `regressor` object will represent multiple timecourses, one for each basis function. `evaluate()` will return a matrix.

```{r basis_set_regressor}
# Use a B-spline basis set
onsets_basis <- seq(0, 10 * 12, length.out = 11)
hrf_basis <- HRF_BSPLINE # Uses N=5 basis functions by default

reg_basis <- regressor(onsets_basis, hrf_basis)
print(reg_basis)
nbasis(reg_basis) # Should be 5

# Evaluate - this returns a matrix
scan_times_basis <- seq(0, max(onsets_basis) + 30, by = TR)
pred_basis_matrix <- evaluate(reg_basis, scan_times_basis)
dim(pred_basis_matrix) # rows = time points, cols = basis functions

# Plot using the plot() method - automatically handles multi-basis
plot(reg_basis, grid = scan_times_basis)
```

## Shifting Regressors

You can temporally shift all onsets within a regressor using the `shift()` method.

```{r shift_regressor}
# Original regressor
reg_orig <- regressor(onsets = c(10, 30, 50), hrf = HRF_SPMG1)

# Shifted regressor (delay by 5 seconds)
reg_shifted <- shift(reg_orig, shift_amount = 5)

onsets(reg_orig)
onsets(reg_shifted) # Onsets are now 15, 35, 55

# Compare original and shifted using plot_regressors()
scan_times_shift <- seq(0, 80, by = TR)
plot_regressors(reg_orig, reg_shifted,
                labels = c("Original", "Shifted +5s"),
                grid = scan_times_shift,
                show_onsets = TRUE,  # Show onsets for both
                title = "Shifting a Regressor")
```