---
title: Using fmridesign with fmrilss
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{Using fmridesign with fmrilss}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup-opts, include = FALSE}
has_fmridesign <- requireNamespace("fmridesign", quietly = TRUE)
if (has_fmridesign &&
    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 = "#>",
  eval = has_fmridesign,
  fig.width = 7,
  fig.height = 5,
  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
))
```

```{r setup}
library(fmrilss)
library(fmridesign)
library(fmrihrf)
```

## Introduction

The `fmridesign` package provides a powerful formula-based interface for creating fMRI design matrices. While `fmrilss` has its own `design_spec` format for the OASIS method, you can now use `event_model` and `baseline_model` objects directly with the new `lss_design()` function.

### When to Use lss_design()

Use `lss_design()` when:

- You have **multi-condition factorial designs**
- You need **parametric modulators** (e.g., RT, difficulty ratings)
- You want **design validation** and diagnostic tools
- You prefer **formula-based specification**
- You need **explicit multi-run handling** with run-relative onsets

### When to Use Traditional lss()

Use the traditional `lss()` interface when:

- You have simple **trial-wise designs** already constructed
- You want **minimal dependencies**
- You're using **internal SBHM pipelines**
- You already have design matrices prepared

## Simulated Data Setup

We'll build a small two-run experiment from scratch so every piece of the workflow is visible. Start with the run-level parameters:

```{r sim-params}
set.seed(123)
n_scans_per_run <- 150
n_runs          <- 2
n_voxels        <- 500
TR              <- 2
```

Next, a trial table with run-relative onsets — six trials per run — and the matching sampling frame:

```{r sim-trials}
trials <- data.frame(
  onset = rep(c(10, 30, 50, 70, 90, 110), times = 2),
  run   = rep(1:2, each = 6)
)
sframe <- sampling_frame(blocklens = rep(n_scans_per_run, n_runs), TR = TR)
```

Build a trial-wise event model, then extract the design matrix used to generate the signal:

```{r sim-design}
emod_sim <- event_model(
  onset ~ trialwise(basis = "spmg1"),
  data = trials, block = ~run, sampling_frame = sframe
)
X_trial <- as.matrix(design_matrix(emod_sim))
```

A per-run intercept + linear drift serves as the baseline:

```{r sim-baseline}
bmodel_sim <- baseline_model(
  basis = "poly", degree = 1, sframe = sframe, intercept = "runwise"
)
Z_baseline <- as.matrix(design_matrix(bmodel_sim))
```

Finally, combine true trial effects, baseline drift, and noise to produce `Y`:

```{r sim-data}
true_betas <- matrix(rnorm(12 * n_voxels, mean = 1.5, sd = 0.8), 12, n_voxels)
Y <- X_trial  %*% true_betas +
     Z_baseline %*% matrix(rnorm(ncol(Z_baseline) * n_voxels, sd = 2), ncol = n_voxels) +
     matrix(rnorm(nrow(X_trial) * n_voxels, sd = 3), nrow(X_trial), n_voxels)
dim(Y)
```

## Quick Start

Here's a minimal example using `lss_design()`:

```{r quickstart}
# Use the first run only for quick start
trials_run1 <- trials[trials$run == 1, ]
sframe_run1 <- sampling_frame(blocklens = n_scans_per_run, TR = TR)
Y_run1 <- Y[1:n_scans_per_run, ]

# Build event model with trialwise design
emod <- event_model(
  onset ~ trialwise(basis = "spmg1"),
  data = trials_run1,
  block = ~run,
  sampling_frame = sframe_run1
)

# Run LSS
beta <- lss_design(Y_run1, emod, method = "oasis")

# Result: 6 trials × 500 voxels
dim(beta)
```

```{r quickstart-check, include = FALSE}
stopifnot(all(dim(beta) == c(6, n_voxels)), all(is.finite(beta)))
```

## Multi-Run Experiments

One of the key advantages of `lss_design()` is automatic handling of multi-run experiments with proper onset timing.

### Onset Convention: Run-Relative

**Important:** When using `fmridesign`, onsets should be **run-relative** (resetting to 0 at the start of each run). This is the standard convention for multi-run fMRI experiments.

```{r multirun}
# Trial data with run-relative onsets (already defined above)
print(trials)

# Create event model - conversion to global time is automatic
emod_multi <- event_model(
  onset ~ trialwise(basis = "spmg1"),
  data = trials,
  block = ~run,
  sampling_frame = sframe
)

