---
title: Running one model across many subjects
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{Running one model across many subjects}
  %\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 = "#>", message = FALSE, warning = TRUE)
```

```{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 load, message=FALSE}
library(fmrireg)
```

## The problem

A common study runs the **same model** across many subjects. The design is
constant -- the same conditions, HRF, baseline, contrasts, and fitting options --
but the data changes per subject: different files, different event order,
different run lengths. Looping `fmri_lm()` by hand is tedious and does not lend
itself to parallel or cluster execution.

`fmrireg` separates the part that is **invariant** (an `fmri_template`) from the
part that **varies per subject** (a *binding*: that subject's data). Combining
the two produces a serializable `fmri_job`. With file-backed bindings, the job
stores paths rather than voxel data and is small enough to ship to a worker. An
inline matrix binding embeds that matrix in the job; this is convenient for the
runnable example below but is not the memory-efficient choice for a real study.

The four pieces:

| Piece | What it is |
|---|---|
| `fmri_template()` | the invariant model: formula, baseline spec, contrasts, fit control, reducer |
| a *binding* | one subject's data: scans, run lengths, TR, events, confounds |
| `fmri_job` | template + one binding -- a serializable recipe (built by `instantiate()`) |
| a *reducer* | what each fitted subject becomes (a tidy table, or files on disk) |

This closes a loop: BIDS in (via `bidser`) → per-subject fit → BIDS-keyed maps
out (`write_results`) → group level (`collect_results` → `fmri_meta`).

## A minimal, runnable example

Define the template once, describe two subjects as inline bindings, then fan
the model out. The resulting jobs contain these small matrices. Real studies
normally use file paths so jobs remain lightweight and workers open the imaging
data lazily.

```{r minimal}
# one subject's data as a binding
make_subject <- function(id, seed) {
  set.seed(seed)
  runs <- c(40L, 40L); TR <- 2
  Y <- matrix(rnorm(sum(runs) * 5), sum(runs), 5)        # timepoints x voxels
  ev <- do.call(rbind, lapply(1:2, function(r)
    data.frame(onset = seq(4, 60, by = 8),
               condition = factor(rep(c("A", "B"), length.out = 8)),
               run = r)))
  list(id = id, scans = Y, TR = TR, run_length = runs, events = ev)
}
subjects <- list(make_subject("sub-01", 1), make_subject("sub-02", 2))

# the invariant model, defined once
tmpl <- fmri_template(
  onset ~ hrf(condition), ~ run,
  baseline = baseline_spec(degree = 3),
  reducer  = reduce_betas()        # each subject -> a tidy beta table
)

# Bind data to jobs, validate locally, and only then fan out.
jobs <- instantiate(tmpl, subjects)
flight <- preflight(jobs)
stopifnot(flight$ok, flight$n_jobs == 2L)

res  <- run_jobs(jobs)
res

