Getting started with fmrilss • fmrilss

You have a preprocessed fMRI run and an event design with one regressor per trial. Your goal is a beta estimate for every trial and voxel, even though nearby haemodynamic responses overlap. Least squares separate (LSS) fits a model for each target trial, combining the other trials into one regressor. fmrilss returns the estimates as a trial-by-voxel matrix for multivariate pattern analysis, representational similarity analysis, connectivity, or reliability analysis.

This article follows one complete workflow: build a trial design, fit LSS, check the result, handle nuisance regressors correctly, and compare LSS with a simultaneous trial-wise GLM. The example is synthetic so that the truth is known; in a real analysis, Y and the design matrices come from your preprocessing and design pipeline.

What are the inputs?

lss() uses matrices with time along the rows:

Y contains the observed data: time points by voxels.
X contains one column per trial, convolved with the haemodynamic response function (HRF): time points by trials.
Z contains regressors retained in every trial-wise model, such as run intercepts or trends.
Nuisance contains regressors to project out, such as motion or physiological traces.

The function returns one row per trial and one column per voxel. Trial and voxel names are preserved when X and Y have column names.

suppressPackageStartupMessages({
  library(fmrilss)
  library(fmrihrf)
})

The following code constructs the trial design used throughout the article. Each event is an impulse convolved with the same canonical HRF, scaled to a unit peak so that trial betas are in units of peak response.

n_time <- 160
n_trials <- 24
n_voxels <- 8
sframe <- sampling_frame(blocklens = n_time, TR = 1)
onsets <- round(seq(8, 112, length.out = n_trials))
rset <- regressor_set(
  onsets,
  factor(seq_len(n_trials)),
  hrf = normalise_hrf(HRF_SPMG1),
  duration = 0,
  span = 30,
  summate = FALSE
)
X <- evaluate(
  rset,
  grid = samples(sframe, global = TRUE),
  precision = 0.1,
  method = "conv"
)
X <- as.matrix(X)
colnames(X) <- sprintf("trial_%02d", seq_len(n_trials))

We add an intercept, a trend, and two centred nuisance traces. The final lines create one synthetic Y with known trial effects. Replace this simulation with your own preprocessed time-by-voxel matrix.

set.seed(2401)
noise_sd <- 2.5
Z <- cbind(
  Intercept = 1,
  Linear_trend = as.numeric(scale(seq_len(n_time), scale = FALSE))
)
Nuisance <- scale(cbind(
  Motion_x = sin(seq_len(n_time) / 9),
  Motion_y = cos(seq_len(n_time) / 13)
), center = TRUE, scale = FALSE)

true_betas <- matrix(
  rnorm(n_trials * n_voxels),
  n_trials,
  n_voxels,
  dimnames = list(colnames(X), sprintf("voxel_%02d", seq_len(n_voxels)))
)
fixed_coef <- matrix(rnorm(ncol(Z) * n_voxels), ncol(Z), n_voxels)
nuisance_coef <- matrix(
  rnorm(ncol(Nuisance) * n_voxels, sd = 1.2),
  ncol(Nuisance),
  n_voxels
)
mean_signal <- X %*% true_betas +
  Z %*% fixed_coef +
  Nuisance %*% nuisance_coef
Y <- mean_signal + matrix(
  rnorm(n_time * n_voxels, sd = noise_sd),
  n_time,
  n_voxels,
  dimnames = list(NULL, colnames(true_betas))
)

Matrix shapes in the example workflow.
Object	Rows	Columns	Role
Y	160	8	observed voxel time series
X	160	24	one HRF-convolved column per trial
Z	160	2	shared modeled regressors
Nuisance	160	2	projected nuisance regressors
LSS result	24	8	one estimate per trial and voxel

How do I estimate the trial betas?

Start with the default optimized R backend and supply all parts of the model in the same call:

beta_lss <- lss(Y, X, Z = Z, Nuisance = Nuisance)
dim(beta_lss)
#> [1] 24  8

beta_lss has 24 trial rows and 8 voxel columns. Inspect the first few estimates:

round(beta_lss[1:4, 1:4], 2)
#>          voxel_01 voxel_02 voxel_03 voxel_04
#> trial_01    -3.41     1.94     0.43    -0.10
#> trial_02     1.56     2.75    -1.40    -0.77
#> trial_03     0.87     1.80    -2.84    -1.15
#> trial_04     2.39     0.02    -1.08    -0.59

