---
title: Functional Connectivity (Seed-Based)
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{Functional Connectivity (Seed-Based)}
  %\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
))
```

## Overview

Seed-based functional connectivity analysis is a fundamental technique in fMRI research for identifying brain regions that show correlated activity with a specific region of interest (the "seed"). This vignette demonstrates how to perform seed-based connectivity analysis using fmrireg's flexible design matrix and GLM framework.

The approach treats the seed time series as a covariate rather than computing
simple correlations. This lets the model adjust for scanner drift and other
nuisance variables while estimating seed coupling. Because the simulated
noise is autocorrelated, we use an AR(1) temporal model. Because we test the
seed coefficient at every voxel, discoveries are defined by Benjamini--Hochberg
(BH) false-discovery-rate correction rather than uncorrected voxelwise p-values.

To illustrate the method clearly, we'll work with simulated data where we know the ground truth. We'll create a dataset with a hidden network of voxels that share a common signal with our seed region, then recover this network using our connectivity analysis.

## Creating a Test Dataset with Known Connectivity

First, we generate the data directly from the model we intend to fit. Every
voxel has independent AR(1) errors with the same `rho = 0.3`; there is no
unmodeled task signal. This matters: calibration cannot be assessed when the
generator contains mean structure that the fitted model silently omits.

```{r simulate-baseline-data}
set.seed(42)

Tlen <- 180L
TR <- 2
V <- 256L
true_rho <- 0.3

simulate_ar1_errors <- function(n_time, n_voxels, rho,
                                innovation_sd = 1) {
  innovations <- matrix(
    rnorm(n_time * n_voxels, sd = innovation_sd),
    nrow = n_time, ncol = n_voxels
  )
  errors <- matrix(0, nrow = n_time, ncol = n_voxels)
  errors[1, ] <- innovations[1, ] / sqrt(1 - rho^2)
  for (time_index in 2:n_time) {
    errors[time_index, ] <-
      rho * errors[time_index - 1L, ] + innovations[time_index, ]
  }
  errors
}

Y <- simulate_ar1_errors(Tlen, V, rho = true_rho)
dim(Y)
```

Now we'll create our ground-truth connectivity pattern. We generate an
external seed time series and inject it into a known subset of voxels. The
series is deliberately external: this simulation tests recovery against known
truth and is not pretending to perform seed extraction. A real analysis would
extract a summary series from a prespecified ROI.

```{r inject-network-signal}
# Generate the signal of interest and an independent negative-control probe.
seed_ts <- arima.sim(model = list(ar = 0.5), n = Tlen)
seed_ts <- as.numeric(base::scale(seed_ts))
null_probe <- arima.sim(model = list(ar = 0.4), n = Tlen)
null_probe <- as.numeric(base::scale(null_probe))

# Define which voxels belong to our injected network
anchor_voxel <- 10
net_idx <- c(anchor_voxel, sample(setdiff(seq_len(V), anchor_voxel), 40))

# Add the external seed signal to the network voxels
Y[, net_idx] <- Y[, net_idx] + 0.7 * seed_ts

# Rebuild the dataset from the modified matrix so the injected network signal
# is the data used by downstream model fitting.
dset_modified <- matrix_frame(
  Y,
  TR = TR,
  run_length = Tlen,
  event_table = data.frame(onset = 0, run = 1)
)
```

## Modeling Scanner Drift

Before we can accurately estimate connectivity, we need to account for low-frequency scanner drift that can create spurious correlations between voxels. The fmrireg package provides flexible tools for modeling these nuisance signals using basis functions.

```{r model-drift}
# Create a sampling frame for our single run
sframe <- sampling_frame(rep(Tlen, 1), TR = TR)

# Model drift using B-splines. This object enters the fitted model below.
bmodel <- baseline_model(basis = "bs", degree = 3, sframe = sframe)
```

## Connectivity Analysis Using fmrireg's GLM Framework

Instead of computing simple correlations, we treat the seed time series as a
GLM covariate. This estimates connectivity while simultaneously controlling
for confounds. `covariate()` adds a regressor without HRF convolution because
the seed signal is already sampled as a BOLD time series.

```{r fit-connectivity-glm}
# Set up the event model structure
# We need a minimal event_data frame to define the model structure
event_data <- data.frame(
  onset = samples(sframe)[1],  # Single onset to define the model
  run = 1                       # Single run indicator
)

# The null probe is included as a negative-control coefficient family: it is
# fitted at every voxel but was injected nowhere.
cov_data <- data.frame(
  seed = seed_ts,
  null_probe = null_probe
)

