---
title: Post-hoc corrections and spatial FDR
output:
  rmarkdown::html_vignette:
    toc: yes
    toc_depth: 2
    css: albers.css
    includes:
      in_header: albers-header.html
params:
  family: red
  preset: interaction
resource_files:
- albers.css
- albers.js
- albers-header.html
- fonts

vignette: |
  %\VignetteIndexEntry{Post-hoc corrections and spatial FDR}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
if (requireNamespace("ragg", quietly = TRUE)) knitr::opts_chunk$set(dev = "ragg_png")
if (requireNamespace("albersdown", quietly = TRUE)) {
  if (requireNamespace("systemfonts", quietly = TRUE)) try(albersdown::albers_register_fonts(), silent = TRUE)
  if (requireNamespace("ggplot2", quietly = TRUE)) try(ggplot2::theme_set(albersdown::theme_albers(family = params$family, preset = params$preset)), silent = TRUE)
}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  message = FALSE,
  warning = FALSE,
  fig.width = 6,
  fig.height = 4
)
library(fmrigds)
set.seed(1)
```

```{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-interaction","preset-study","preset-structural","preset-adobe","preset-midnight");',
    'document.body.classList.add("palette-%s","preset-%s");',
    '});</script>'
  ),
  params$family,
  params$preset
))
```

## When do you need this?

After running a group-level analysis, you typically have one p-value per
spatial unit (voxel, parcel, or ROI). Testing thousands of locations
simultaneously inflates false positives, so you need a multiple comparison
correction. Standard FDR (Benjamini--Hochberg) treats every test
independently. **Spatial FDR** goes further by grouping spatially related
tests---for instance, voxels within the same parcel---and adjusting within
and across groups to better control the false discovery rate in structured
neuroimaging data.

This vignette shows how to apply both standard and spatial FDR corrections in
fmrigds, including how to chain operations on realized GDS objects.

## Sample data

We'll use a small dataset with 6 ROIs, 4 subjects, and 1 contrast. ROIs 1--3
represent a "frontal" region with a genuine effect; ROIs 4--6 represent an
"occipital" region with little activation.

```{r sample-data, include = FALSE}
# Build a small CSV with known signal structure
tmp <- tempfile(fileext = ".csv")
lines <- "sample,subject,contrast,beta,var"
rois <- paste0("ROI_", 1:6)
subs <- paste0("s", 1:4)

for (s in subs) {
  for (r in seq_along(rois)) {
    mu <- if (r <= 3) 2.0 else 0.1
    b <- round(mu + rnorm(1, 0, 0.3), 3)
    lines <- c(lines, paste(rois[r], s, "c1", b, 0.04, sep = ","))
  }
}

writeLines(lines, tmp)
```

```{r load-data}
g0 <- compute(gds(tmp))
explain(g0)
```

```{r data-overview, echo = FALSE, fig.cap = "Per-subject betas for 6 ROIs. ROIs 1--3 (frontal) have a true effect around 2; ROIs 4--6 (occipital) are near zero.", fig.height = 3.5}
mat <- assay(g0, "beta")[, , 1]
cols <- rep(c("#5B8DBE", "#D4694A"), each = 3)
par(mar = c(5, 4, 3, 1))
matplot(t(mat), type = "b", pch = 19, lty = 1, col = cols,
        xlab = "Subject", ylab = "Beta",
        main = "Effect estimates by ROI", xaxt = "n")
axis(1, at = 1:4, labels = paste0("s", 1:4))
legend("topright", legend = c("Frontal (ROI 1-3)", "Occipital (ROI 4-6)"),
       col = c("#5B8DBE", "#D4694A"), pch = 19, lty = 1, bty = "n", cex = 0.85)
```

The frontal ROIs cluster around 2 while the occipital ROIs hover near zero,
just as we designed. A group model should find the frontal effect significant
and the occipital effect not.

## Chaining operations on realized objects

All fmrigds verbs accept realized GDS objects, not just plans. Internally,
`as_plan()` wraps the in-memory data so you can keep chaining without going
back to disk:

```{r realized-verbs}
# Reduce directly on a realized GDS
g_fe <- reduce(g0, method = "fixed") |> compute()
names(assays(g_fe))

# Chain reduce and post-hoc in one pipeline
g_bh <- reduce(g0, method = "fixed") |>
  posthoc("fdr:bh") |>
  compute()

