Least Squares Separate (LSS) Analysis
lss.RdComputes trial-wise beta estimates using the Least Squares Separate approach of Mumford et al. (2012). This method fits a separate GLM for each trial, with the trial of interest and all other trials as separate regressors.
Arguments
- Y
A numeric matrix of size n × V where n is the number of timepoints and V is the number of voxels/variables
- X
A numeric matrix of size n × T where T is the number of trials. Each column represents the design for one trial
- Z
A numeric matrix of size n × F representing experimental regressors to include in all trial-wise models. These are regressors we want to model and get beta estimates for, but not trial-wise (e.g., intercept, condition effects, block effects). If NULL, an intercept-only design is used. Defaults to NULL
- Nuisance
A numeric matrix of size n × N representing nuisance regressors to be projected out before LSS analysis (e.g., motion parameters, physiological noise). If NULL, no nuisance projection is performed. Defaults to NULL
- method
Character string specifying which implementation to use. Options are:
"r_optimized" - Optimized R implementation (recommended, default)
"cpp_optimized" - Optimized C++ implementation with parallel support
"r_vectorized" - Standard R vectorized implementation
"cpp" - Standard C++ implementation
"naive" - Simple loop-based R implementation (for testing)
"oasis" - OASIS method with HRF support and ridge regularization
- block_size
An integer specifying the voxel block size for parallel processing, only applicable when `method = "cpp_optimized"`. Defaults to 96.
- oasis
A list of options for the OASIS method. See Details for available options.
- prewhiten
A list of prewhitening options using fmriAR. See Details for available options.
Value
A numeric matrix of size T × V containing the trial-wise beta estimates. Note: Currently only returns estimates for the trial regressors (X). Beta estimates for the experimental regressors (Z) are computed but not returned.
Details
The LSS approach fits a separate GLM for each trial, where each model includes:
The trial of interest (from column i of X)
All other trials combined (sum of all other columns of X)
Experimental regressors (Z matrix) - these are modeled to get beta estimates but not trial-wise
If Nuisance regressors are provided, they are first projected out from both Y and X using standard linear regression residualization.
When using method="oasis", the following options are available in the oasis list:
design_spec: A list for building trial-wise designs from event onsets using fmrihrf. Must contain:sframe(sampling frame),cond(list withonsets,hrf, and optionallyspan), and optionallyothers(list of other conditions to be modeled as nuisances). When provided, X can be NULL and will be constructed automatically.K: Explicit basis dimension for multi-basis HRF models (e.g., 3 for SPMG3). If not provided, it's auto-detected from X dimensions or defaults to 1 for single-basis HRFs.ridge_mode: Either "absolute" (default) or "fractional". In absolute mode, ridge_x and ridge_b are used directly as regularization parameters. In fractional mode, they represent fractions of the maximum eigenvalue for adaptive regularization.ridge_x: Ridge parameter for trial-specific regressors (default 0). Controls regularization strength for individual trial estimates.ridge_b: Ridge parameter for the aggregator regressor (default 0). Controls regularization strength for the sum of all other trials.return_se: Logical, whether to return standard errors (default FALSE). When TRUE, returns a list withbeta(trial estimates) andse(standard errors) components.return_diag: Logical, whether to return design diagnostics (default FALSE). When TRUE, includes diagnostic information about the design matrix structure.block_cols: Integer, voxel block size for memory-efficient processing (default 4096). Larger values use more memory but may be faster for systems with sufficient RAM.whiten: Logical, whether to apply AR(1) whitening (default FALSE). When TRUE, estimates AR(1) coefficients and pre-whitens data to account for temporal autocorrelation.ntrials: Explicit number of trials (used when K > 1 to determine output dimensions). If not provided, calculated as ncol(X) / K.hrf_grid: Vector of HRF indices for grid-based HRF selection (advanced use). Allows testing multiple HRF shapes simultaneously.
When using the prewhiten parameter, the following options are available:
method: Character, "ar" (default), "arma", or "none" for the noise model type.p: Integer or "auto" for AR order (default "auto").q: Integer for MA order in ARMA models (default 0).p_max: Maximum AR order when p="auto" (default 6).pooling: Character, "global" (default), "voxel", "run", or "parcel" for parameter estimation strategy.runs: Integer vector of run identifiers for run-aware estimation.parcels: Integer vector of parcel memberships for parcel-based pooling.exact_first: Character, "ar1" or "none" for exact AR(1) scaling at segment starts.
Prewhitening is applied before the LSS analysis to account for temporal autocorrelation in the
fMRI time series. The fmriAR package provides flexible AR/ARMA modeling with various pooling
strategies. For backward compatibility, the old oasis$whiten = "ar1" syntax is still
supported and will be converted to the equivalent prewhiten settings.
The OASIS method provides a mathematically equivalent but computationally optimized version of standard LSS. It reformulates the per-trial GLM fitting as a single matrix operation, eliminating redundant computations. This is particularly beneficial for designs with many trials or when processing large datasets. When K > 1 (multi-basis HRFs), the output will have K*ntrials rows, with basis functions for each trial arranged sequentially.
References
Mumford, J. A., Turner, B. O., Ashby, F. G., & Poldrack, R. A. (2012). Deconvolving BOLD activation in event-related designs for multivoxel pattern classification analyses. NeuroImage, 59(3), 2636-2643.
Examples
# Generate example data
n_timepoints <- 100
n_trials <- 10
n_voxels <- 50
# Create trial design matrix
X <- matrix(0, n_timepoints, n_trials)
for(i in 1:n_trials) {
start <- (i-1) * 8 + 1
if(start + 5 <= n_timepoints) {
X[start:(start+5), i] <- 1
}
}
# Create data with some signal
Y <- matrix(rnorm(n_timepoints * n_voxels), n_timepoints, n_voxels)
true_betas <- matrix(rnorm(n_trials * n_voxels, 0, 0.5), n_trials, n_voxels)
for(i in 1:n_trials) {
Y <- Y + X[, i] %*% matrix(true_betas[i, ], 1, n_voxels)
}
# Run LSS analysis
beta_estimates <- lss(Y, X)
# With experimental regressors (intercept + condition effects)
Z <- cbind(1, scale(1:n_timepoints))
beta_estimates_with_regressors <- lss(Y, X, Z = Z)
# With nuisance regression (motion parameters)
Nuisance <- matrix(rnorm(n_timepoints * 6), n_timepoints, 6)
beta_estimates_clean <- lss(Y, X, Z = Z, Nuisance = Nuisance)
if (FALSE) { # \dontrun{
# Using OASIS method with ridge regularization
beta_oasis <- lss(Y, X, method = "oasis",
oasis = list(ridge_x = 0.1, ridge_b = 0.1,
ridge_mode = "fractional"))
# OASIS with standard errors
result_with_se <- lss(Y, X, method = "oasis",
oasis = list(return_se = TRUE))
beta_estimates <- result_with_se$beta
standard_errors <- result_with_se$se
# Building design from event onsets using fmrihrf
sframe <- sampling_frame(blocklens = 200, TR = 1.0)
# OASIS with automatic design construction
beta_auto <- lss(Y, X = NULL, method = "oasis",
oasis = list(
design_spec = list(
sframe = sframe,
cond = list(
onsets = c(10, 30, 50, 70, 90, 110, 130, 150),
hrf = HRF_SPMG1,
span = 25
),
others = list(
list(onsets = c(20, 40, 60, 80, 100, 120, 140))
)
)
))
# Multi-basis HRF example (3 basis functions per trial)
beta_multibasis <- lss(Y, X = NULL, method = "oasis",
oasis = list(
design_spec = list(
sframe = sframe,
cond = list(
onsets = c(10, 30, 50, 70, 90),
hrf = HRF_SPMG3, # 3-basis HRF
span = 30
)
),
K = 3 # Explicit basis dimension
))
# Returns 15 rows (5 trials * 3 basis functions)
} # }