# Build the event model with both sampled time series as covariates.
emodel <- event_model(
  onset ~
    covariate(seed, data = cov_data, prefix = "seed") +
    covariate(null_probe, data = cov_data, prefix = "null"),
  data = event_data,
  block = ~ run,
  sampling_frame = sframe
)

# Combine event and baseline models with the dataset
fmodel <- fmri_model(emodel, bmodel, dset_modified)

# Fit the connectivity GLM across all voxels. The simulated errors are AR(1),
# so inference uses the corresponding global temporal-noise model.
fit <- fmri_lm(
  fmodel,
  dataset = dset_modified,
  control = fmri_lm_control(
    noise = noise_spec("ar1", pooling = "global")
  )
)
```

The seed-regressor statistic describes coupling after accounting for the
baseline model and the fitted AR(1) covariance. It is not yet a whole-brain
discovery decision; that requires multiplicity correction.

```{r extract-connectivity-stats}
# t-statistics of the model coefficients (coef() gives the estimates)
all_stats <- as.matrix(stats(fit, type = "betas"))
seed_cols <- grep("seed", colnames(all_stats), value = TRUE)
if (length(seed_cols) == 0) {
  stop("No seed estimate found in fitted model output")
}
seed_col_name <- seed_cols[1]
t_seed <- as.numeric(all_stats[, seed_col_name])

# Correct the seed-coefficient p-values across all fitted voxels.
all_pvals <- as.matrix(p_values(fit, type = "betas"))
p_seed <- as.numeric(all_pvals[, seed_col_name])
q_seed <- stats::p.adjust(p_seed, method = "BH")
discovered <- is.finite(q_seed) & q_seed < 0.05

# The negative-control coefficient is null at every voxel.
null_col_name <- grep("^null", colnames(all_pvals), value = TRUE)
stopifnot(length(null_col_name) == 1L)
q_null <- stats::p.adjust(
  as.numeric(all_pvals[, null_col_name]), method = "BH"
)
null_discovered <- is.finite(q_null) & q_null < 0.05
fitted_rho <- as.numeric(ar_parameters(fit, scope = "raw")[[1]])

# Check the distribution of our connectivity map
summary(t_seed)
```

## Validating the Results

Since the network membership is known, we can report both recovery and error.
Sensitivity is the fraction of injected network voxels discovered after BH
correction. Background false-positive rate (FPR) is the fraction of null
voxels called significant. False-discovery proportion (FDP) is the fraction
of all discoveries that are background voxels. A single FDP is a diagnostic,
not an estimate of FDR; aggregate calibration is assessed in the next section.

```{r validate-network-recovery}
mean_abs_t_network    <- mean(abs(t_seed[net_idx]), na.rm = TRUE)
mean_abs_t_background <- mean(abs(t_seed[-net_idx]), na.rm = TRUE)
sensitivity <- mean(discovered[net_idx])
background_fpr <- mean(discovered[-net_idx])
n_discoveries <- sum(discovered)
n_false_discoveries <- sum(discovered[-net_idx])
false_discovery_proportion <- n_false_discoveries / max(1L, n_discoveries)

top_ranked <- order(t_seed, decreasing = TRUE)[seq_along(net_idx)]
top_rank_enrichment <- mean(top_ranked %in% net_idx)

stopifnot(
  is.finite(mean_abs_t_network),
  is.finite(mean_abs_t_background),
  is.finite(sensitivity),
  is.finite(background_fpr),
  is.finite(false_discovery_proportion),
  is.finite(top_rank_enrichment),
  mean_abs_t_network > 5 * mean_abs_t_background,
  sensitivity >= 0.9,
  background_fpr <= 0.02,
  abs(fitted_rho - true_rho) < 0.12,
  top_rank_enrichment > 0.9
)