# Run LSS
beta_multi <- lss_design(Y, emod_multi, method = "oasis")

# Result: 12 trials × 500 voxels
dim(beta_multi)
```

### Why Run-Relative Onsets?

- **Standard convention**: Most experiment software (E-Prime, PsychoPy) logs onsets relative to run start
- **Easier data management**: No manual offset calculations needed
- **Automatic conversion**: `fmridesign` handles conversion to global time internally
- **Less error-prone**: Reduces risk of incorrect timing specifications

## Adding Baseline Correction

The `baseline_model` allows you to specify drift correction, block intercepts, and nuisance regressors in a structured way.

```{r baseline}
# Create baseline model with B-spline drift correction
bmodel <- baseline_model(
  basis = "bs",
  degree = 5,
  sframe = sframe,
  intercept = "runwise"
)

# LSS with baseline correction
beta_baseline <- lss_design(Y, emod_multi, bmodel, method = "oasis")

dim(beta_baseline)
```

### Adding Motion Parameters

For demonstration, we'll create synthetic motion regressors.

```{r motion-demo}
# Simulate motion parameters (6 motion parameters per run)
motion_run1 <- matrix(rnorm(n_scans_per_run * 6, sd = 0.5),
                      nrow = n_scans_per_run, ncol = 6)
motion_run2 <- matrix(rnorm(n_scans_per_run * 6, sd = 0.5),
                      nrow = n_scans_per_run, ncol = 6)

# Create baseline model with motion as nuisance
bmodel_motion <- baseline_model(
  basis = "bs",
  degree = 5,
  sframe = sframe,
  intercept = "runwise",
  nuisance_list = list(motion_run1, motion_run2)
)

beta_motion <- lss_design(Y, emod_multi, bmodel_motion, method = "oasis")
dim(beta_motion)
```

```{r motion-real, eval = FALSE}
# In real analysis, load motion from files:
motion_run1 <- as.matrix(read.table("motion_run1.txt"))
motion_run2 <- as.matrix(read.table("motion_run2.txt"))

bmodel <- baseline_model(
  basis = "bs",
  degree = 5,
  sframe = sframe,
  intercept = "runwise",
  nuisance_list = list(motion_run1, motion_run2)
)
```

## Multi-Basis HRFs

For multi-basis HRF models (e.g., canonical + temporal + dispersion derivatives), use `nbasis`:

```{r multibasis}
# Create event model with SPMG3 (3 basis functions)
emod_3basis <- event_model(
  onset ~ trialwise(basis = "spmg3", nbasis = 3),
  data = trials,
  block = ~run,
  sampling_frame = sframe
)

# LSS will auto-detect K = 3
beta_3basis <- lss_design(Y, emod_3basis, method = "oasis")

# Output: (12 trials × 3 basis) × 500 voxels = 36 × 500
dim(beta_3basis)

# Extract canonical basis estimates (every 3rd row starting at 1)
beta_canonical <- beta_3basis[seq(1, nrow(beta_3basis), by = 3), ]
dim(beta_canonical)
```

## Ridge Regularization

For designs with potential collinearity, use ridge regularization:

```{r ridge}
beta_ridge <- lss_design(
  Y, emod_multi, bmodel,
  method = "oasis",
  oasis = list(
    ridge_mode = "fractional",
    ridge_x = 0.02,
    ridge_b = 0.02
  )
)

dim(beta_ridge)
```

## Comparison with design_spec

### Using design_spec (old approach)

```{r designspec, eval = FALSE}
# Manual design_spec construction requires global/absolute onsets
# For run-relative onsets (10, 30, 50, 70, 90, 110) in each of 2 runs,
# global onsets would be: 10, 30, 50, 70, 90, 110, 310, 330, 350, 370, 390, 410
# (second run starts at 150 scans × 2s = 300s)

beta_old <- lss(
  Y,
  method = "oasis",
  oasis = list(
    design_spec = list(
      sframe = sframe,
      cond = list(
        onsets = c(10, 30, 50, 70, 90, 110, 310, 330, 350, 370, 390, 410),
        hrf = HRF_SPMG1,
        span = 30
      )
    )
  )
)
```

**Limitations:**
- Requires **global/absolute onsets** (not run-relative)
- No structured baseline handling
- No multi-condition support
- No parametric modulators
- Less validation

### Using lss_design() (new approach)

```{r lssdesign}
# Formula-based design with run-relative onsets
emod_new <- event_model(
  onset ~ trialwise(basis = "spmg1"),
  data = trials,
  block = ~run,
  sampling_frame = sframe
)

