Fit an AR/ARMA noise model (run-aware) and return a whitening plan

Usage

fit_noise(
  resid = NULL,
  Y = NULL,
  X = NULL,
  runs = NULL,
  method = c("ar", "arma"),
  p = "auto",
  q = 0L,
  p_max = 6L,
  exact_first = c("ar1", "none"),
  pooling = c("global", "run", "parcel"),
  parcels = NULL,
  parcel_sets = NULL,
  multiscale = c("pacf_weighted", "acvf_pooled"),
  ms_mode = NULL,
  p_target = NULL,
  beta = 0.5,
  hr_iter = 0L,
  step1 = c("burg", "yw"),
  parallel = FALSE
)

Arguments

resid: Numeric matrix (time x voxels) of residuals from an initial OLS fit.
Y: Optional data matrix used to compute residuals when resid is omitted.
X: Optional design matrix used with Y to compute residuals.
runs: Optional integer vector of run identifiers.
method: Either "ar" or "arma".
p: AR order (integer or "auto" if method == "ar").
q: MA order (integer).
p_max: Maximum AR order when p = "auto".
exact_first: Apply exact AR(1) scaling at segment starts ("ar1" or "none").
pooling: Combine parameters across runs or parcels ("global", "run", "parcel").
parcels: Integer vector (length = ncol(resid)) giving fine parcel memberships when pooling = "parcel".
parcel_sets: Optional named list with entries coarse, medium, fine of equal length specifying nested parcel labels for multi-scale pooling.
multiscale: Multi-scale pooling mode when parcel_sets is supplied ("pacf_weighted" or "acvf_pooled"), or TRUE/FALSE to toggle pooling.
ms_mode: Explicit multiscale mode when multiscale is logical.
p_target: Target AR order for multi-scale pooling (defaults to p_max).
beta: Size exponent for multi-scale weights (default 0.5).
hr_iter: Number of Hannan–Rissanen refinement iterations for ARMA.
step1: Preliminary high-order AR fit method for HR ("burg" or "yw").
parallel: Reserved for future parallel estimation (logical).

Value

An object of class fmriAR_plan used by whiten_apply().

Examples

# Generate example data with AR(1) structure
n_time <- 200
n_voxels <- 50
phi_true <- 0.5

# Simulate residuals with AR(1) structure
resid <- matrix(0, n_time, n_voxels)
for (v in 1:n_voxels) {
  e <- rnorm(n_time)
  resid[1, v] <- e[1]
  for (t in 2:n_time) {
    resid[t, v] <- phi_true * resid[t-1, v] + e[t]
  }
}

# Fit AR model
plan <- fit_noise(resid, method = "ar", p = 1)

# With multiple runs
runs <- rep(1:2, each = 100)
plan_runs <- fit_noise(resid, runs = runs, method = "ar", pooling = "run")