---
title: Voxel-wise HRF Modeling with fmrilss
output:
  rmarkdown::html_vignette:
    mathjax: default
    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{Voxel-wise HRF Modeling with fmrilss}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup-opts, 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,
  dev = "png",
  dpi = 150,
  warning = FALSE,
  message = FALSE
)
```

```{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
))
```

## Why HRF Variability Matters

Vascular properties, neurovascular coupling, and acquisition protocols all influence the hemodynamic response. Relying on a single canonical HRF shape can bias trial-wise estimates, especially across brain regions or clinical populations.

This vignette shows how to estimate voxel-specific HRFs and incorporate them into LSS. You should be familiar with the core LSS workflow (`vignette("fmrilss")`) and `fmrihrf` basics.

**Alternative approach:** For a library-constrained method that is usually faster and more stable, see `vignette("sbhm")`.

## Setup

```{r setup}
library(fmrilss)
library(fmrihrf)
set.seed(123)
```

## Simulate Data with Variable HRFs

### Experiment parameters

You need a sampling frame, jittered onsets, and a small set of voxels whose HRFs differ.

```{r sim-params}
n_time <- 200
n_vox  <- 5
TR     <- 1.0
sframe <- fmrihrf::sampling_frame(blocklens = n_time, TR = TR)
grid   <- fmrihrf::samples(sframe, global = TRUE)
```

```{r sim-onsets}
isi    <- runif(200, min = 3, max = 9)
onsets <- cumsum(c(10, isi))
onsets <- onsets[onsets < (n_time - 20)]
n_trials <- length(onsets)
```

### Voxel-specific HRFs

Each voxel gets a slightly shifted and scaled version of the canonical HRF. This mimics spatial variation in vascular properties.

We create each voxel's HRF by wrapping the canonical shape with a peak shift and width scaling.

```{r sim-hrfs}
voxel_hrfs <- lapply(1:n_vox, function(v) {
  peak_shift  <- (v - 3) * 0.5        # -1 to +1 s
  width_scale <- 1 + (v - 3) * 0.1    #  0.8 to 1.2
  fmrihrf::HRF(
    fun   = function(t) HRF_SPMG1(t - peak_shift) * width_scale,
    name  = paste0("voxel_", v),
    span  = attr(HRF_SPMG1, "span"),
    nbasis = 1L
  )
})
```

### Generate signal and noise

For each voxel, you convolve the trial onsets with that voxel's HRF, scale by true betas, then add AR(1) noise.

```{r sim-betas}
true_betas <- matrix(rnorm(n_trials * n_vox, mean = 1, sd = 0.3),
                     nrow = n_trials, ncol = n_vox)
```

```{r sim-signal}
Y <- matrix(0, n_time, n_vox)
for (v in 1:n_vox) {
  rset <- fmrihrf::regressor_set(
    onsets = onsets, fac = factor(seq_len(n_trials)),
    hrf = voxel_hrfs[[v]], duration = 0, span = 30, summate = FALSE
  )
  Xv <- fmrihrf::evaluate(rset, grid = grid, precision = 0.1, method = "conv")
  Y[, v] <- as.matrix(Xv) %*% true_betas[, v]
}
```

```{r sim-noise}
noise_sd <- 0.5; ar_coef <- 0.3
for (v in 1:n_vox) {
  e <- rnorm(n_time, sd = noise_sd)
  noise <- as.numeric(stats::filter(e, filter = ar_coef, method = "recursive"))
  Y[, v] <- Y[, v] + noise
}
colnames(Y) <- paste0("V", 1:n_vox)
```

### Visualise the design

```{r design-heatmap, fig.width=8, fig.height=4, fig.alt="Trial-wise design heatmap."}
rset_vis <- fmrihrf::regressor_set(
  onsets = onsets, fac = factor(seq_len(n_trials)),
  hrf = HRF_SPMG1, duration = 0, span = 30, summate = FALSE)
X_vis <- as.matrix(fmrihrf::evaluate(rset_vis, grid = grid, precision = 0.1, method = "conv"))
image(seq_len(nrow(X_vis)), seq_len(ncol(X_vis)), X_vis,
      col = hcl.colors(64, "BluGrn"), xlab = "Time (TR)", ylab = "Trial",
      main = "Trial-wise design (jittered ISIs)")
```

## Standard LSS with Canonical HRF

Standard LSS assumes every voxel shares the same canonical HRF. When that assumption is wrong, you get biased betas.

```{r standard-lss}
rset_can <- fmrihrf::regressor_set(
  onsets = onsets, fac = factor(seq_len(n_trials)),
  hrf = HRF_SPMG1, duration = 0, span = 30, summate = FALSE)
X_can <- as.matrix(fmrihrf::evaluate(rset_can, grid = grid, precision = 0.1, method = "conv"))
standard_betas <- lss(Y, X_can, method = "r_optimized")
```

## Estimate Voxel-Specific HRFs

### Multi-basis GLM

The SPMG3 basis set includes the canonical HRF plus its temporal and dispersion derivatives. Fitting a GLM with this basis set lets you estimate how each voxel's HRF deviates from canonical.

```{r multibasis-design}
rset_mb <- fmrihrf::regressor_set(
  onsets = onsets, fac = factor(rep(1, n_trials)),
  hrf = HRF_SPMG3, duration = 0, span = 30, summate = TRUE)