Each row comes from a model with three parts: the target trial, the sum of all other trials, and the shared regressors. The target coefficient becomes the reported beta. Because overlapping trial regressors share variance, fitting them all at once makes their individual estimates compete for that variance. Repeating the three-part model for every trial reduces this competition among individual trial columns, but it does not make a poor or rank-deficient design informative.

Heatmap with 24 trial rows and 8 voxel columns on a colour scale centred at zero. Dark red cells are strongly positive and dark blue cells are strongly negative.

Estimated single-trial responses for the eight simulated voxels. Red is positive, blue is negative, and white is zero; deeper colour indicates larger magnitude, not statistical significance.

How should I handle nuisance regressors?

Pass Nuisance directly to lss(), as above. Before fitting the trial models, the package removes the shared effects of Z and Nuisance from both the response and trial design. This Frisch–Waugh–Lovell projection accounts for linearly dependent nuisance columns.

If you deliberately cache a projection for repeated analyses, apply it to every matrix that remains in the model. Here we project only Nuisance, so the same projection must transform Y, X, and Z:

Q_nuisance <- project_confounds(Nuisance)
beta_projected <- lss(
  Q_nuisance %*% Y,
  Q_nuisance %*% X,
  Z = Q_nuisance %*% Z
)
c(maximum_absolute_difference = max(abs(beta_lss - beta_projected)))
#> maximum_absolute_difference 
#>                1.993961e-13

The difference is at the level of numerical rounding: when the projection is applied to all three matrices, the two routes give the same estimates.

Leaving Z unprojected changes the model and can change the trial estimates. project_confounds() also materializes an $n \times n$ matrix; passing Nuisance to lss() avoids that storage and should be your default.

When does LSS differ from LSA?

Least Squares All (LSA) estimates every trial column simultaneously. LSS replaces the non-target columns with their sum in each trial-specific model. The choice involves a bias–variance trade-off.

To make the comparison reproducible, the next experiment holds X, true_betas, and all nuisance effects fixed. Because both estimators are linear here, we can calculate conditional squared bias and noise variance directly. We also generate 100 independent noise realizations as a check on that calculation.

The first table below checks that the trial design remains estimable after the shared Z and Nuisance columns are removed (the “residualized” design): its rank equals the number of trials. The second table gives the exact decomposition of mean squared error (MSE) at the noise level used to simulate Y.

Diagnostics for the residualized trial design used in the comparison.
Diagnostic	Value
Residualized trial-design rank	24.000
Largest absolute trial correlation	0.378
Scaled condition number	8.706

Exact conditional decomposition at noise SD = 2.5.
Estimator	Squared bias (beta^2)	Variance (beta^2)	MSE (beta^2)	RMSE (beta)
LSS	0.374	1.927	2.301	1.517
LSA	0.000	5.423	5.423	2.329

The 100-repetition check gave MSE = 2.285 for LSS and 5.296 for LSA, close to the exact values in the table.

Grouped bar chart comparing LSS and LSA squared bias, variance, and mean squared error on a common beta-squared scale. LSS has higher squared bias and lower variance and mean squared error.

At noise SD 2.5, LSS has more squared bias but sufficiently lower variance to produce lower mean squared error than LSA.

Which estimator wins depends on the noise level. As noise grows, the variance term grows while the squared bias stays fixed, so lower variance matters more. The crossover is the noise SD at which LSS and LSA have equal MSE; for this fixed design and truth it is SD = 0.818:

Conditional RMSE on both sides of the design-specific crossover.
Noise SD	LSS RMSE	LSA RMSE
0.410	0.653	0.382
0.818	0.762	0.762
2.500	1.517	2.329

At noise SD 0.41, below the crossover, LSA has lower RMSE; at the illustrated SD 2.5, LSS has lower RMSE. The preferred estimator therefore depends on the noise level even in this one design. Event spacing, noise, HRF mismatch, and the distribution of true trial effects can all move the crossover. Simulate conditions that resemble your own acquisition.

What should I change for real data?

Real fMRI errors are temporally correlated. Once the design is correct, pass a prewhitening specification so the response and every design matrix receive the same filter. This call also selects the optimized C++ backend; prewhitening has the same model meaning across backends:

