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

Usage

fit_noise(
  resid = NULL,
  Y = NULL,
  X = NULL,
  runs = NULL,
  censor = NULL,
  method = c("ar", "arma"),
  p = "auto",
  q = 0L,
  p_max = 6L,
  q_max = 2L,
  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"),
  design = NULL,
  acvf_correction = NULL,
  correction_max_lag = 25L,
  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 run labels, one per timepoint. Each label must occupy one contiguous block and may not be missing.

censor

Optional integer vector of 1-based timepoint indices to exclude from AR parameter estimation, or a logical vector of length nrow(resid) where TRUE

indicates censored timepoints. Censored frames (e.g., motion-corrupted) are excluded when computing autocorrelations. Each run's estimation uses only its own valid (non-censored) segments.

method

Either "ar" or "arma". ARMA models are estimated by a Hannan–Rissanen regression pooled over voxels (every voxel is treated as a replicate series sharing the coefficients), with all lags kept inside contiguous run/censor segments.

p

AR order (integer or "auto" if method == "ar"). An explicit order is honoured even when it exceeds p_max.

q

MA order (integer), or "auto" (ARMA only) to choose it by BIC from 0:q_max.

p_max

Maximum AR order when p = "auto". The order is chosen by BIC whose sample size is the number of frames times the effective number of independent voxels (a design effect from the mean inter-voxel correlation), so pooling many independent voxels can support a higher order than one series would. For method = "arma", p = "auto" searches 0:p_max jointly with the MA order by BIC on the pooled Hannan–Rissanen regression.

q_max

Maximum MA order when q = "auto".

exact_first

How whiten_apply() treats the start of each run or post-censoring segment. "ar1" (default; the name is historical) applies the exact stationary initialisation of the fitted AR/ARMA model, i.e. exact GLS whitening within each segment; for AR(1) this is the familiar sqrt(1 - phi^2) scaling of the first sample. "none" uses the conditional filter with a zero pre-sample.

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").

design

Optional design matrix (timepoints x regressors) whose projection produced resid. Supplying it corrects the downward bias that projecting a design out of the data puts into the autocovariance, and hence into phi. Opt-in, because it changes estimates and needs the design to be the one that actually formed the OLS residuals. A numerical orthogonality check rejects detectable mismatches, including typical raw, GLS, robust, or unrelated residuals. Orthogonality cannot prove the exact residual-forming provenance, so the caller remains responsible for supplying the matching design. Supported for every pooling mode with method = "ar"; under pooling = "parcel" the per-run bias maps are combined in proportion to the lag pairs each run contributes.

acvf_correction

Precomputed bias matrices from acvf_bias_matrix(), as an alternative to design when many datasets share one design. A single matrix is applied to every run; a list is matched against the runs in order. Mutually exclusive with design.

correction_max_lag

Lag budget for the bias correction (default 25). The correction solves a system truncated at this lag, so too small a budget leaves bias behind; too large a one approaches the run length and makes the system ill-conditioned, which is refused with a warning rather than solved. A design leaving fewer residual degrees of freedom than the budget also cannot support it, and the budget is reduced accordingly, again with a warning. Ignored unless design is supplied.

parallel

Reserved for future parallel estimation (logical).

Value

An object of class fmriAR_plan used by whiten_apply(). Besides the AR/MA coefficients the plan carries the noise scale and shape it was fitted from, so consumers can reconstruct the covariance it implies rather than only its correlation structure:

gamma: list of autocovariance vectors, one per pooling unit – a single entry for pooling = "global", one per run for pooling = "run". Lags run 0 to the highest the data supported, which is governed by p_max and the run length rather than by p, so fit_noise(p = 1, p_max = 6) returns seven values, not two. Under global pooling every run is truncated to the shortest available length before averaging, since a zero-padded autocovariance is not a valid covariance. Global AR coefficients are fitted by Yule–Walker on the frame-weighted average of the runs' autocorrelations.
sigma2: list of innovation variances, matching gamma. For method = "ar" it is derived as gamma_0 - sum_k phi_k gamma_k from the coefficients stored on the plan so the two are always mutually consistent, and is NA whenever gamma does not reach lag length(phi) – heavy censoring can truncate it that far, and a partial sum would overstate the innovation variance rather than report that it is unavailable. For method = "arma" it is the mean squared exact-whitened residual over all voxels and valid frames.
gamma_by_parcel, sigma2_by_parcel: the same quantities per parcel when pooling = "parcel", keyed like phi_by_parcel.

For a run-stationary noise process with autocovariance gamma, the covariance of the data within a run is the Toeplitz matrix built from it, which is what makes design-specific variance calculations possible downstream without refitting.

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")