values <- batch_values(res)        # named by job id
stopifnot(
  length(batch_errors(res)) == 0L,
  identical(names(values), c("sub-01", "sub-02")),
  all(vapply(values, nrow, integer(1)) == 10L)
)
head(values[["sub-01"]])
```

`run_jobs()` isolates per-job failures: a subject that errors is recorded in the
result rather than aborting the batch (see `batch_errors()`).

## Reducers: what each subject becomes

A reducer runs on the worker, right after the fit, so only the reduced output
crosses the worker→driver (or worker→disk) boundary -- not a whole fitted model.

| Reducer | Output |
|---|---|
| `reduce_identity()` | the entire `fmri_lm` object (largest) |
| `reduce_betas()` | tidy data frame: `job_id, term, voxel, estimate, se, stat` |
| `reduce_contrasts()` | tidy data frame of fitted contrasts |
| `reduce_write_results()` | writes BIDS-keyed maps to disk, returns the paths |

You can also pass your own `function(fit, job)`; keep it a top-level/package
function so it serializes to a worker.

## Discovering subjects from BIDS

For BIDS-formatted data, `from_bids()` (which uses the `bidser` package) builds
the per-subject bindings -- scans, events, confounds, TR, and run lengths -- so
you do not assemble them by hand.

```{r from-bids, eval=FALSE}
proj <- bidser::bids_project("study/", fmriprep = TRUE)
mani <- from_bids(
  proj, task = "stroop", space = "MNI152NLin2009cAsym",
  confounds = bidser::confound_set("motion6"),
  mask = "study/derivatives/.../space-MNI152NLin2009cAsym_desc-brain_mask.nii.gz"
)
jobs <- instantiate(tmpl, mani)
```

The design formula's variables must match the columns of the BIDS
`events.tsv` (e.g. `trial_type`); a `run` column is added for the block
structure.

## Preflight before you fan out

`preflight()` validates jobs on the driver -- design columns present, TR and run
lengths consistent, confound dimensions correct -- so problems surface in
seconds rather than on a compute node after the queue drains.
The runnable workflow called it immediately after `instantiate()` and asserted
`flight$ok` before `run_jobs()`. Keep that order in production code.

## Running: locally, in parallel, or on a cluster

Sequential is the default. For parallelism, set a `future` plan; `run_jobs()`
dispatches through it -- including `future.batchtools` cluster plans -- with **no
scheduler-specific code** in `fmrireg`.

```{r parallel, eval=FALSE}
library(future)
plan(multisession, workers = 4)        # or future.batchtools::batchtools_slurm
res <- run_jobs(jobs, parallel = TRUE)
```

For an array scheduler, export the jobs to disk and drive them with whatever
system you have:

```{r export, eval=FALSE}
export_jobs(jobs, "study/jobs")
# study/jobs/manifest.rds + a backend-agnostic run_one.R
# then, per array task:
#   SLURM: Rscript study/jobs/run_one.R $SLURM_ARRAY_TASK_ID
#   local: for i in $(seq 1 N); do Rscript study/jobs/run_one.R $i; done
```

`run_one.R` reads the manifest, reconstructs job `i`, fits, and writes its
reduced output. Each worker needs `fmrireg`, packages referenced by the template
or custom reducer, access to every file-backed input path, and write access to
the chosen output location. A serialized recipe does not make local paths or
package dependencies portable by itself.

## Closing the loop: the group level

The inline example can close the statistical loop immediately. Here we collect
the condition-B estimate and SE for teaching voxel 1 from each completed job,
then compute their inverse-variance fixed-effect mean. With only two simulated
subjects this is a pipeline check, not a scientifically powered group study.

```{r executable-group-handoff}
reduced <- do.call(rbind, values)
teaching_effects <- subset(
  reduced,
  term == "condition_condition.B" & voxel == 1L,
  select = c(job_id, estimate, se)
)
teaching_effects$roi <- "teaching_voxel_1"
teaching_effects$contrast <- "condition_B"

gd_inline <- fmrireg::group_data(
  teaching_effects,
  format = "csv",
  effect_cols = c(beta = "estimate", se = "se"),
  subject_col = "job_id",
  roi_col = "roi",
  contrast_col = "contrast"
)
fm_inline <- fmri_meta(gd_inline, ~ 1, method = "fe", verbose = FALSE)
group_summary <- data.frame(
  target = "Condition B at teaching voxel 1",
  subjects = nrow(teaching_effects),
  estimate = as.numeric(fm_inline$coefficients[1, "(Intercept)"]),
  std_error = as.numeric(fm_inline$se[1, "(Intercept)"])
)
stopifnot(
  group_summary$subjects == 2L,
  all(is.finite(as.matrix(group_summary[c("estimate", "std_error")])))
)
knitr::kable(group_summary, digits = 3,
             caption = "Executable subject-to-group handoff")
```

With `reduce_write_results()`, each worker writes BIDS-keyed statistical maps.
`collect_results()` gathers them back into a `group_data` object that
`fmri_meta()` meta-analyses. (Writing both `beta` and `se` requires a contrast in
the model, so the group step has per-subject variance.)

```{r group, eval=FALSE}
con <- contrast_set(pair_contrast(~ trial_type == "incongruent",
                                  ~ trial_type == "congruent",
                                  name = "incong_gt_cong"))
tmpl <- fmri_template(
  onset ~ hrf(trial_type, contrasts = con), ~ run,
  baseline = baseline_spec(degree = 3, confounds = bidser::confound_set("motion6")),
  reducer  = reduce_write_results(format = "nifti", stats = c("beta", "se"),
                                  path = "study/glm")
)
run_jobs(instantiate(tmpl, mani))

gd <- collect_results("study/glm", space = "MNI152NLin2009cAsym")
fm <- fmri_meta(gd, ~ 1, method = "fe")
```

## Summary

- Define the model once with `fmri_template()`; vary only the data.
- `instantiate()` produces serializable `fmri_job` recipes; `from_bids()`
  populates them from a BIDS dataset.
- `preflight()` catches problems before fan-out.
- `run_jobs()` runs sequentially, in parallel via `future`, or through a custom
  backend; `export_jobs()` emits a scheduler-agnostic runner for array jobs.
- Reducers keep per-subject output compact; `collect_results()` →
  `fmri_meta()` close the loop at the group level.

## Next

- `vignette("group_analysis", package = "fmrireg")` — Combine subject-level estimates
- `vignette("functional_connectivity", package = "fmrireg")` — Use a time series as a regressor