X_mb <- as.matrix(fmrihrf::evaluate(rset_mb, grid = grid, precision = 0.1, method = "conv"))
```

```{r estimate-weights}
hrf_weights <- sapply(1:n_vox, function(v) coef(lm(Y[, v] ~ X_mb - 1)))
cat("Basis weights (3 x", n_vox, "voxels):\n")
print(round(hrf_weights, 2))
```

The first row captures the canonical amplitude; rows 2--3 capture latency and width shifts.

## Apply Voxel-Specific HRFs in LSS

With basis weights in hand, you can build a voxel-specific design by weighting the three basis columns for each trial. Notice that this is a per-voxel loop: each voxel gets its own tailored design matrix.

```{r voxel-lss-setup}
voxel_betas <- matrix(NA, n_trials, n_vox)
```

For each voxel, build a tailored design matrix by weighting the three basis columns with that voxel's estimated HRF coefficients, then run LSS.

```{r voxel-lss-loop}
for (v in 1:n_vox) {
  Xv <- X_can * 0
  for (tr in 1:n_trials) {
    rset_tr <- fmrihrf::regressor_set(onsets = onsets[tr], fac = factor(1),
      hrf = HRF_SPMG3, duration = 0, span = 30, summate = FALSE)
    cols <- as.matrix(fmrihrf::evaluate(rset_tr, grid = grid, precision = 0.1, method = "conv"))
    Xv[, tr] <- cols %*% hrf_weights[, v]
  }
  voxel_betas[, v] <- lss(Y[, v, drop = FALSE], Xv, method = "r_optimized")
}
```

For production analyses, `lss_with_hrf()` wraps this loop with optional C++ acceleration. See `?lss_with_hrf`.

## OASIS Alternative

OASIS handles multi-basis HRFs and LSS in a single call. You pass `HRF_SPMG3` and it solves for all basis components simultaneously, with optional ridge regularization for stability.

```{r oasis-method}
oasis_betas <- lss(
  Y, X = NULL, method = "oasis",
  oasis = list(
    design_spec = list(
      sframe = sframe,
      cond = list(onsets = onsets, hrf = HRF_SPMG3, span = 30)
    ),
    ridge_mode = "fractional", ridge_x = 0.01, ridge_b = 0.01
  )
)
oasis_canonical <- oasis_betas[seq(1, nrow(oasis_betas), by = 3), ]
```

See `vignette("oasis_method")` for ridge tuning and standard-error computation.

## Compare Methods

### Compute accuracy metrics

```{r compare-metrics}
cors  <- c(Standard  = cor(as.vector(standard_betas), as.vector(true_betas)),
           VoxelHRF  = cor(as.vector(voxel_betas),    as.vector(true_betas)),
           OASIS     = cor(as.vector(oasis_canonical), as.vector(true_betas)))
rmses <- c(Standard  = sqrt(mean((standard_betas - true_betas)^2)),
           VoxelHRF  = sqrt(mean((voxel_betas    - true_betas)^2)),
           OASIS     = sqrt(mean((oasis_canonical - true_betas)^2)))
comparison_summary <- data.frame(
  Correlation = round(cors, 3),
  RMSE = round(rmses, 3)
)
comparison_summary
```

```{r compare-metrics-check, include = FALSE}
stopifnot(all(is.finite(cors)), all(is.finite(rmses)))
```

These metrics give you a direct accuracy check on the simulated data rather than relying on a visual impression alone.

### Scatter plots

Points closer to the diagonal mean better recovery.

```{r compare-setup}
rng <- range(true_betas)
cls <- rep(1:n_vox, each = n_trials)
```

Each panel shows estimated versus true betas; colour distinguishes voxels.

```{r compare-plots, fig.width=9, fig.height=3.5, fig.alt="Scatter plots comparing true vs estimated betas for each method."}
par(mfrow = c(1, 3), mar = c(4, 4, 3, 1))
for (i in seq_along(cors)) {
  est <- list(standard_betas, voxel_betas, oasis_canonical)[[i]]
  plot(true_betas, est, pch = 19, col = cls, xlim = rng, ylim = rng,
       xlab = "True", ylab = "Estimated",
       main = paste0(names(cors)[i], " (r=", round(cors[i], 2), ")"))
  abline(0, 1, lty = 2, col = "gray")
}
```

## Next Steps

**When to use voxel-wise HRFs.** You should consider this approach when regions differ in vascular architecture (motor vs. visual cortex), when studying populations with altered neurovascular coupling (aging, clinical), or when high-resolution acquisitions expose laminar-level variation.

**Choosing a method.** Standard LSS with a canonical HRF is sufficient for homogeneous responses. Voxel-wise HRF LSS improves accuracy when HRF heterogeneity is expected. OASIS is preferred for rapid-event designs or when you want a single-step solve with HRF flexibility and ridge control.

**Further reading:**

- `vignette("fmrilss")` -- LSS basics
- `vignette("oasis_method")` -- OASIS solver with ridge, multi-basis HRFs, and standard errors
- `vignette("sbhm")` -- Shared-Basis HRF Matching for efficient voxel-specific HRFs
- `?estimate_voxel_hrf` and `?lss_with_hrf` -- production-ready voxel-wise HRF workflow