beta_ar1 <- lss(
  Y,
  X,
  Z = Z,
  Nuisance = Nuisance,
  method = "cpp_optimized",
  prewhiten = prewhiten_options(method = "ar", p = 1)
)

The simulated errors above are independent, so this call demonstrates the API and backend agreement, not a benefit from whitening. Treat AR(1) as an example starting model: choose an order and pooling strategy using residual diagnostics and acquisition-aware validation.

How the noise model is estimated matters for single-trial designs: see vignette("prewhitening") for the residual_model choices and when each is appropriate.

For multiple runs, use run-aware intercepts and supply run labels to prewhiten_options(method = "ar", p = 1, pooling = "run", runs = ...). Do not concatenate runs and let filtering cross a run boundary.

For a large voxel matrix, the optimized C++ backend changes computation, not the estimand:

beta_cpp <- lss(
  Y,
  X,
  Z = Z,
  Nuisance = Nuisance,
  method = "cpp_optimized"
)
c(maximum_absolute_difference = max(abs(beta_lss - beta_cpp)))
#> maximum_absolute_difference 
#>                 4.72955e-14

Start with the default backend while developing the analysis, then switch when runtime or memory justifies it. Backend agreement is a useful regression check, but agreement alone does not validate the scientific model.

Three options refine the trial-wise model itself. They are supported by the r_optimized, cpp_optimized, and cpp backends (and the naive reference for trial_groups and ridge), and can be combined in one call:

trial_groups gives one summed “other trials” regressor per condition (LSS-N; Turner et al., 2012) instead of pooling every other trial. Use it when conditions are expected to evoke different responses.
ridge adds a fractional ridge penalty to each trial model, shrinking trial estimates toward zero. Give one value, or c(trial, others) to penalize the target-trial and other-trials coefficients separately; each value is a fraction of the mean design energy. In rapid designs, where neighbouring trials overlap, this can trade a little bias for a large reduction in variance; choose the amount by validation on your design.
prewhiten_options(pooling = "voxel") fits a voxel-adaptive noise model: voxels are binned by residual autocorrelation and each bin is fitted with its own filtered design.

conditions <- rep(c("A", "B"), length.out = ncol(X))
beta_lss_n <- lss(Y, X, Z = Z, Nuisance = Nuisance, trial_groups = conditions)
beta_ridge <- lss(
  Y, X, Z = Z, Nuisance = Nuisance,
  trial_groups = conditions,
  ridge = c(0.5, 0),
  prewhiten = prewhiten_options(method = "ar", p = 1, pooling = "voxel",
                                voxel_bins = 4)
)
dim(beta_ridge)
#> [1] 24  8

What are the important limits?

X must already encode the HRF-convolved response for each trial. If you begin from event tables, build X with a design-aware workflow such as vignette("lss_with_fmridesign").
LSS reduces competition among individual trial columns; it cannot recover information absent from the acquisition or rescue a singular design.
The output contains trial coefficients, not significance tests or an automatic decision about whether a trial is usable.
Nuisance projection and prewhitening change the model. Apply them consistently to the response and every relevant design matrix.
LSS versus LSA is a design- and noise-dependent choice. Validate it with repeated simulations, not a single attractive plot.

Where should I go next?

Read ?lss and ?prewhiten_options for the complete core interface.
Continue to vignette("oasis_method") for the OASIS solver, explicit ridge regularization, diagnostics, and multi-basis designs.
Then read vignette("oasis_theory") for the derivation and computational calculations behind that practical guide.
Next, use vignette("lss_with_fmridesign") when your inputs are event tables, formulas, and multiple runs.
Continue to vignette("voxel-wise-hrf") for estimated voxel-wise HRF shapes and trial coefficients.
Use vignette("sbhm") for library-constrained voxel-specific shapes, score margins, and trial coefficients.
Use vignette("glmsingle") for library HRF selection, GLMdenoise, and voxel-wise ridge selection when conditions repeat across runs.

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. https://doi.org/10.1016/j.neuroimage.2011.08.076

Turner, B. O., Mumford, J. A., Poldrack, R. A., & Ashby, F. G. (2012). Spatiotemporal activity estimation for multivoxel pattern analysis with rapid event-related designs. NeuroImage, 62(3), 1429–1438. https://doi.org/10.1016/j.neuroimage.2012.05.057