c(
  mean_abs_t_network = mean_abs_t_network,
  mean_abs_t_background = mean_abs_t_background,
  BH_discoveries = n_discoveries,
  sensitivity = sensitivity,
  background_FPR = background_fpr,
  false_discovery_proportion = false_discovery_proportion,
  complete_null_discoveries = sum(null_discovered),
  fitted_rho = fitted_rho,
  top_rank_enrichment = top_rank_enrichment
)
```

The effect-size separation, BH discoveries, sensitivity, FPR, and FDP answer
different questions. Reporting them together prevents a strong ranking result
from concealing null discoveries. The fitted AR coefficient also checks that
the temporal model is recovering the generator's `rho = 0.3`. None of these
single-run quantities, by itself, establishes repeated-sampling calibration.

## Repeated-sampling calibration

We therefore repeat a smaller, correctly specified simulation 100 times. Each
fit contains two coefficient families: the seed coefficient has a known
12-voxel network, while the independent negative-control coefficient is null
at all 64 voxels. BH correction is applied separately to each family.

```{r calibration-helpers, include=FALSE}
run_connectivity_fixture <- function(simulation_seed, n_time = 96L,
                                     n_voxels = 64L, n_active = 12L,
                                     effect = 0.8, rho = 0.3) {
  set.seed(simulation_seed)
  signal <- as.numeric(base::scale(arima.sim(
    model = list(ar = 0.5), n = n_time
  )))
  negative_control <- as.numeric(base::scale(arima.sim(
    model = list(ar = 0.4), n = n_time
  )))
  network <- sample(seq_len(n_voxels), n_active)
  fixture_data <- simulate_ar1_errors(n_time, n_voxels, rho)
  fixture_data[, network] <- fixture_data[, network] + effect * signal

  fixture_dataset <- matrix_frame(
    fixture_data, TR = 2, run_length = n_time,
    event_table = data.frame(onset = 0, run = 1L)
  )
  fixture_frame <- sampling_frame(n_time, TR = 2)
  fixture_covariates <- data.frame(
    seed = signal, null_probe = negative_control
  )
  fixture_event <- event_model(
    onset ~
      covariate(seed, data = fixture_covariates, prefix = "seed") +
      covariate(null_probe, data = fixture_covariates, prefix = "null"),
    data = data.frame(onset = samples(fixture_frame)[1], run = 1L),
    block = ~ run, sampling_frame = fixture_frame
  )
  fixture_baseline <- baseline_model(
    basis = "constant", sframe = fixture_frame
  )
  fixture_model <- fmri_model(
    fixture_event, fixture_baseline, fixture_dataset
  )
  fixture_fit <- fmri_lm(
    fixture_model, dataset = fixture_dataset,
    control = fmri_lm_control(
      noise = noise_spec("ar1", pooling = "global")
    )
  )

  fixture_p <- as.matrix(p_values(fixture_fit, type = "betas"))
  fixture_seed_col <- grep("^seed", colnames(fixture_p), value = TRUE)
  fixture_null_col <- grep("^null", colnames(fixture_p), value = TRUE)
  seed_discovery <- p.adjust(
    fixture_p[, fixture_seed_col], method = "BH"
  ) < 0.05
  null_discovery <- p.adjust(
    fixture_p[, fixture_null_col], method = "BH"
  ) < 0.05

  c(
    sensitivity = mean(seed_discovery[network]),
    background_fpr = mean(seed_discovery[-network]),
    fdp = sum(seed_discovery[-network]) / max(1L, sum(seed_discovery)),
    complete_null_rejection = as.numeric(any(null_discovery)),
    fitted_rho = as.numeric(
      ar_parameters(fixture_fit, scope = "raw")[[1]]
    )
  )
}

one_sided_mean_bounds <- function(x, level = 0.95) {
  critical <- qt(level, df = length(x) - 1L)
  standard_error <- sd(x) / sqrt(length(x))
  c(
    lower = mean(x) - critical * standard_error,
    upper = mean(x) + critical * standard_error
  )
}

one_sided_binomial_upper <- function(successes, trials, level = 0.95) {
  if (successes == trials) return(1)
  qbeta(level, successes + 1, trials - successes)
}
```

```{r calibrate-connectivity}
n_calibration <- 100L
calibration <- vapply(
  1000L + seq_len(n_calibration),
  run_connectivity_fixture,
  numeric(5)
)

fdp_bound <- one_sided_mean_bounds(calibration["fdp", ])
sensitivity_bound <- one_sided_mean_bounds(
  calibration["sensitivity", ]
)
background_fpr_bound <- one_sided_mean_bounds(
  calibration["background_fpr", ]
)
n_null_rejections <- sum(calibration["complete_null_rejection", ])
complete_null_upper <- one_sided_binomial_upper(
  n_null_rejections, n_calibration
)

