---
title: Overview and Workflow
author: Bradley R. Buchsbaum
date: '`r Sys.Date()`'
output:
  rmarkdown::html_vignette:
    toc: yes
    toc_depth: 2.0
    css: albers.css
    includes:
      in_header: albers-header.html
params:
  family: red
  preset: homage
resource_files:
- albers.css
- albers.js
- albers-header.html

vignette: |
  %\VignetteIndexEntry{Overview and Workflow}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
if (requireNamespace("ggplot2", quietly = TRUE) && requireNamespace("albersdown", quietly = TRUE)) ggplot2::theme_set(albersdown::theme_albers(family = params$family, preset = params$preset))
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width=7, fig.height=5, message=FALSE, warning=TRUE)
library(fmrireg)
library(fmrihrf)
```

```{r albers-classes, echo=FALSE, results='asis'}
cat(sprintf(
  paste0(
    '<script>document.addEventListener("DOMContentLoaded",function(){',
    'document.body.classList.remove("palette-red","palette-lapis","palette-ochre","palette-teal","palette-green","palette-violet","preset-homage","preset-study","preset-structural","preset-adobe","preset-midnight");',
    'document.body.classList.add("palette-%s","preset-%s");',
    '});</script>'
  ),
  params$family,
  params$preset
))
```

## 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.

```{r quick-fit-data}
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
)
```

```{r quick-fit}
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")
```

```{r check-quick-fit, include=FALSE}
stopifnot(
  all(is.finite(coef(fit))),
  max(abs(t(coef(fit)) - quick_truth)) < 0.2,
  all(is.finite(as.matrix(quick_contrast_numeric[-1]))),
  !any(quick_contrast$p_value == "0")
)
```

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.

### 1. Data Representation

Before modeling, you need to represent your fMRI data and its structure. `fmrireg` uses objects like:

*   `fmri_frame` (from the `fmridataset` package): an observations-by-features container that carries the data (in memory or file-backed), its feature space (voxel grid, parcels, or latent components), the TR and run structure as observation metadata, and the experimental event table. Build one with `matrix_frame()`, `neurovec_frame()`, `nifti_frame()`, or `latent_frame()`.
*   `sampling_frame`: Defines the temporal structure – the number of scans per run (`blocklens`) and the repetition time (`TR`).

```{r sampling_frame_example}
# Example: one run, 120 scans, TR=2s
sframe <- sampling_frame(blocklens = 120, TR = 2)
sframe
```

### 2. Model Building

#### a) Hemodynamic Response Function (HRF)

The HRF models the BOLD signal response to a brief neural event. `fmrireg` provides various pre-defined HRFs (e.g., `HRF_SPMG1`, `HRF_GAUSSIAN`) and functions to create or modify them (`gen_hrf`, `hrf_blocked`, `hrf_lagged`).

```{r hrf_example, fig.alt = "Line plot of the SPM Canonical HRF (BOLD response over time in seconds)"}
# 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)")
```

*   **Related Topic:** Hemodynamic Response Functions

#### b) Event Model (`event_model`)

This defines the task-related part of your GLM. You specify how experimental events (onsets, conditions, durations) map to regressors using a formula syntax. The `hrf()` or `afni_hrf()` functions within the formula link event variables to specific HRF shapes.

```{r event_design_table}
# 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)
)
```

```{r event_model_example}
# 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
```

*   **Related Topic:** Event Models

#### c) Baseline Model (`baseline_model`)

This models non-task-related components like scanner drift, run intercepts, and other nuisance factors (e.g., motion parameters).

```{r baseline_model_example}
# 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)
```

*   **Related Topic:** Baseline Models

#### d) Contrasts (`pair_contrast`, `unit_contrast`, etc.)

Contrasts define specific comparisons between condition estimates.

```{r contrast_example}
# 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)
```

*   **See Vignette:** `vignette("a_05_contrasts", package = "fmrireg")`

### 3. Model Estimation (`fmri_model`, `fmri_lm`)

You combine the event and baseline models into a full `fmri_model` and then fit it to your `fmri_frame` using estimation functions like `fmri_lm` (for standard GLM) or `estimate_betas` (for single-trial betas).

```{r fmri_model_example}
# 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"
)
```

```{r check-object-fit, include=FALSE}
stopifnot(
  all(is.finite(coef(overview_fit))),
  identical(overview_coef_table$condition, c("face", "scene")),
  ncol(overview_coef_matrix) == length(levels(design_table$stim_type)),
  nrow(overview_coef_matrix) == 3L
)
```

This second fit uses the object-building route rather than the compact formula
route shown first. See the linear-model vignette for AR prewhitening, robust
row weighting, uncertainty estimates, and model diagnostics.

*   **See Vignettes:**
    *   `vignette("a_09_linear_model", package = "fmrireg")`
    *   Linear Models and Trialwise Analysis

### 4. Advanced Topics

`fmrireg` also supports:

*   **Simulation:** Creating artificial fMRI data (`simulate_fmri_matrix`).
    *   `vignette("a_08_simulation", package = "fmrireg")`
*   **Low-level Regressors:** Manually creating and manipulating regressors (`regressor`).

#### Acceleration: Sketched / Nyström Analysis

For large problems, you can accelerate voxelwise GLM estimation with:
- Temporal sketching (CountSketch, Gaussian or SRHT) to reduce timepoints, with standard errors that account for the sketch.
- Nyström/landmark spatial extension to solve on a subset of voxels and extend with parcel‑aware kernels.

See `vignette("sketched-ar", package = "fmrireg")` for guidance, tradeoffs, and practical defaults.

#### Group-Level Meta-Analysis

Combine subject‑level maps or ROI summaries using fixed/random‑effects meta‑regression:
- Inputs: HDF5/NIfTI (beta+SE or t+df) or ROI CSV.
- Random‑effects: Paule–Mandel/DerSimonian–Laird/REML alias; optional robust Huber.
- t‑only combine: Stouffer/Fisher/Lancaster (supports equal/IVW/custom weights where applicable).
- Exact contrasts: at fit‑time (contrasts=...) or post‑hoc via `return_cov = "tri"` + `contrast()`.
- Multiple comparisons: BH FDR or spatial FDR (block/parcel aware).

See `vignette("group_analysis", package = "fmrireg")` for end‑to‑end examples.

## 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