names(assays(g_bh))
```

The group result now includes FDR-adjusted q-values alongside the standard
meta-analytic outputs.

## Standard FDR correction

Benjamini--Hochberg (`"fdr:bh"`) controls the expected proportion of false
discoveries. It treats each sample location independently:

```{r standard-fdr}
assay(g_bh, "p")
assay(g_bh, "q")
```

```{r fdr-comparison, echo = FALSE, fig.cap = "Raw p-values (left) and FDR-adjusted q-values (right). The dashed line marks 0.05.", fig.height = 3.5}
p_vals <- drop(assay(g_bh, "p"))
q_vals <- drop(assay(g_bh, "q"))
nms <- rownames(assay(g_bh, "p"))
if (is.null(nms)) nms <- paste0("ROI_", seq_along(p_vals))
cols <- ifelse(seq_along(p_vals) <= 3, "#5B8DBE", "#D4694A")

par(mfrow = c(1, 2), mar = c(5, 4, 3, 1))

bp1 <- barplot(-log10(p_vals), names.arg = nms, col = cols, border = NA,
        ylab = expression(-log[10](p)), main = "Raw p-values", las = 2,
        cex.names = 0.8)
abline(h = -log10(0.05), lty = 2, col = "grey40")

bp2 <- barplot(-log10(q_vals), names.arg = nms, col = cols, border = NA,
        ylab = expression(-log[10](q)), main = "BH-adjusted q-values", las = 2,
        cex.names = 0.8)
abline(h = -log10(0.05), lty = 2, col = "grey40")
```

ROIs with strong effects (frontal, blue bars) remain significant after
correction---their bars extend well above the 0.05 threshold (dashed line).
The occipital ROIs (orange) have much weaker evidence and fall below the
threshold after adjustment.

## Spatial FDR

When your samples have spatial structure---parcels belonging to brain regions,
or voxels within an atlas---`"fdr:spatial"` can exploit that grouping for a
more powerful correction. The algorithm:

1. **Aggregates** p-values within each group using Simes' method
2. **Applies** weighted Benjamini--Hochberg across groups (weights proportional
   to group size)
3. **Broadcasts** the group-level q-values back to individual samples

You supply the grouping via `options$group`, a vector with one entry per
sample. Here we assign ROIs 1--3 to region 1 ("frontal") and ROIs 4--6 to
region 2 ("occipital"):

```{r spatial-fdr}
region <- c(1, 1, 1, 2, 2, 2)

g_sfdr <- posthoc(g_fe, method = "fdr:spatial",
                  options = list(group = region)) |>
  compute()

assay(g_sfdr, "q")
```

Samples within the same region share a group-level q-value. This is useful
when you expect effects to cluster spatially---as they typically do in
neuroimaging.

```{r spatial-vs-standard, echo = FALSE, fig.cap = "Standard BH correction (left) versus spatial FDR (right). Spatial FDR assigns a single q-value per region.", fig.height = 3.5}
q_bh <- drop(assay(g_bh, "q"))
q_sp <- drop(assay(g_sfdr, "q"))

par(mfrow = c(1, 2), mar = c(5, 4, 3, 1))
bp1 <- barplot(-log10(q_bh), names.arg = nms, col = cols, border = NA,
        ylab = expression(-log[10](q)), main = "Standard BH", las = 2,
        cex.names = 0.8)
abline(h = -log10(0.05), lty = 2, col = "grey40")

bp2 <- barplot(-log10(q_sp), names.arg = nms, col = cols, border = NA,
        ylab = expression(-log[10](q)), main = "Spatial FDR", las = 2,
        cex.names = 0.8)
abline(h = -log10(0.05), lty = 2, col = "grey40")
```

Notice that within each region, the spatial FDR q-values are identical---the
algorithm operates on groups, not individual tests. Spatial FDR can be more
powerful than standard BH when the spatial grouping reflects real biological
structure, because combining evidence across neighbouring locations
strengthens the signal.

### Requirements for spatial FDR

- The GDS must contain `p` values (or `z`/`t` from which p can be derived).
  Running `reduce()` first produces these automatically.
- `options$group` must be a vector of length equal to the number of samples
  (the first dimension of the assay arrays). Values can be numeric, character,
  or factor. `NA` entries are excluded from grouping.

## Next steps

- `vignette("fmrigds")` --- getting started with the core pipeline
- `vignette("spatial-operations")` --- masking, alignment, and space
  transformations
- `vignette("fmristore-ingestion")` --- reading fmristore HDF5 files
- `?posthoc` for the full list of correction methods

```{r cleanup, include = FALSE}
unlink(tmp)
```