calibration_summary <- data.frame(
  metric = c(
    "Mixed-family FDP",
    "Network sensitivity",
    "Background FPR",
    "Complete-null family rejection rate",
    "Fitted AR(1) coefficient (truth 0.3)"
  ),
  estimate = c(
    mean(calibration["fdp", ]),
    mean(calibration["sensitivity", ]),
    mean(calibration["background_fpr", ]),
    mean(calibration["complete_null_rejection", ]),
    mean(calibration["fitted_rho", ])
  ),
  one_sided_95_bound = c(
    fdp_bound[["upper"]],
    sensitivity_bound[["lower"]],
    background_fpr_bound[["upper"]],
    complete_null_upper,
    NA_real_
  ),
  bound_direction = c("upper", "lower", "upper", "upper", "-")
)

stopifnot(
  fdp_bound[["upper"]] <= 0.08,
  sensitivity_bound[["lower"]] >= 0.8,
  background_fpr_bound[["upper"]] <= 0.025,
  complete_null_upper <= 0.15,
  abs(mean(calibration["fitted_rho", ]) - true_rho) < 0.08
)

knitr::kable(calibration_summary, digits = 3)
```

The first three bounds use one-sided 95% t bounds for Monte Carlo means
across independent fixtures. For the complete-null family, any BH discovery
makes that replicate's FDP equal to one, so the table reports the exact
one-sided 95% Clopper--Pearson upper bound on its rejection probability. This
finite simulation is a regression stress test, not a proof of universal FDR
control. The package test repeats the same contract over 200 fixtures.

## Visualizing the Connectivity Map

A rank-ordered view shows both ground truth and the corrected decision. Red
circles are true discoveries, amber triangles are false discoveries, blue open
circles would mark missed network voxels, and gray open circles are correctly
rejected background voxels.

```{r plot-connectivity-map, fig.alt="Ranked seed-regressor t-statistics classified as true discoveries, false discoveries, missed network voxels, or correct background rejections."}
keep <- which(is.finite(t_seed))
ord <- keep[order(t_seed[keep], decreasing = TRUE)]
is_network <- ord %in% net_idx
is_discovered <- discovered[ord]
point_class <- ifelse(
  is_network & is_discovered, "true discovery",
  ifelse(
    !is_network & is_discovered, "false discovery",
    ifelse(is_network, "missed network voxel", "correct rejection")
  )
)
point_col <- c(
  "true discovery" = "firebrick",
  "false discovery" = "#E69F00",
  "missed network voxel" = "#0072B2",
  "correct rejection" = "gray70"
)
point_pch <- c(
  "true discovery" = 16,
  "false discovery" = 17,
  "missed network voxel" = 1,
  "correct rejection" = 1
)

stopifnot(length(ord) > 0)

plot(
  seq_along(ord),
  t_seed[ord],
  pch = unname(point_pch[point_class]),
  col = unname(point_col[point_class]),
  main = "Known-Truth Recovery after AR(1) GLS and BH-FDR",
  xlab = "Voxel rank",
  ylab = "AR(1) GLS seed-regressor t-statistic"
)
abline(h = 0, lty = 2, col = "gray40")
legend(
  "topright",
  legend = c(
    "True discovery",
    "False discovery",
    "Missed network voxel",
    "Correct rejection"
  ),
  col = unname(point_col),
  pch = unname(point_pch),
  cex = 0.8,
  bty = "n"
)
mtext(
  sprintf(
    "Single simulation: sensitivity %.2f | background FPR %.3f | FDP %.3f",
    sensitivity, background_fpr, false_discovery_proportion
  ),
  side = 3, line = 0.25, cex = 0.75
)
```

## Extending to Real Data

With real data, define the seed independently of the tested effect, then
extract a prespecified summary such as the mean time series across an ROI.
Avoid selecting the seed from the same connectivity statistic being tested;
that circular choice invalidates the inferential interpretation.

For a more complete analysis, prespecify nuisance regressors such as motion or
physiological signals; `baseline_model()` accepts them through
`nuisance_list`. If global signal regression or temporal filtering is used,
declare it explicitly because either choice changes the estimand and its
interpretation.

For whole-brain connectivity mapping, declare the temporal covariance model,
nuisance strategy, multiplicity procedure, and voxel inclusion mask before
interpreting discoveries. ROI-to-ROI analysis requires its own correction
across tested edges.

## Summary

This vignette used a seed covariate, a spline baseline, AR(1) GLS, and BH-FDR
correction to recover a known injected network. A negative-control coefficient
family and repeated simulations distinguish an appealing single map from
calibration evidence. The figure displays the corrected decision, while the
stress table reports sensitivity, null errors, and fitted AR behavior with
Monte Carlo bounds. Real-data validity still depends on defensible seed
definition, nuisance modeling, temporal covariance, and multiplicity control.