bmodel_new <- baseline_model(basis = "bs", degree = 5, sframe = sframe)

beta_new <- lss_design(Y, emod_new, bmodel_new, method = "oasis")
```

**Advantages:**
- **Run-relative onsets** (standard convention)
- **Structured baseline** (drift + nuisance)
- **Automatic validation**
- **Richer metadata**
- **Formula-based DSL**

## Validation: Estimated vs True Betas

Let's visualize how well LSS recovers the true trial effects from our simulated data.

```{r validation-plot, fig.alt = "Scatter plot showing estimated LSS betas versus true simulated betas with near-diagonal alignment"}
# Compare estimated betas (with baseline correction) to true betas
# Flatten matrices for plotting
est_vec <- as.vector(beta_baseline)
true_vec <- as.vector(true_betas)

recovery_summary <- data.frame(
  Correlation = cor(est_vec, true_vec),
  RMSE = sqrt(mean((est_vec - true_vec)^2))
)
recovery_summary

# Plot
plot(true_vec, est_vec,
     pch = 16, cex = 0.3, col = rgb(0, 0, 0, 0.3),
     xlab = "True Beta", ylab = "Estimated Beta",
     main = sprintf("LSS Recovery (r = %.3f)", recovery_summary$Correlation))
abline(0, 1, col = "red", lwd = 2, lty = 2)
grid()
```

```{r validation-check, include = FALSE}
stopifnot(all(is.finite(as.matrix(recovery_summary))))
```

The fit is not perfect because the simulation adds substantial baseline structure and noise, but the recovered betas still track the ground truth in a way that is easy to diagnose numerically and visually.

## Advanced: Parametric Modulators

`event_model` supports parametric modulators for trial-by-trial amplitude modulation. This example demonstrates the syntax but requires special setup.

```{r parametric, eval = FALSE}
# Trial data with reaction times
trials_rt <- data.frame(
  onset = c(10, 30, 50, 70, 90, 110),
  RT = c(0.5, 0.7, 0.6, 0.8, 0.5, 0.9),
  run = 1
)

# Center RT
trials_rt$RT_c <- scale(trials_rt$RT, center = TRUE, scale = FALSE)[, 1]

# Model: trial effects + RT modulation
emod_rt <- event_model(
  onset ~ trialwise() + hrf(RT_c),
  data = trials_rt,
  block = ~run,
  sampling_frame = sframe
)

# Note: This creates trial-wise regressors PLUS an RT amplitude modulator
# May require special handling in OASIS for proper separation
```

## Troubleshooting

### Error: "Y has X rows but sampling_frame expects Y scans"

**Cause:** Mismatch between data dimensions and sampling_frame specification.

**Solution:** Check that `sum(blocklens)` matches `nrow(Y)`:

```{r troubleshoot1}
sframe_check <- sampling_frame(blocklens = c(150, 150), TR = 2)
sum(fmrihrf::blocklens(sframe_check))
nrow(Y)
```

### Error: "event_model and baseline_model have different sampling_frames"

**Cause:** The two models were created with different `sampling_frame` objects.

**Solution:** Use the same `sframe` object for both:

```{r troubleshoot2, eval = FALSE}
sframe <- sampling_frame(blocklens = c(150, 150), TR = 2)

emod <- event_model(..., sampling_frame = sframe)
bmodel <- baseline_model(..., sframe = sframe)
```

### Warning: "High collinearity detected"

**Cause:** Events are too close together or design is ill-conditioned.

**Solution:** Use ridge regularization:

```{r troubleshoot3, eval = FALSE}
beta <- lss_design(
  Y, emod, bmodel,
  oasis = list(ridge_mode = "fractional", ridge_x = 0.02)
)
```

## Summary

The `lss_design()` function provides a modern, formula-based interface for LSS analysis that:

- Handles **multi-run experiments** correctly with run-relative onsets
- Provides **structured baseline** specification
- Supports **multi-basis HRFs** with automatic detection
- Validates **design compatibility** automatically
- Integrates with the **fmridesign ecosystem**

For simple designs or when you already have design matrices prepared, the traditional `lss()` interface remains fully supported and unchanged.

## Further Reading

- `vignette("fmrilss")` - Traditional LSS interface
- `vignette("oasis_method")` - OASIS method details
- `vignette("a_04_event_models", package = "fmridesign")` - Event model tutorial
- `vignette("a_03_baseline_model", package = "fmridesign")` - Baseline model tutorial