| Title: | Regression Analysis of Functional Magnetic Resonance Imaging Data |
|---|---|
| Description: | Provides tools for the analysis of functional magnetic resonance imaging (fMRI) data. Facilities are included to construct flexible hemodynamic response functions, experimental regressors, and conduct univariate fMRI regression analyses. |
| Authors: | Bradley Buchsbaum [aut, cre] |
| Maintainer: | Bradley Buchsbaum <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 0.3.0 |
| Built: | 2026-10-11 23:59:21 UTC |
| Source: | https://github.com/bbuchsbaum/fmrireg |
Plugin registration and helpers for fmrireg
.fmrireg_engine_registry.fmrireg_engine_registry
An object of class environment of length 2.
Helper for drawing per-event amplitudes/durations around base values.
.resample_param( base, sd, dist = c("lognormal", "gamma", "gaussian"), allow_negative = FALSE ).resample_param( base, sd, dist = c("lognormal", "gamma", "gaussian"), allow_negative = FALSE )
base |
Numeric vector of base values to jitter |
sd |
Numeric standard deviation for the sampling distribution |
dist |
Distribution name: "lognormal", "gamma", or "gaussian" |
allow_negative |
Logical; if TRUE, allow negative draws (only used for gaussian) |
A numeric vector of resampled values with the same length as base
Applies the soft projection to both the data matrix Y and design matrix X. Both must be projected to avoid bias in coefficient estimates.
apply_soft_projection(proj, Y, X)apply_soft_projection(proj, Y, X)
proj |
A soft_projection object from |
Y |
Data matrix (time x voxels). |
X |
Design matrix (time x predictors). |
List with projected Y and X.
set.seed(123) N <- matrix(rnorm(100 * 20), nrow = 100, ncol = 20) Y <- matrix(rnorm(100 * 50), nrow = 100, ncol = 50) X <- cbind(1, rnorm(100)) proj <- soft_projection(N, lambda = "auto") cleaned <- apply_soft_projection(proj, Y, X) # Use cleaned$Y and cleaned$X for GLM fittingset.seed(123) N <- matrix(rnorm(100 * 20), nrow = 100, ncol = 20) Y <- matrix(rnorm(100 * 50), nrow = 100, ncol = 50) X <- cbind(1, rnorm(100)) proj <- soft_projection(N, lambda = "auto") cleaned <- apply_soft_projection(proj, Y, X) # Use cleaned$Y and cleaned$X for GLM fitting
Retrieves the estimated autoregressive parameters from a fitted fMRI linear model that used AR error modeling.
ar_parameters(object, ...) ## S3 method for class 'fmri_lm' ar_parameters(object, scope = c("average", "per_run", "raw"), ...)ar_parameters(object, ...) ## S3 method for class 'fmri_lm' ar_parameters(object, scope = c("average", "per_run", "raw"), ...)
object |
An object of class |
... |
Additional arguments (currently unused) |
scope |
Character; |
Depending on scope, either a numeric vector of averaged AR
coefficients, a list of per-run coefficient vectors, or the raw stored
structure. Returns NULL when no AR modeling was performed.
## Not run: # Fit model with AR(1) errors fit <- fmri_lm(onset ~ hrf(cond), dataset = dset, cor_struct = "ar1") ar_parameters(fit) # Extract estimated AR(1) coefficient ## End(Not run)## Not run: # Fit model with AR(1) errors fit <- fmri_lm(onset ~ hrf(cond), dataset = dset, cor_struct = "ar1") ar_parameters(fit) # Extract estimated AR(1) coefficient ## End(Not run)
Wraps a list of bindings (or a manifest data.frame) into an
fmri_manifest for instantiate(). This is the generic, BIDS-free way
to describe subjects; see from_bids() for the BIDS populator.
as_manifest(x)as_manifest(x)
x |
A list of bindings (named lists, each with at least |
An object of class fmri_manifest.
m <- as_manifest(list( list(id = "sub-01", scans = matrix(rnorm(80 * 2), 80, 2), TR = 2, run_length = c(40, 40), events = data.frame(onset = c(5, 45), condition = factor(c("A", "B")), run = c(1, 2))))) length(m)m <- as_manifest(list( list(id = "sub-01", scans = matrix(rnorm(80 * 2), 80, 2), TR = 2, run_length = c(40, 40), events = data.frame(onset = c(5, 45), condition = factor(c("A", "B")), run = c(1, 2))))) length(m)
This method allows NeuroVec objects from neuroim2 to be converted to base R arrays using the standard as.array() function. This is particularly useful in testing and data manipulation contexts.
## S3 method for class 'NeuroVec' as.array(x, ...)## S3 method for class 'NeuroVec' as.array(x, ...)
x |
A NeuroVec object from neuroim2 |
... |
Additional arguments (currently unused) |
A base R array containing the data from the NeuroVec object
Visual summaries of an fmri_lm fit. Five views are available:
## S3 method for class 'fmri_lm' autoplot( object, type = c("betas", "contrasts", "hrf", "timecourse", "residuals"), voxel = NULL, level = 0.95, sample_at = seq(0, 20, by = 0.25), direct_labels = TRUE, ... ) ## S3 method for class 'fmri_lm' plot( x, type = c("betas", "contrasts", "hrf", "timecourse", "residuals"), voxel = NULL, ... )## S3 method for class 'fmri_lm' autoplot( object, type = c("betas", "contrasts", "hrf", "timecourse", "residuals"), voxel = NULL, level = 0.95, sample_at = seq(0, 20, by = 0.25), direct_labels = TRUE, ... ) ## S3 method for class 'fmri_lm' plot( x, type = c("betas", "contrasts", "hrf", "timecourse", "residuals"), voxel = NULL, ... )
object, x
|
An |
type |
One of |
voxel |
Optional integer voxel index (or indices, for
|
level |
Confidence level for intervals and bands. |
sample_at |
Time points (s) at which to evaluate the HRF, for
|
direct_labels |
Logical; for |
... |
Passed on from |
"betas"With no voxel, the distribution of
t-statistics across voxels for each event regressor, drawn as a sina
plot (points spread in proportion to their density), with reference
lines at the two-sided p < .001 (uncorrected) critical value and the
share of voxels beyond it in each tail. With voxel, a coefficient
plot of estimates with confidence intervals, one panel per event term
so that regressors in different units (a condition and a parametric
modulator) do not share an axis.
"contrasts"As "betas", for the fitted contrasts.
"hrf"The estimated hemodynamic response for each condition: the term's HRF basis weighted by the fitted coefficients. For one voxel, the band is a pointwise confidence interval from the coefficient covariance, and a curve whose interval includes zero at its peak is drawn faded and dashed and labelled "n.s.". For several voxels, the line is the mean across voxels, the band a t-based confidence interval of that mean (voxels treated as independent), and up to 200 voxels are drawn faintly underneath. Peak latency is reported only for terms with more than one basis function; with a single basis the shape, and so the latency, is fixed by the HRF rather than estimated.
"timecourse"Observed signal and model fit for one voxel, with event onsets in a lane per condition and residuals in a lower panel drawn on the same vertical scale.
"residuals"Autocorrelation of the residuals at lags 1 to
10 scans. By default this summarises up to 200 voxels spread evenly
through the data (median and 10th-90th percentile per lag); with
voxel, it shows those voxels. Least-squares residuals are
negatively autocorrelated even when the noise is white, so the 95\
white-noise range is centred on the autocorrelation this design
induces, not on zero. For fits with an AR noise model the residuals
are also shown after prewhitening with the fitted AR coefficients,
which is the check that the noise model removed the serial
correlation. The subtitle reports a portmanteau test of lags 1-10 that
accounts for both the design-induced mean and the correlation between
lags that regression introduces (a generalised Ljung-Box test).
When voxel is NULL for the "hrf" and
"timecourse" views, the voxel with the largest absolute event
t-statistic is shown and named in the subtitle.
Critical values and intervals use the degrees of freedom the fit used for
inference (fit$result$df$inference), which differ from the nominal
residual df for AR, robust, or effective-df corrected fits.
In the "timecourse" and "residuals" views the fitted values
are an ordinary least squares fit of the full design (events plus baseline)
to the voxel: run by run for runwise fits, jointly otherwise. For OLS fits
these are the model's own fitted values; for AR or robust fits they are the
unweighted counterpart, on the scale of the data. For runwise fits the
fitted values use run-specific coefficients, whereas the estimates view
shows coefficients pooled across runs. The reported partial R-squared is
: the share of the variance left by
the baseline model that the event regressors explain.
A ggplot2 object. plot() prints it and returns it invisibly.
set.seed(1) ev <- data.frame(onset = seq(10, 180, by = 12), run = 1, condition = factor(rep(c("a", "b"), length.out = 15))) Y <- matrix(rnorm(100 * 5), 100, 5) dset <- matrix_frame(Y, TR = 2, run_length = 100, event_table = ev) fit <- fmri_lm(onset ~ hrf(condition), block = ~ run, dataset = dset) ggplot2::autoplot(fit) ggplot2::autoplot(fit, type = "timecourse", voxel = 1)set.seed(1) ev <- data.frame(onset = seq(10, 180, by = 12), run = 1, condition = factor(rep(c("a", "b"), length.out = 15))) Y <- matrix(rnorm(100 * 5), 100, 5) dset <- matrix_frame(Y, TR = 2, run_length = 100, event_table = ev) fit <- fmri_lm(onset ~ hrf(condition), block = ~ run, dataset = dset) ggplot2::autoplot(fit) ggplot2::autoplot(fit, type = "timecourse", voxel = 1)
Plots the predicted BOLD time course of an fMRI regressor: the event train convolved with its HRF. Event onsets are marked as ticks beneath the curve, so the hemodynamic lag is visible. When responses to nearby events overlap, each event's own response is drawn faintly under the summed curve. Multi-basis HRFs are drawn as stacked panels that share the time axis, each with its own vertical scale.
## S3 method for class 'Reg' autoplot( object, grid = NULL, precision = 0.1, method = "conv", title = NULL, ... )## S3 method for class 'Reg' autoplot( object, grid = NULL, precision = 0.1, method = "conv", title = NULL, ... )
object |
A |
grid |
Optional numeric vector specifying time points (seconds) for evaluation. If NULL, a default grid spans the first onset to the end of the last event's response (including its duration). |
precision |
Numeric precision for HRF evaluation if |
method |
Evaluation method passed to |
title |
Optional plot title, e.g. the condition this regressor models. |
... |
Additional arguments (currently unused). |
A ggplot object.
Describes how the nuisance / baseline model should be built for every
subject, without binding any particular subject's confound values. The
concrete baseline_model is assembled per subject at instantiation
time from this spec, the subject's sampling_frame, and the subject's
confound matrix.
baseline_spec( degree = 3, basis = c("bs", "poly", "ns", "constant"), confounds = NULL, intercept = c("runwise", "global", "none"), nuisance_check = c("warn", "error", "drop", "none") )baseline_spec( degree = 3, basis = c("bs", "poly", "ns", "constant"), confounds = NULL, intercept = c("runwise", "global", "none"), nuisance_check = c("warn", "error", "drop", "none") )
degree |
Integer drift degree passed to |
basis |
Drift basis: one of |
confounds |
Optional confound selection used to populate the per-subject
nuisance regressors. Either |
intercept |
Intercept handling passed to |
nuisance_check |
How to handle problematic nuisance regressors, passed
to |
An object of class baseline_spec.
fmri_template(), baseline_model()
baseline_spec(degree = 3, basis = "bs")baseline_spec(degree = 3, basis = "bs")
Extract per-job errors from a batch result
batch_errors(x)batch_errors(x)
x |
An |
A named character vector (by job id) of error messages for jobs that failed; empty if all succeeded.
Extract successful values from a batch result
batch_values(x)batch_values(x)
x |
An |
A named list (by job id) of reducer outputs for jobs that succeeded.
Builds the per-subject baseline_model (from the template's
baseline_spec, the dataset's sampling frame, and any nuisance regressors)
and combines it with the event model into an fmri_model.
build_model(job, dataset = realize_dataset(job))build_model(job, dataset = realize_dataset(job))
job |
An fmri_job. |
dataset |
The realized dataset (defaults to |
An fmri_model.
Extract Image/Volume for Coefficient
## S3 method for class 'fmri_lm' coef_image( object, coef = 1, statistic = c("estimate", "se", "tstat", "prob"), type = c("estimates", "contrasts", "F"), ... ) coef_image(object, coef = 1, statistic = c("estimate", "se", "z", "p"), ...)## S3 method for class 'fmri_lm' coef_image( object, coef = 1, statistic = c("estimate", "se", "tstat", "prob"), type = c("estimates", "contrasts", "F"), ... ) coef_image(object, coef = 1, statistic = c("estimate", "se", "z", "p"), ...)
object |
An fmri_meta object |
coef |
Coefficient name or index |
statistic |
Type of statistic to extract ("estimate", "se", "z", "p") |
type |
For |
... |
Additional arguments (currently unused). |
NeuroVol object or matrix
# Create a small example X <- matrix(rnorm(50 * 4), 50, 4) edata <- data.frame( condition = factor(c("A", "B", "A", "B")), onsets = c(1, 12, 25, 38), run = c(1, 1, 1, 1) ) dset <- matrix_frame(X, TR = 2, run_length = 50, event_table = edata) fit <- fmri_lm(onsets ~ hrf(condition), block = ~run, dataset = dset) # Get coefficient estimates as a numeric vector coef_image(fit, coef = 1) toy_meta <- structure( list( coefficients = matrix(c(0.3, 0.1), nrow = 1, dimnames = list(NULL, c("A", "B"))), se = matrix(c(0.05, 0.06), nrow = 1) ), class = "fmri_meta" ) coef_image(toy_meta, coef = "A")# Create a small example X <- matrix(rnorm(50 * 4), 50, 4) edata <- data.frame( condition = factor(c("A", "B", "A", "B")), onsets = c(1, 12, 25, 38), run = c(1, 1, 1, 1) ) dset <- matrix_frame(X, TR = 2, run_length = 50, event_table = edata) fit <- fmri_lm(onsets ~ hrf(condition), block = ~run, dataset = dset) # Get coefficient estimates as a numeric vector coef_image(fit, coef = 1) toy_meta <- structure( list( coefficients = matrix(c(0.3, 0.1), nrow = 1, dimnames = list(NULL, c("A", "B"))), se = matrix(c(0.05, 0.06), nrow = 1) ), class = "fmri_meta" ) coef_image(toy_meta, coef = "A")
Creates a NeuroVol image from coefficients in an fmri_ttest_fit object.
## S3 method for class 'fmri_ttest_fit' coef_image(object, coef = 1, statistic = c("estimate", "se", "z", "p"), ...)## S3 method for class 'fmri_ttest_fit' coef_image(object, coef = 1, statistic = c("estimate", "se", "z", "p"), ...)
object |
An fmri_ttest_fit object |
coef |
Character or integer; coefficient to extract |
statistic |
Character string; type of statistic to extract:
|
... |
Additional arguments (e.g., mask to apply) |
NeuroVol object or numeric vector
A plural companion to coef_image. Returns a named list of
volumes (or numeric vectors for non-spatial datasets), one per coefficient,
so callers can iterate over or write out every coefficient without
reimplementing the loop over coef_names.
## S3 method for class 'fmri_lm' coef_images( object, statistic = c("estimate", "se", "tstat", "prob"), type = c("estimates", "contrasts", "F"), coefs = NULL, ... ) coef_images(object, ...)## S3 method for class 'fmri_lm' coef_images( object, statistic = c("estimate", "se", "tstat", "prob"), type = c("estimates", "contrasts", "F"), coefs = NULL, ... ) coef_images(object, ...)
object |
A fitted model object. |
statistic |
For |
type |
For |
coefs |
Optional character vector of coefficient names (a subset of
|
... |
Additional arguments passed to methods (and on to
|
A named list of NeuroVol objects (or numeric vectors when the
dataset is non-spatial), named by coefficient.
# Create a small example X <- matrix(rnorm(50 * 4), 50, 4) edata <- data.frame( condition = factor(c("A", "B", "A", "B")), onsets = c(1, 12, 25, 38), run = c(1, 1, 1, 1) ) dset <- matrix_frame(X, TR = 2, run_length = 50, event_table = edata) fit <- fmri_lm(onsets ~ hrf(condition), block = ~run, dataset = dset) # Named list of one volume per event regressor imgs <- coef_images(fit, statistic = "estimate", type = "estimates") names(imgs)# Create a small example X <- matrix(rnorm(50 * 4), 50, 4) edata <- data.frame( condition = factor(c("A", "B", "A", "B")), onsets = c(1, 12, 25, 38), run = c(1, 1, 1, 1) ) dset <- matrix_frame(X, TR = 2, run_length = 50, event_table = edata) fit <- fmri_lm(onsets ~ hrf(condition), block = ~run, dataset = dset) # Named list of one volume per event regressor imgs <- coef_images(fit, statistic = "estimate", type = "estimates") names(imgs)
Return the names of available coefficients from a fitted model object.
This helps users discover which coefficient names can be passed to
coef_image or other extraction functions.
coef_names(x, ...) ## S3 method for class 'fmri_lm' coef_names(x, type = c("estimates", "contrasts", "F", "all"), ...)coef_names(x, ...) ## S3 method for class 'fmri_lm' coef_names(x, type = c("estimates", "contrasts", "F", "all"), ...)
x |
A fitted model object |
... |
Additional arguments passed to methods |
type |
Which set of names to return: |
A character vector of coefficient names
Other statistical_measures:
p_values(),
standard_error(),
stats()
# Create a small example X <- matrix(rnorm(50 * 4), 50, 4) edata <- data.frame( condition = factor(c("A", "B", "A", "B")), onsets = c(1, 12, 25, 38), run = c(1, 1, 1, 1) ) dset <- matrix_frame(X, TR = 2, run_length = 50, event_table = edata) fit <- fmri_lm(onsets ~ hrf(condition), block = ~run, dataset = dset) coef_names(fit)# Create a small example X <- matrix(rnorm(50 * 4), 50, 4) edata <- data.frame( condition = factor(c("A", "B", "A", "B")), onsets = c(1, 12, 25, 38), run = c(1, 1, 1, 1) ) dset <- matrix_frame(X, TR = 2, run_length = 50, event_table = edata) fit <- fmri_lm(onsets ~ hrf(condition), block = ~run, dataset = dset) coef_names(fit)
Returns estimated coefficients from a fitted fmri_lm model. Every
form of the result has one orientation: voxels x terms, with one
row per voxel (or latent component) and one column per coefficient or
contrast, and the term names on the column margin. This matches
stats(), standard_error() and
p_values(), so results from the accessor family line up
element-wise (for example coef(fit) / standard_error(fit)).
## S3 method for class 'fmri_lm' coef( object, type = c("betas", "contrasts"), include_baseline = FALSE, recon = FALSE, ... )## S3 method for class 'fmri_lm' coef( object, type = c("betas", "contrasts"), include_baseline = FALSE, recon = FALSE, ... )
object |
A fitted |
type |
|
include_baseline |
Logical. For |
recon |
Logical; used by latent-space fits to reconstruct voxel-space coefficients. Ignored for ordinary fits. |
... |
Unused. |
For type = "betas", a numeric matrix with one row per voxel
and one column per coefficient, columns named after the design matrix
(see coef_names()). For type = "contrasts", a tibble
with one row per voxel and one column per contrast.
Before 0.3.0 the default call (type = "betas",
include_baseline = FALSE) returned the transpose, terms x voxels,
while the other forms were voxels x terms. Code written for the old
default that indexed coef(fit)[term, ] must now use
coef(fit)[, term].
coef_names(), stats() (t-statistics),
standard_error(), p_values(),
coef_image().
X <- matrix(rnorm(50 * 3), 50, 3) edata <- data.frame( condition = factor(c("A", "B", "A", "B")), onsets = c(1, 12, 25, 38), run = 1 ) dset <- matrix_frame(X, TR = 2, run_length = 50, event_table = edata) fit <- fmri_lm(onsets ~ hrf(condition), block = ~run, dataset = dset) b <- coef(fit) # 3 voxels x 2 event coefficients dim(b) colnames(b) dim(coef(fit, include_baseline = TRUE))X <- matrix(rnorm(50 * 3), 50, 3) edata <- data.frame( condition = factor(c("A", "B", "A", "B")), onsets = c(1, 12, 25, 38), run = 1 ) dset <- matrix_frame(X, TR = 2, run_length = 50, event_table = edata) fit <- fmri_lm(onsets ~ hrf(condition), block = ~run, dataset = dset) b <- coef(fit) # 3 voxels x 2 event coefficients dim(b) colnames(b) dim(coef(fit, include_baseline = TRUE))
Extract Coefficients from Meta-Analysis
## S3 method for class 'fmri_meta' coef(object, ...)## S3 method for class 'fmri_meta' coef(object, ...)
object |
An fmri_meta object |
... |
Additional arguments |
Matrix of coefficients
toy_meta <- structure( list(coefficients = matrix(c(0.2, -0.1), nrow = 1)), class = "fmri_meta" ) coef(toy_meta)toy_meta <- structure( list(coefficients = matrix(c(0.2, -0.1), nrow = 1)), class = "fmri_meta" ) coef(toy_meta)
Globs the per-subject statistical maps written by reduce_write_results()
under dir, pairs them by subject, and builds a group_data() object
for fmri_meta(). When both beta and se maps are present the result carries
variance (suitable for inverse-variance meta-analysis); a beta-only result is
returned otherwise.
collect_results( dir, space = NULL, beta_desc = "beta", se_desc = "se", format = c("nifti") )collect_results( dir, space = NULL, beta_desc = "beta", se_desc = "se", format = c("nifti") )
dir |
Directory (searched recursively) containing the written maps. |
space |
Optional BIDS space label to filter on (e.g.
|
beta_desc, se_desc
|
The BIDS |
format |
Group-data backend format (currently |
A group_data object.
Pairing by BIDS subject requires fmrigds with BIDS-aware subject
keying. Multi-file "h5" aggregation is not yet supported upstream.
reduce_write_results(), run_jobs(), fmri_meta()
## Not run: tmpl <- fmri_template(onset ~ hrf(trial_type, contrasts = con), ~ run, reducer = reduce_write_results(format = "nifti", stats = c("beta", "se"), path = "study/glm")) run_jobs(instantiate(tmpl, manifest)) gd <- collect_results("study/glm", space = "MNI152NLin2009cAsym") fm <- fmri_meta(gd, ~ 1, method = "fe") ## End(Not run)## Not run: tmpl <- fmri_template(onset ~ hrf(trial_type, contrasts = con), ~ run, reducer = reduce_write_results(format = "nifti", stats = c("beta", "se"), path = "study/glm")) run_jobs(instantiate(tmpl, manifest)) gd <- collect_results("study/glm", space = "MNI152NLin2009cAsym") fm <- fmri_meta(gd, ~ 1, method = "fe") ## End(Not run)
DVARS measures the root mean square of the temporal derivative across all voxels, providing a single quality metric per volume. High DVARS indicates rapid signal changes, often associated with motion artifacts.
compute_dvars(Y, normalize = TRUE)compute_dvars(Y, normalize = TRUE)
Y |
Numeric matrix of fMRI data (time x voxels). |
normalize |
Logical. If TRUE, normalize by median DVARS. Default TRUE. |
Numeric vector of length nrow(Y) with DVARS values. The first timepoint is set to the median of subsequent values.
# Simulate fMRI data with a motion spike set.seed(123) Y <- matrix(rnorm(100 * 50), nrow = 100, ncol = 50) Y[50, ] <- Y[50, ] + 5 # Add spike at volume 50 dvars <- compute_dvars(Y) # dvars[50] will be elevated compared to other volumes # Identify suspect volumes (normalized DVARS > 1.5) suspect <- which(dvars > 1.5) cat("Suspect volumes:", suspect, "\n")# Simulate fMRI data with a motion spike set.seed(123) Y <- matrix(rnorm(100 * 50), nrow = 100, ncol = 50) Y[50, ] <- Y[50, ] + 5 # Add spike at volume 50 dvars <- compute_dvars(Y) # dvars[50] will be elevated compared to other volumes # Identify suspect volumes (normalized DVARS > 1.5) suspect <- which(dvars > 1.5) cat("Suspect volumes:", suspect, "\n")
A polished wrapper around the internal fast contrast engine. Accepts name-based contrast specifications and returns a consistent tibble by default.
compute_lm_contrasts( B, XtXinv, df, sigma = NULL, sigma2 = NULL, contrasts = NULL, t_contrasts = NULL, f_contrasts = NULL, columns = NULL, output = c("stacked", "list"), robust_weights = NULL, ar_order = 0, drop_failed = TRUE )compute_lm_contrasts( B, XtXinv, df, sigma = NULL, sigma2 = NULL, contrasts = NULL, t_contrasts = NULL, f_contrasts = NULL, columns = NULL, output = c("stacked", "list"), robust_weights = NULL, ar_order = 0, drop_failed = TRUE )
B |
Numeric matrix (p x V) of coefficients. |
XtXinv |
Numeric matrix (p x p), inverse cross-product of the design. |
df |
Residual degrees of freedom. |
sigma |
Optional numeric vector (length V) of residual std. dev.; ignored if |
sigma2 |
Optional numeric vector (length V) of residual variances. |
contrasts |
Optional named list mixing t- and F-contrasts; vectors are t, matrices are F. |
t_contrasts |
Optional named list of numeric vectors (t-contrasts). |
f_contrasts |
Optional named list of numeric matrices (F-contrasts). |
columns |
Optional character vector (length p) naming coefficients; used for name matching. |
output |
Either "stacked" (default; tibble) or "list" (raw list of tibbles). |
robust_weights |
Optional numeric vector of robust weights or NULL. Used only for df adjustments. |
ar_order |
Integer AR order; used only for effective df adjustments. |
drop_failed |
Logical; drop contrasts that fail validation (default TRUE). |
Inputs mirror a standard GLM contrast computation: provide B (p x V),
XtXinv (p x p), residual degrees of freedom df, and either sigma or
sigma2 (per-voxel noise). Contrasts can be specified as a single mixed
contrasts list or separately as t_contrasts (vectors) and f_contrasts
(matrices). When contrast vectors/matrices are named (for vectors) or have
column names (for matrices), names are matched to columns to determine the
appropriate design indices; otherwise lengths must match the number of
parameters p.
A tibble with rows per contrast (default) or a named list if output = "list".
Convenience wrapper that accepts design/data sufficient statistics, computes
betas and residual variance, and delegates to compute_lm_contrasts().
compute_lm_contrasts_from_suffstats( XtX, XtS, StS, df, sigma = NULL, sigma2 = NULL, contrasts = NULL, t_contrasts = NULL, f_contrasts = NULL, columns = NULL, output = c("stacked", "list"), robust_weights = NULL, ar_order = 0, drop_failed = TRUE )compute_lm_contrasts_from_suffstats( XtX, XtS, StS, df, sigma = NULL, sigma2 = NULL, contrasts = NULL, t_contrasts = NULL, f_contrasts = NULL, columns = NULL, output = c("stacked", "list"), robust_weights = NULL, ar_order = 0, drop_failed = TRUE )
XtX |
Numeric (p x p) cross-product of the design. |
XtS |
Numeric (p x V) cross-product of design with data. |
StS |
Numeric length-V vector of sum of squares per voxel. |
df |
Residual degrees of freedom. |
sigma |
Optional numeric vector (length V) of residual std. dev.; ignored if |
sigma2 |
Optional numeric vector (length V) of residual variances. |
contrasts |
Optional named list mixing t- and F-contrasts; vectors are t, matrices are F. |
t_contrasts |
Optional named list of numeric vectors (t-contrasts). |
f_contrasts |
Optional named list of numeric matrices (F-contrasts). |
columns |
Optional character vector (length p) naming coefficients; used for name matching. |
output |
Either "stacked" (default; tibble) or "list" (raw list of tibbles). |
robust_weights |
Optional numeric vector of robust weights or NULL. Used only for df adjustments. |
ar_order |
Integer AR order; used only for effective df adjustments. |
drop_failed |
Logical; drop contrasts that fail validation (default TRUE). |
A tibble with rows per contrast (default) or a named list if output = "list".
Mechanical execution specification
compute_spec( voxel_chunks = 10L, backend = c("matrix", "reference"), parallel = c("none", "voxels", "chunks"), progress = FALSE )compute_spec( voxel_chunks = 10L, backend = c("matrix", "reference"), parallel = c("none", "voxels", "chunks"), progress = FALSE )
voxel_chunks |
Number of voxel partitions for joint fitting. |
backend |
Matrix or formula/reference implementation. |
parallel |
Parallelization dimension. |
progress |
Display progress. |
Provide a contrast generic that dispatches on the first argument. Falls back to fmridesign::contrast for non-fmri_meta classes.
contrast(x, ...)contrast(x, ...)
x |
object |
... |
passed to methods |
A contrast object with computed contrast weights and statistics
meta <- fmrireg:::.demo_fmri_meta() contrast(meta, c("(Intercept)" = 1))meta <- fmrireg:::.demo_fmri_meta() contrast(meta, c("(Intercept)" = 1))
Create a contrast set
contrast_set(...)contrast_set(...)
... |
contrast specifications |
A list of class "contrast_set" containing the specified contrasts
cs <- contrast_set( fmridesign::pair_contrast(~condition == "A", ~condition == "B", name = "A_vs_B") )cs <- contrast_set( fmridesign::pair_contrast(~condition == "A", ~condition == "B", name = "A_vs_B") )
The contrast standard error is , using the
coefficient covariance that fmri_meta() stores with every fit. For
objects without a stored covariance only single-coefficient contrasts are
exact; other contrasts fall back to a diagonal approximation with a warning.
## S3 method for class 'fmri_meta' contrast(x, contrast, ...)## S3 method for class 'fmri_meta' contrast(x, contrast, ...)
x |
An fmri_meta object |
contrast |
Contrast specification. Can be:
|
... |
Additional arguments |
An fmri_meta_contrast object with contrast results
toy_meta <- structure( list( coefficients = matrix(c(0.3, 0.1), nrow = 1, dimnames = list(NULL, c("A", "B"))), se = matrix(c(0.05, 0.06), nrow = 1), robust = "none" ), class = "fmri_meta" ) contrast(toy_meta, c(1, -1))toy_meta <- structure( list( coefficients = matrix(c(0.3, 0.1), nrow = 1, dimnames = list(NULL, c("A", "B"))), se = matrix(c(0.05, 0.06), nrow = 1), robust = "none" ), class = "fmri_meta" ) contrast(toy_meta, c(1, -1))
Heatmap of pairwise correlations between the columns of an fmri_model's
combined event + baseline design matrix, grouped into panels by model term
(event terms, drift, nuisance).
By default (within_run = TRUE) each column is centred within each
run and the run-intercept columns are dropped before correlating. This
shows the collinearity that matters for estimation once run means are
modelled, and removes spurious correlations that zero-padding creates
between columns belonging to different runs. Cells for pairs of columns
that live in different runs are left blank.
## S3 method for class 'fmri_model' correlation_map( x, method = c("pearson", "spearman"), half_matrix = TRUE, absolute_limits = TRUE, within_run = TRUE, label_values = NULL, ... )## S3 method for class 'fmri_model' correlation_map( x, method = c("pearson", "spearman"), half_matrix = TRUE, absolute_limits = TRUE, within_run = TRUE, label_values = NULL, ... )
x |
An |
method |
Correlation method (e.g., "pearson", "spearman"). |
half_matrix |
Logical; if TRUE (default), display only the lower triangle. |
absolute_limits |
Logical; if TRUE, the colour scale spans -1 to 1; otherwise it spans plus or minus the largest absolute correlation shown. |
within_run |
Logical; centre columns within runs and drop run
intercepts before correlating (default |
label_values |
Logical or |
... |
Additional arguments passed to |
A ggplot2 object.
Creates spatial blocks from a 3D mask for use with spatial_fdr. This is useful for voxelwise analyses where you want to group nearby voxels together for more powerful multiple comparisons correction.
create_3d_blocks(mask, block_size = c(10, 10, 10), connectivity = 26)create_3d_blocks(mask, block_size = c(10, 10, 10), connectivity = 26)
mask |
Numeric 3D array or NeuroVol object defining the brain mask. Non-zero values indicate voxels to include. |
block_size |
Integer vector of length 3 specifying block dimensions in voxels (default: c(10, 10, 10)) |
connectivity |
Integer scalar; type of connectivity for neighbors: 6 (face connectivity) or 26 (face, edge, and corner connectivity). Default: 26 |
List with components:
group_id |
Integer vector of group IDs for each voxel in the mask |
neighbors |
List of length n_groups where element i contains integer vector of 1-based neighbor IDs for group i |
n_groups |
Integer scalar; total number of groups created |
block_size |
Block size used |
mask <- array(c(1L, 1L, 0L, 1L, 1L, 0L, 0L, 0L, 0L), dim = c(3, 3, 1)) blocks <- create_3d_blocks(mask, block_size = c(2, 2, 1)) blocks$n_groupsmask <- array(c(1L, 1L, 0L, 1L, 1L, 0L, 0L, 0L, 0L), dim = c(3, 3, 1)) blocks <- create_3d_blocks(mask, block_size = c(2, 2, 1)) blocks$n_groups
Helper function to create a design matrix from a benchmark dataset using a specified HRF. This is useful for testing different HRF assumptions against the ground truth.
create_design_matrix_from_benchmark( dataset_name, hrf, include_intercept = TRUE )create_design_matrix_from_benchmark( dataset_name, hrf, include_intercept = TRUE )
dataset_name |
Character string specifying which dataset to use |
hrf |
HRF object to use for convolution (e.g., fmrihrf::HRF_SPMG1) |
include_intercept |
Logical, whether to include an intercept column |
A matrix with the design matrix (time x conditions)
# Create design matrix using canonical HRF X <- create_design_matrix_from_benchmark("BM_Canonical_HighSNR", fmrihrf::HRF_SPMG1) # Test with a different HRF X_wrong <- create_design_matrix_from_benchmark("BM_Canonical_HighSNR", fmrihrf::HRF_SPMG2)# Create design matrix using canonical HRF X <- create_design_matrix_from_benchmark("BM_Canonical_HighSNR", fmrihrf::HRF_SPMG1) # Test with a different HRF X_wrong <- create_design_matrix_from_benchmark("BM_Canonical_HighSNR", fmrihrf::HRF_SPMG2)
A small, serializable recipe for a dataset: the name of a dataset constructor plus the arguments to call it with. For file-backed data the arguments are paths and run lengths (no voxel data), which keeps the enclosing fmri_job tiny and portable. The dataset is realized lazily on the worker.
dataset_spec(constructor, args = list(), source = c("file", "inline"))dataset_spec(constructor, args = list(), source = c("file", "inline"))
constructor |
Name of a frame constructor (a string): one of
|
args |
A named list of arguments passed to |
source |
Either |
An object of class dataset_spec.
fmri_job(), nifti_frame(), matrix_frame()
dataset_spec("nifti_frame", args = list(scans = c("run-1_bold.nii.gz", "run-2_bold.nii.gz"), TR = 2, run_length = c(200, 200)), source = "file")dataset_spec("nifti_frame", args = list(scans = c("run-1_bold.nii.gz", "run-2_bold.nii.gz"), TR = 2, run_length = c(200, 200)), source = "file")
Produces a single heatmap of all columns in the design matrix of an
fmri_model, which merges the event_model and baseline_model
regressors. Rows are scans (time runs downward); columns are regressors,
kept in model order and grouped into panels by term (event terms, drift,
intercept, nuisance). Run boundaries are marked with horizontal rules, and
the subtitle reports the rank and condition number of the design.
By default each column is divided by its maximum absolute value, so
regressors on very different scales (an intercept, a motion trace, a
convolved event train) are all legible on one colour scale. Set
scale_columns = FALSE to show raw values. Run-intercept columns
are drawn in grey, since they carry no information beyond run membership.
The condition number is computed after scaling every column to unit length, so it reflects collinearity rather than the units of the columns.
## S3 method for class 'fmri_model' design_map( x, block_separators = TRUE, rotate_x_text = TRUE, fill_midpoint = NULL, fill_limits = NULL, scale_columns = TRUE, ... )## S3 method for class 'fmri_model' design_map( x, block_separators = TRUE, rotate_x_text = TRUE, fill_midpoint = NULL, fill_limits = NULL, scale_columns = TRUE, ... )
x |
An |
block_separators |
Logical; if |
rotate_x_text |
Logical; if |
fill_midpoint |
Numeric or |
fill_limits |
Numeric vector of length 2 or |
scale_columns |
Logical; if |
... |
Additional arguments passed to |
A ggplot2 plot object.
Extract the design matrix from an fmri_lm object by delegating to its model component.
## S3 method for class 'fmri_lm' design_matrix(x, ...)## S3 method for class 'fmri_lm' design_matrix(x, ...)
x |
An fmri_lm object |
... |
Additional arguments passed to the design_matrix method for the model |
The design matrix from the fmri_lm object's model
Extract the combined design matrix from an fMRI model containing both event and baseline terms.
## S3 method for class 'fmri_model' design_matrix(x, blockid = NULL, ...)## S3 method for class 'fmri_model' design_matrix(x, blockid = NULL, ...)
x |
An fmri_model object |
blockid |
Optional numeric vector specifying which blocks/runs to include |
... |
Additional arguments (not used) |
A tibble containing the combined design matrix with event and baseline terms
fm <- fmrireg:::.demo_fmri_model() head(design_matrix(fm))fm <- fmrireg:::.demo_fmri_model() head(design_matrix(fm))
Generates an interactive Shiny app that plots the design matrix for a given fMRI model. The design matrix is first converted into a long-format tibble and then plotted over time, faceted by block. Several customization options allow you to adjust the title, axis labels, line thickness, color palette, and more.
design_plot( fmrimod, term_name = NULL, longnames = FALSE, plot_title = NULL, x_label = "Time (s)", y_label = "Amplitude", line_size = 1, color_palette = "viridis", facet_ncol = 2, theme_custom = ggplot2::theme_minimal(base_size = 15) + ggplot2::theme(panel.spacing = ggplot2::unit(1, "lines")), legend_threshold = 30, ... )design_plot( fmrimod, term_name = NULL, longnames = FALSE, plot_title = NULL, x_label = "Time (s)", y_label = "Amplitude", line_size = 1, color_palette = "viridis", facet_ncol = 2, theme_custom = ggplot2::theme_minimal(base_size = 15) + ggplot2::theme(panel.spacing = ggplot2::unit(1, "lines")), legend_threshold = 30, ... )
fmrimod |
An |
term_name |
Optional: Name of the term to plot. If |
longnames |
Logical; if TRUE, use long condition names in the legend. Default is FALSE. |
plot_title |
Optional plot title. If |
x_label |
Label for the x-axis. Default is "Time". |
y_label |
Label for the y-axis. Default is "Value". |
line_size |
Numeric; line thickness for the plot. Default is 1. |
color_palette |
Character; name of a ColorBrewer palette to use (e.g., "Set1"). Default is "Set1". |
facet_ncol |
Number of columns for facet_wrap. Default is 1. |
theme_custom |
A ggplot2 theme to apply. Default is |
legend_threshold |
Numeric; if the number of unique conditions exceeds this value, the legend is hidden. Default is 25. |
... |
Additional arguments passed to ggplot2::geom_line(). |
A Shiny app that displays the design plot.
if (interactive()) { ## --- Construct a sampling frame --- sframe <- fmrihrf::sampling_frame(blocklens = c(100, 100), TR = 2, precision = 0.5) ## --- Create a dummy event table --- set.seed(123) event_table <- data.frame( onset = seq(10, 190, length.out = 20), x = rnorm(20), y = rnorm(20), run = rep(1:2, each = 10) ) ## --- Construct a baseline model --- base_mod <- baseline_model(basis = "bs", degree = 3, sframe = sframe, intercept = "runwise") ## --- Construct an event model using a formula --- ev_mod <- event_model(x = onset ~ hrf(x) + hrf(y), data = event_table, block = ~ run, sampling_frame = sframe, drop_empty = TRUE, durations = rep(0, nrow(event_table))) ## --- Combine into an fMRI model --- fmri_mod <- fmri_model(ev_mod, base_mod) ## --- Launch the interactive design plot --- design_plot(fmrimod = fmri_mod, term_name = NULL, longnames = TRUE, plot_title = "fMRI Design Matrix", x_label = "Time (s)", y_label = "Signal", line_size = 1.5, color_palette = "Set2", facet_ncol = 1) }if (interactive()) { ## --- Construct a sampling frame --- sframe <- fmrihrf::sampling_frame(blocklens = c(100, 100), TR = 2, precision = 0.5) ## --- Create a dummy event table --- set.seed(123) event_table <- data.frame( onset = seq(10, 190, length.out = 20), x = rnorm(20), y = rnorm(20), run = rep(1:2, each = 10) ) ## --- Construct a baseline model --- base_mod <- baseline_model(basis = "bs", degree = 3, sframe = sframe, intercept = "runwise") ## --- Construct an event model using a formula --- ev_mod <- event_model(x = onset ~ hrf(x) + hrf(y), data = event_table, block = ~ run, sampling_frame = sframe, drop_empty = TRUE, durations = rep(0, nrow(event_table))) ## --- Combine into an fMRI model --- fmri_mod <- fmri_model(ev_mod, base_mod) ## --- Launch the interactive design plot --- design_plot(fmrimod = fmri_mod, term_name = NULL, longnames = TRUE, plot_title = "fMRI Design Matrix", x_label = "Time (s)", y_label = "Signal", line_size = 1.5, color_palette = "Set2", facet_ncol = 1) }
Transforms DVARS quality metrics into weights for weighted least squares fitting. Volumes with high DVARS receive lower weights, implementing "soft scrubbing" without hard censoring thresholds.
dvars_to_weights( dvars, method = c("inverse_squared", "soft_threshold", "tukey"), threshold = 1.5, steepness = 2 )dvars_to_weights( dvars, method = c("inverse_squared", "soft_threshold", "tukey"), threshold = 1.5, steepness = 2 )
dvars |
Numeric vector of DVARS values (from |
method |
Character. Weighting method: "inverse_squared" (default), "soft_threshold", or "tukey". |
threshold |
Numeric. For "soft_threshold", the DVARS value above which weights decay. Default 1.5 (1.5x median if normalized). |
steepness |
Numeric. For "soft_threshold", controls decay rate. Default 2. |
Numeric vector of weights in the interval from 0 to 1, same length as dvars.
set.seed(123) Y <- matrix(rnorm(100 * 50), nrow = 100, ncol = 50) Y[50, ] <- Y[50, ] + 5 dvars <- compute_dvars(Y) # Compare different weighting methods w_inv <- dvars_to_weights(dvars, method = "inverse_squared") w_soft <- dvars_to_weights(dvars, method = "soft_threshold") w_tukey <- dvars_to_weights(dvars, method = "tukey") # Check weight at the artifact volume cat("Weights at volume 50:\n") cat(" inverse_squared:", round(w_inv[50], 3), "\n") cat(" soft_threshold:", round(w_soft[50], 3), "\n") cat(" tukey:", round(w_tukey[50], 3), "\n")set.seed(123) Y <- matrix(rnorm(100 * 50), nrow = 100, ncol = 50) Y[50, ] <- Y[50, ] + 5 dvars <- compute_dvars(Y) # Compare different weighting methods w_inv <- dvars_to_weights(dvars, method = "inverse_squared") w_soft <- dvars_to_weights(dvars, method = "soft_threshold") w_tukey <- dvars_to_weights(dvars, method = "tukey") # Check weight at the artifact volume cat("Weights at volume 50:\n") cat(" inverse_squared:", round(w_inv[50], 3), "\n") cat(" soft_threshold:", round(w_soft[50], 3), "\n") cat(" tukey:", round(w_tukey[50], 3), "\n")
Returns a read-only description of a registered engine, including its
normalized capabilities, source ("builtin" vs "plugin"), aliases, and
dispatch strategy. This is intended for extension authors and diagnostic
tooling; it does not expose the underlying fit/preflight functions.
engine_spec(name)engine_spec(name)
name |
Engine name, such as |
A list of class fmrireg_engine_spec, or NULL if the engine is not
registered.
Returns the read-only specifications for all currently registered engines.
engine_specs()engine_specs()
A named list of fmrireg_engine_spec objects.
This function is deprecated. Use estimate_betas() instead.
estimate(...)estimate(...)
... |
Ignored. |
No return value; always errors with a deprecation message.
Estimate beta coefficients (regression parameters) from fMRI data using various methods. This function supports different estimation approaches for:
Single-trial beta estimation
Fixed and random effects
Various regularization techniques
Optional HRF estimation
This function estimates betas (regression coefficients) for fixed and random effects using various regression methods including mixed models, least squares, and PLS.
estimate_betas(x, ...) ## S3 method for class 'fmri_frame' estimate_betas( x, fixed = NULL, ran, block, method = c("mixed", "lss", "ols"), basemod = NULL, maxit = 1000, fracs = 0.5, progress = TRUE, ... )estimate_betas(x, ...) ## S3 method for class 'fmri_frame' estimate_betas( x, fixed = NULL, ran, block, method = c("mixed", "lss", "ols"), basemod = NULL, maxit = 1000, fracs = 0.5, progress = TRUE, ... )
x |
An |
... |
Additional arguments passed to the estimation method. |
fixed |
A formula specifying the fixed regressors that model constant effects (i.e., non-varying over trials). |
ran |
A formula specifying the random (trialwise) regressors that model single trial effects. |
block |
A formula specifying the block factor. |
method |
The regression method for estimating trialwise betas; one of "mixed", "lss", or "ols". |
basemod |
A |
maxit |
Maximum number of iterations for optimization methods (default: 1000). |
fracs |
Fraction of voxels used for prewhitening. |
progress |
Logical; show progress bar. |
This is a generic function whose fmri_frame method adapts to the frame's
feature space: volumetric frames (volume_space) return NeuroVec betas,
while matrix-format (index_space) and latent (basis_space) frames return
coefficient matrices.
Available estimation methods include:
Mixed-effects model using ridge/BLUP estimation
Rank-1 GLM with joint HRF estimation
Least-squares separate estimation
Partial least squares regression
Ordinary least squares
A list of class "fmri_betas" containing:
Fixed effect coefficients
Random (trial-wise) coefficients
Design matrix for random effects
Design matrix for fixed effects
Design matrix for baseline model
Additional components specific to the estimation method used
A list of class "fmri_betas" containing the following components:
betas_fixed: fixed effect betas (NeuroVec for volumetric frames, matrix otherwise).
betas_ran: random effect betas (NeuroVec for volumetric frames, matrix otherwise).
design_ran: Design matrix for random effects.
design_fixed: Design matrix for fixed effects.
design_base: Design matrix for baseline model.
basemod: Baseline model object.
fixed_model: Fixed effect model object.
ran_model: Random effect model object.
estimated_hrf: The estimated HRF vector (NULL for most methods).
Mumford, J. A., et al. (2012). Deconvolving BOLD activation in event-related designs for multivoxel pattern classification analyses. NeuroImage, 59(3), 2636-2643.
Pedregosa, F., et al. (2015). Data-driven HRF estimation for encoding and decoding models. NeuroImage, 104, 209-220.
matrix_frame, neurovec_frame, latent_frame
matrix_frame, baseline_model, event_model
# Create example data event_data <- data.frame( condition = factor(c("A", "B", "A", "B")), onset = c(1, 10, 20, 30), run = c(1, 1, 1, 1) ) # Create sampling frame and dataset sframe <- sampling_frame(blocklens = 100, TR = 2) dset <- matrix_frame( matrix(rnorm(100 * 2), 100, 2), TR = 2, run_length = 100, event_table = event_data ) # Estimate betas using mixed-effects model betas <- estimate_betas( dset, fixed = onset ~ hrf(condition), ran = onset ~ trialwise(), block = ~run, method = "mixed" ) ## Not run: facedes <- read.table(system.file("extdata", "face_design.txt", package = "fmrireg"), header=TRUE) facedes$frun <- factor(facedes$run) scans <- paste0("rscan0", 1:6, ".nii") dset <- nifti_frame(scans=scans, mask="mask.nii", TR=1.5, run_length=rep(436,6), event_table=facedes) fixed = onset ~ hrf(run) ran = onset ~ trialwise() block = ~ run betas <- estimate_betas(dset, fixed=fixed, ran=ran, block=block, method="mixed") ## End(Not run)# Create example data event_data <- data.frame( condition = factor(c("A", "B", "A", "B")), onset = c(1, 10, 20, 30), run = c(1, 1, 1, 1) ) # Create sampling frame and dataset sframe <- sampling_frame(blocklens = 100, TR = 2) dset <- matrix_frame( matrix(rnorm(100 * 2), 100, 2), TR = 2, run_length = 100, event_table = event_data ) # Estimate betas using mixed-effects model betas <- estimate_betas( dset, fixed = onset ~ hrf(condition), ran = onset ~ trialwise(), block = ~run, method = "mixed" ) ## Not run: facedes <- read.table(system.file("extdata", "face_design.txt", package = "fmrireg"), header=TRUE) facedes$frun <- factor(facedes$run) scans <- paste0("rscan0", 1:6, ".nii") dset <- nifti_frame(scans=scans, mask="mask.nii", TR=1.5, run_length=rep(436,6), event_table=facedes) fixed = onset ~ hrf(run) ran = onset ~ trialwise() block = ~ run betas <- estimate_betas(dset, fixed=fixed, ran=ran, block=block, method="mixed") ## End(Not run)
estimate_hrf() estimates condition-level HRF curves directly from an fMRI
dataset. It builds an event-aligned spline basis, removes baseline and fixed
nuisance effects once, and solves every voxel together with one penalized
multiresponse linear system. This replaces the former voxel-by-voxel GAM
implementation, whose predictor did not represent post-stimulus time.
estimate_hrf( form, fixed = NULL, block, dataset, bs = NULL, rsam = seq(0, 20, by = 1), basemod = NULL, k = 8L, fx = NULL, progress = FALSE, basis = c("bspline", "tent"), lambda = "gcv", lambda_grid = c(0, 10^seq(-4, 4, length.out = 25L)), ci_level = 0.95 )estimate_hrf( form, fixed = NULL, block, dataset, bs = NULL, rsam = seq(0, 20, by = 1), basemod = NULL, k = 8L, fx = NULL, progress = FALSE, basis = c("bspline", "tent"), lambda = "gcv", lambda_grid = c(0, 10^seq(-4, 4, length.out = 25L)), ci_level = 0.95 )
form |
A two-sided event-model formula containing one or more |
fixed |
Optional event-model formula whose design is treated as nuisance. |
block |
Formula identifying acquisition runs or blocks. |
dataset |
An |
bs |
Deprecated legacy basis selector. Values from the former GAM API are
mapped to |
rsam |
Strictly increasing, finite post-stimulus times beginning at zero. |
basemod |
Optional baseline model. The default is a constant baseline. |
k |
Number of free HRF basis functions per curve. At least four for the default cubic B-spline basis. |
fx |
Deprecated legacy selector. |
progress |
Show progress while scanning the GCV grid. |
basis |
Either |
lambda |
Non-negative smoothing strength or |
lambda_grid |
Candidate smoothing strengths used for GCV. |
ci_level |
Confidence level in |
A single smoothing parameter is shared across voxels. With lambda = "gcv",
it is selected by scale-normalized generalized cross-validation so high-
variance voxels do not dominate the choice. Curves are constrained to zero at
the beginning and end of rsam. Standard errors and confidence intervals use
the fitted penalized-linear-model covariance under independent, homoscedastic
time-point errors; they are not autocorrelation-robust.
An object of class fmri_hrf_estimate. Its estimate and
std.error arrays have dimensions time by curve by voxel. The object also
contains labeled curve metadata, basis coefficients, smoothing diagnostics,
and the designs used for estimation. Use tidy(), predict(), coef(), or
as.matrix() for common downstream representations.
fitted_hrf(), event_model(), baseline_model()
set.seed(18) n <- 80L events <- data.frame( onset = seq(6, 62, by = 8), condition = factor(rep(c("A", "B"), 4)), run = 1L ) dataset <- matrix_frame( matrix(rnorm(n * 2), nrow = n), TR = 1, run_length = n, event_table = events ) fit <- estimate_hrf( onset ~ hrf(condition), block = ~run, dataset = dataset, rsam = 0:12, k = 6, lambda = 1 ) tidy(fit, voxel = 1)set.seed(18) n <- 80L events <- data.frame( onset = seq(6, 62, by = 8), condition = factor(rep(c("A", "B"), 4)), run = 1L ) dataset <- matrix_frame( matrix(rnorm(n * 2), nrow = n), TR = 1, run_length = n, event_table = events ) fit <- estimate_hrf( onset ~ hrf(condition), block = ~run, dataset = dataset, rsam = 0:12, k = 6, lambda = 1 ) tidy(fit, voxel = 1)
Statistical fitting scope
estimation_spec( scope = c("joint", "runwise_meta"), run_pooling = "fixed_effect" )estimation_spec( scope = c("joint", "runwise_meta"), run_pooling = "fixed_effect" )
scope |
Fit a single joint model or fit runs separately and combine them by fixed-effect meta-analysis. |
run_pooling |
Runwise combination rule. Only fixed-effect pooling is currently supported. |
Helper function to evaluate the performance of beta estimation methods on benchmark datasets by comparing estimated betas to ground truth.
evaluate_method_performance( dataset_name, estimated_betas, method_name = "Unknown" )evaluate_method_performance( dataset_name, estimated_betas, method_name = "Unknown" )
dataset_name |
Character string specifying which dataset to use |
estimated_betas |
Matrix of estimated beta values (conditions x voxels) |
method_name |
Character string describing the method used |
A list with performance metrics
## Not run: # Load dataset and create design matrix dataset <- load_benchmark_dataset("BM_Canonical_HighSNR") X <- create_design_matrix_from_benchmark("BM_Canonical_HighSNR", fmrihrf::HRF_SPMG1) # Fit simple linear model betas <- solve(t(X) %*% X) %*% t(X) %*% dataset$Y_noisy # Evaluate performance performance <- evaluate_method_performance("BM_Canonical_HighSNR", betas[-1, ], # Remove intercept "OLS") ## End(Not run)## Not run: # Load dataset and create design matrix dataset <- load_benchmark_dataset("BM_Canonical_HighSNR") X <- create_design_matrix_from_benchmark("BM_Canonical_HighSNR", fmrihrf::HRF_SPMG1) # Fit simple linear model betas <- solve(t(X) %*% X) %*% t(X) %*% dataset$Y_noisy # Evaluate performance performance <- evaluate_method_performance("BM_Canonical_HighSNR", betas[-1, ], # Remove intercept "OLS") ## End(Not run)
Writes a serialized job manifest plus a minimal, backend-agnostic
run_one.R runner to dir. Deliberately emits no scheduler code:
drive it with whatever array system you have, passing the 1-based job index
as the first argument, e.g. Rscript run_one.R $SLURM_ARRAY_TASK_ID.
export_jobs( jobs, dir, overwrite = FALSE, setup = "library(fmrireg)", results_dir = "results" )export_jobs( jobs, dir, overwrite = FALSE, setup = "library(fmrireg)", results_dir = "results" )
jobs |
A single fmri_job or a list of them (e.g. from |
dir |
Output directory (created if needed). |
overwrite |
Overwrite an existing manifest. Default |
setup |
Character vector of R lines run at the top of |
results_dir |
Default output subdirectory written by |
Invisibly, a list describing what was written (dir,
manifest, runner, n, ids).
read_jobs(), run_jobs(), instantiate()
## Not run: jobs <- instantiate(tmpl, manifest) export_jobs(jobs, "study/jobs") # then, per array task: Rscript study/jobs/run_one.R $SLURM_ARRAY_TASK_ID ## End(Not run)## Not run: jobs <- instantiate(tmpl, manifest) export_jobs(jobs, "study/jobs") # then, per array task: Rscript study/jobs/run_one.R $SLURM_ARRAY_TASK_ID ## End(Not run)
Extract Data for Meta-Analysis from CSV
extract_csv_data(gd, roi = NULL, contrast = NULL)extract_csv_data(gd, roi = NULL, contrast = NULL)
gd |
A group_data_csv object |
roi |
Optional ROI name to extract |
contrast |
Optional contrast name to extract |
List with effect sizes and variances
gd <- fmrireg:::.demo_group_data_csv() extract_csv_data(gd, roi = "ROI1")gd <- fmrireg:::.demo_group_data_csv() extract_csv_data(gd, roi = "ROI1")
Extracts voxel timeseries from regions defined by a mask (e.g., WM/CSF). This is the typical input for soft subspace projection.
extract_nuisance_timeseries(dataset, mask, run = NULL)extract_nuisance_timeseries(dataset, mask, run = NULL)
dataset |
An fmri_frame object. |
mask |
A binary mask (logical vector or 3D array) indicating nuisance voxels, or a file path to a NIfTI mask. |
run |
Optional run index to extract data from a specific run. |
Matrix of nuisance timeseries (time x nuisance_voxels).
Generic function for fitting contrasts to model objects.
fit_contrasts(object, ...)fit_contrasts(object, ...)
object |
A fitted model object |
... |
Additional arguments passed to methods |
Contrast results (format depends on method)
This function calculates contrasts for a fitted linear model.
## Default S3 method: fit_contrasts(object, conmat, colind, se = TRUE, cov.unscaled = NULL, ...)## Default S3 method: fit_contrasts(object, conmat, colind, se = TRUE, cov.unscaled = NULL, ...)
object |
The fitted linear model object. |
conmat |
The contrast matrix or contrast vector. |
colind |
The subset column indices in the design associated with the contrast. |
se |
Whether to compute standard errors, t-statistics, and p-values (default: TRUE). |
cov.unscaled |
Optional unscaled coefficient covariance |
... |
Additional arguments (unused) |
A list containing the following elements:
conmat: Contrast matrix.
sigma: Residual standard error.
df.residual: Degrees of freedom for residuals.
estimate: Estimated contrasts.
se: Standard errors of the contrasts (if se = TRUE).
stat: t-statistics for the contrasts (if se = TRUE).
prob: Probabilities associated with the t-statistics (if se = TRUE).
stat_type: Type of the statistics calculated.
S3 method for computing contrasts on fitted fmri_lm objects.
## S3 method for class 'fmri_lm' fit_contrasts(object, contrasts, ...)## S3 method for class 'fmri_lm' fit_contrasts(object, contrasts, ...)
object |
An fmri_lm object |
contrasts |
A list of contrast specifications |
... |
Additional arguments (unused) |
Multi-run fits estimated run by run (runwise_meta) with temporal
noise models, robust weighting, volume weights or projection whiten or
weight each run differently, so their pooled betas have no single stored
covariance; post-hoc contrasts are refused for them. Specify such
contrasts in the model formula, or fit with
estimation_spec(scope = "joint").
A list of contrast results
Computes OLS/GLS-equivalent estimates from cross-products without materializing
the transformed series. Engines can stream XtX, XtS, and StS and call this
to obtain a standard fmri_lm object. AR/robust are not estimated here; this is
a low-level helper for OLS-equivalent inference from suffstats.
fit_glm_from_suffstats( model, XtX, XtS, StS, df, cfg = NULL, dataset = NULL, strategy = "external", engine = "external" )fit_glm_from_suffstats( model, XtX, XtS, StS, df, cfg = NULL, dataset = NULL, strategy = "external", engine = "external" )
model |
An |
XtX |
p×p cross-product of the design matrix. |
XtS |
p×V cross-product of design with data. |
StS |
length-V vector of sum of squares per voxel. |
df |
Residual degrees of freedom. |
cfg |
Optional |
dataset |
Optional dataset whose feature axis describes the columns of
|
strategy |
Character label for the returned object. |
engine |
Character label for the returned object. |
An object of class fmri_lm.
This helper performs IID, non-robust OLS inference. It rejects controls that
request AR or robust fitting; use fit_glm_with_config() for those models.
fit_glm_on_transformed_series( model, Y, cfg = NULL, dataset = NULL, strategy = "external", engine = "external" )fit_glm_on_transformed_series( model, Y, cfg = NULL, dataset = NULL, strategy = "external", engine = "external" )
model |
An |
Y |
Numeric matrix with |
cfg |
Optional |
dataset |
Optional dataset whose feature axis describes the columns of
|
strategy |
Character label recorded on the returned object. Defaults to "external". |
engine |
Character label indicating the engine that produced the fit. Defaults to "external". |
An object of class fmri_lm.
Runs the same integrated solver used by fmri_lm, honoring AR/robust options
from cfg, but on an externally provided response matrix Y (T×V). This is
intended for engines that transform the time-series before inference.
fit_glm_with_config( model, Y, cfg = NULL, dataset = NULL, strategy = "external", engine = "external" )fit_glm_with_config( model, Y, cfg = NULL, dataset = NULL, strategy = "external", engine = "external" )
model |
An |
Y |
Numeric matrix with |
cfg |
Optional |
dataset |
Optional dataset whose feature axis describes the columns of
|
strategy |
Character label recorded on the returned object. Defaults to "external". |
engine |
Character label indicating the engine that produced the fit. Defaults to "external". |
An object of class fmri_lm.
Compute and return the fitted hemodynamic response function (HRF) for a model object. The HRF represents the expected BOLD response to neural activity. For models with multiple basis functions, this returns the combined HRF shape.
fitted_hrf(x, sample_at, ...)fitted_hrf(x, sample_at, ...)
x |
An object for which the fitted HRF should be computed |
sample_at |
A vector of time points at which the HRF should be sampled |
... |
Additional arguments passed to methods |
This generic function computes the fitted hemodynamic response function (HRF) for an object. The method needs to be implemented for specific object types.
A numeric vector containing the fitted HRF values at the requested time points
fmrihrf::HRF_SPMG1(), fmri_lm()
# Create a simple dataset with two conditions X <- matrix(rnorm(100 * 100), 100, 100) # 100 timepoints, 100 voxels event_data <- data.frame( condition = factor(c("A", "B", "A", "B")), onsets = c(1, 25, 50, 75), run = c(1, 1, 1, 1) ) # Create dataset and sampling frame dset <- matrix_frame(X, TR = 2, run_length = 100, event_table = event_data) sframe <- sampling_frame(blocklens = 100, TR = 2) # Create event model with canonical HRF evmodel <- event_model( onsets ~ hrf(condition), data = event_data, block = ~run, sampling_frame = sframe ) # Fit model fit <- fmri_lm( onsets ~ hrf(condition), block = ~run, dataset = dset ) # Get fitted HRF at specific timepoints times <- seq(0, 20, by = 0.5) # Sample from 0-20s every 0.5s hrf_values <- fitted_hrf(fit, sample_at = times)# Create a simple dataset with two conditions X <- matrix(rnorm(100 * 100), 100, 100) # 100 timepoints, 100 voxels event_data <- data.frame( condition = factor(c("A", "B", "A", "B")), onsets = c(1, 25, 50, 75), run = c(1, 1, 1, 1) ) # Create dataset and sampling frame dset <- matrix_frame(X, TR = 2, run_length = 100, event_table = event_data) sframe <- sampling_frame(blocklens = 100, TR = 2) # Create event model with canonical HRF evmodel <- event_model( onsets ~ hrf(condition), data = event_data, block = ~run, sampling_frame = sframe ) # Fit model fit <- fmri_lm( onsets ~ hrf(condition), block = ~run, dataset = dset ) # Get fitted HRF at specific timepoints times <- seq(0, 20, by = 0.5) # Sample from 0-20s every 0.5s hrf_values <- fitted_hrf(fit, sample_at = times)
This method computes the fitted hemodynamic response functions (HRFs) for an fmri_lm object.
## S3 method for class 'fmri_lm' fitted_hrf(x, sample_at = seq(0, 24, by = 1), ...)## S3 method for class 'fmri_lm' fitted_hrf(x, sample_at = seq(0, 24, by = 1), ...)
x |
An |
sample_at |
A numeric vector of time points at which the HRFs should be sampled. Default is |
... |
Additional arguments (currently unused). |
A list where each element corresponds to an event term in the fmri_lm object. Each element contains:
predA matrix of predicted HRF values.
designA tibble containing the design matrix for the HRFs.
Reverses the sign of coefficient-like outputs in a fit object. Useful for switching between A-B and B-A conventions.
flip_sign(fit, coef = NULL)flip_sign(fit, coef = NULL)
fit |
An fmri_ttest_fit or similar object |
coef |
Character vector of coefficient names to flip (default: all) |
Modified fit object with flipped signs
A list of simulated datasets used for testing the package.
A list with elements BM_Canonical_HighSNR,
BM_Canonical_LowSNR,
BM_HRF_Variability_AcrossVoxels,
BM_Trial_Amplitude_Variability,
BM_Complex_Realistic, and metadata.
Each element contains simulated BOLD data and ground truth information.
Generated by data-raw/generate_benchmark_datasets.R
Methods for smooth HRF estimates
## S3 method for class 'fmri_hrf_estimate' print(x, ...) ## S3 method for class 'fmri_hrf_estimate' coef(object, ...) ## S3 method for class 'fmri_hrf_estimate' as.matrix(x, curve = NULL, what = c("estimate", "std.error"), ...) ## S3 method for class 'fmri_hrf_estimate' predict(object, newdata = object$time, se.fit = FALSE, ...)## S3 method for class 'fmri_hrf_estimate' print(x, ...) ## S3 method for class 'fmri_hrf_estimate' coef(object, ...) ## S3 method for class 'fmri_hrf_estimate' as.matrix(x, curve = NULL, what = c("estimate", "std.error"), ...) ## S3 method for class 'fmri_hrf_estimate' predict(object, newdata = object$time, se.fit = FALSE, ...)
x, object
|
An |
... |
Unused. |
curve |
A curve name or one-based curve index. Required by |
what |
Return estimated curves or their standard errors. |
newdata |
Numeric post-stimulus times within the fitted span. |
se.fit |
Return prediction standard errors with fitted curves. |
Binds an fmri_template to one subject's dataset_spec. The result is a
fully serializable recipe (no voxel data, no captured environment) that
run_job() turns into a fitted model and reduced output. Build jobs with
instantiate() rather than by hand in the common case.
fmri_job(id, template, dataset_spec, meta = list(), nuisance = NULL)fmri_job(id, template, dataset_spec, meta = list(), nuisance = NULL)
id |
A unique job identifier (e.g. a subject label). |
template |
An fmri_template. |
dataset_spec |
A dataset_spec describing this subject's data. |
meta |
Optional named list of metadata used for output keying
(e.g. |
nuisance |
Optional per-subject nuisance / confound regressors fed to the
baseline model at run time: |
An object of class fmri_job.
instantiate(), run_job(), export_jobs()
tmpl <- fmri_template(onset ~ hrf(condition), ~ run) ds <- dataset_spec("nifti_frame", args = list(scans = "run-1_bold.nii.gz", TR = 2, run_length = 200), source = "file") fmri_job("sub-01", tmpl, ds, meta = list(subject = "01"))tmpl <- fmri_template(onset ~ hrf(condition), ~ run) ds <- dataset_spec("nifti_frame", args = list(scans = "run-1_bold.nii.gz", TR = 2, run_length = 200), source = "file") fmri_job("sub-01", tmpl, ds, meta = list(subject = "01"))
This function estimates a regression model for fMRI data using a latent component dataset.
The dataset must be an fmri_frame whose feature space is a
fmridataset::basis_space (see latent_frame()), so that the assay holds
component scores and the space carries the spatial loadings.
fmri_latent_lm( formula, block, baseline_model = NULL, dataset, durations, drop_empty = TRUE, robust = FALSE, autocor = c("none", "auto", "ar1", "ar2", "arma"), bootstrap = FALSE, nboot = 1000, ... )fmri_latent_lm( formula, block, baseline_model = NULL, dataset, durations, drop_empty = TRUE, robust = FALSE, autocor = c("none", "auto", "ar1", "ar2", "arma"), bootstrap = FALSE, nboot = 1000, ... )
formula |
A formula specifying the regression model. |
block |
A factor indicating the block structure of the data. |
baseline_model |
An optional baseline model. |
dataset |
A latent |
durations |
The duration of events in the dataset. |
drop_empty |
Whether to drop empty events from the model. Default is TRUE. |
robust |
Must be FALSE. Robust weights are voxel-specific and cannot be computed on component scores, so any robust request is an error. |
autocor |
The autocorrelation model. 'none' (OLS), 'ar1' or 'ar2'. AR models use one set of AR coefficients per run shared across voxels, estimated so that the fit equals a dense fit of the reconstructed voxel data. 'auto' and 'arma' are not implemented and are an error. |
bootstrap |
Must be FALSE; bootstrapping is not implemented for latent datasets and TRUE is an error. |
nboot |
Unused (see |
... |
Additional arguments. |
An object of class 'fmri_latent_lm' containing the regression model and dataset.
This method is currently experimental.
# Estimate the fMRI regression model using the latent dataset #result <- fmri_latent_lm(formula = formula, block = block, dataset = dset, # durations = NULL, drop_empty = TRUE, robust = FALSE) # Print the result #print(result)# Estimate the fMRI regression model using the latent dataset #result <- fmri_latent_lm(formula = formula, block = block, dataset = dset, # durations = NULL, drop_empty = TRUE, robust = FALSE) # Print the result #print(result)
fmri_lm is a generic for fitting fMRI regression models. The
default interface accepts a model formula and dataset. An
alternative method can be used with a preconstructed
fmri_model object that already contains the design and data.
The public fitting boundary separates statistical choices in fmri_lm_control()
from mechanical execution choices in compute_spec(). Optional engines use
the same control object and receive their private arguments through
engine_args.
fmri_lm(formula, ...) ## S3 method for class 'formula' fmri_lm( formula, block, baseline_model = NULL, dataset, durations = 0, drop_empty = TRUE, control = fmri_lm_control(estimation = estimation_spec("runwise_meta")), compute = compute_spec(), engine = NULL, engine_args = list(), lowrank = NULL, ... ) ## S3 method for class 'fmri_model' fmri_lm( formula, dataset = NULL, control = fmri_lm_control(estimation = estimation_spec("runwise_meta")), compute = compute_spec(), engine = NULL, engine_args = list(), lowrank = NULL, ... )fmri_lm(formula, ...) ## S3 method for class 'formula' fmri_lm( formula, block, baseline_model = NULL, dataset, durations = 0, drop_empty = TRUE, control = fmri_lm_control(estimation = estimation_spec("runwise_meta")), compute = compute_spec(), engine = NULL, engine_args = list(), lowrank = NULL, ... ) ## S3 method for class 'fmri_model' fmri_lm( formula, dataset = NULL, control = fmri_lm_control(estimation = estimation_spec("runwise_meta")), compute = compute_spec(), engine = NULL, engine_args = list(), lowrank = NULL, ... )
formula |
A model formula describing the event structure or an
|
... |
Transitional legacy arguments retained for one compatibility
window. New code must express statistical choices in |
block |
Formula describing run/block structure. |
baseline_model |
Optional baseline/nuisance model. |
dataset |
An |
durations |
Event durations passed to model construction. |
drop_empty |
Remove empty factor levels during model construction. |
control |
A validated |
compute |
A validated |
engine |
Optional registered engine name. |
engine_args |
Named list passed only to |
lowrank |
Optional |
An object of class fmri_lm.
fmri_lm_control() is the single statistical configuration boundary for
fmri_lm(). Each section is a validated typed specification; mechanical
choices such as chunking and parallelism belong in compute_spec().
fmri_lm_control( estimation = NULL, noise = NULL, robust = NULL, variance = NULL, weights = NULL, projection = NULL, na_action = c("error", "propagate"), ... )fmri_lm_control( estimation = NULL, noise = NULL, robust = NULL, variance = NULL, weights = NULL, projection = NULL, na_action = c("error", "propagate"), ... )
estimation |
|
noise |
A |
robust |
|
variance |
|
weights |
|
projection |
|
na_action |
Either |
... |
Transitional legacy option lists retained for one compatibility window. New code should use the typed sections above. |
An object of class fmri_lm_control.
Performs voxelwise or ROI-based meta-analysis on group fMRI data using fixed-effects, random-effects, or robust methods. Supports meta-regression with covariates for group comparisons and other moderator analyses.
fmri_meta( data, formula = ~1, method = c("pm", "fe", "dl", "reml"), robust = c("none", "huber", "t"), weights = c("ivw", "equal", "custom"), weights_custom = NULL, combine = NULL, contrasts = NULL, return_cov = NULL, chunk_size = 10000, n_threads = getOption("fmrireg.num_threads", 0), verbose = TRUE )fmri_meta( data, formula = ~1, method = c("pm", "fe", "dl", "reml"), robust = c("none", "huber", "t"), weights = c("ivw", "equal", "custom"), weights_custom = NULL, combine = NULL, contrasts = NULL, return_cov = NULL, chunk_size = 10000, n_threads = getOption("fmrireg.num_threads", 0), verbose = TRUE )
data |
A group_data object created by |
formula |
Formula specifying the meta-regression model. Default is ~ 1 (intercept only). Use ~ 1 + group for group comparisons, or include continuous covariates. |
method |
Character string specifying the meta-analysis method:
|
robust |
Character string specifying robust estimation:
|
weights |
Character string specifying weighting scheme:
"equal" and "custom" weights are relative: the model is weighted
least squares, |
weights_custom |
Numeric vector or matrix of positive custom weights
(required if |
combine |
For t-statistic-only data, combination method ("stouffer",
"fisher", or "lancaster"). Stouffer combines z-scores and supports
equal or custom weighting (via |
contrasts |
Optional numeric vector or matrix specifying fit-time exact contrasts. If a vector is provided, its names must match the column names of the design matrix X. A matrix should have columns corresponding to predictors and rows corresponding to contrasts. |
return_cov |
Optional; kept for backward compatibility. The packed
upper-triangular Var(beta) per feature is now always returned under
|
chunk_size |
Number of voxels to process at once (default: 10000) |
n_threads |
Number of parallel threads to use. Defaults to fmrireg.num_threads option. |
verbose |
Logical. Print progress messages (default: TRUE) |
An fmri_meta object containing:
coefficients: Meta-regression coefficients
se: Standard errors
tau2: Between-study variance (for random-effects)
I2: I-squared heterogeneity statistic
Q: Cochran's Q statistic
model: Model specification
data: Input group_data object
## Not run: # Simple fixed-effects meta-analysis fit <- fmri_meta(gd, method = "fe") # Random-effects with group comparison fit <- fmri_meta(gd, formula = ~ 1 + group, method = "pm") # Robust meta-regression with continuous covariate fit <- fmri_meta(gd, formula = ~ 1 + age + sex, method = "reml", robust = "huber") # Stouffer's Z for t-statistics only fit <- fmri_meta(gd_tstat, combine = "stouffer") # Exact post-hoc contrasts by storing covariance fit_cov <- fmri_meta(gd, formula = ~ 1 + group, method = "pm", return_cov = "tri") con <- contrast(fit_cov, c("(Intercept)" = 0, groupB = 1)) # Exact fit-time contrast without storing covariance fit_con <- fmri_meta(gd, formula = ~ 1 + group, method = "pm", contrasts = c("(Intercept)" = 0, groupB = 1)) ## End(Not run)## Not run: # Simple fixed-effects meta-analysis fit <- fmri_meta(gd, method = "fe") # Random-effects with group comparison fit <- fmri_meta(gd, formula = ~ 1 + group, method = "pm") # Robust meta-regression with continuous covariate fit <- fmri_meta(gd, formula = ~ 1 + age + sex, method = "reml", robust = "huber") # Stouffer's Z for t-statistics only fit <- fmri_meta(gd_tstat, combine = "stouffer") # Exact post-hoc contrasts by storing covariance fit_cov <- fmri_meta(gd, formula = ~ 1 + group, method = "pm", return_cov = "tri") con <- contrast(fit_cov, c("(Intercept)" = 0, groupB = 1)) # Exact fit-time contrast without storing covariance fit_con <- fmri_meta(gd, formula = ~ 1 + group, method = "pm", contrasts = c("(Intercept)" = 0, groupB = 1)) ## End(Not run)
Low-level function that calls the C++ meta-analysis implementation. This is typically called internally by higher-level functions like fmri_meta().
fmri_meta_fit( Y, V, X, method = c("pm", "dl", "fe", "reml"), robust = c("none", "huber"), huber_c = 1.345, robust_iter = 2, n_threads = getOption("fmrireg.num_threads", 0) )fmri_meta_fit( Y, V, X, method = c("pm", "dl", "fe", "reml"), robust = c("none", "huber"), huber_c = 1.345, robust_iter = 2, n_threads = getOption("fmrireg.num_threads", 0) )
Y |
Numeric matrix of effect sizes (subjects x features) |
V |
Numeric matrix of variances (subjects x features) |
X |
Numeric matrix; design matrix (subjects x predictors), including intercept |
method |
Character scalar; meta-analysis method: "pm" (Paule-Mandel), "dl" (DerSimonian-Laird),
"fe" (fixed-effects). "reml" is not implemented: it is accepted for
backward compatibility, runs Paule-Mandel with a message, and the result
records |
robust |
Character scalar; robust estimation method: "none" or "huber" |
huber_c |
Numeric scalar; tuning constant for Huber M-estimator (default: 1.345). Smaller values provide more robust estimates but may reduce efficiency. |
robust_iter |
Integer scalar; number of IRLS iterations for robust estimation (default: 2) |
n_threads |
Integer scalar; number of OpenMP threads (0 = use all available) |
List with components:
beta |
Numeric matrix of coefficients (predictors x features) |
se |
Numeric matrix of standard errors (predictors x features) |
z |
Numeric matrix of z-scores (predictors x features) |
tau2 |
Numeric vector of between-study variance estimates |
Q_fe |
Numeric vector of Q statistics from fixed-effects model |
I2_fe |
Numeric vector of I² statistics from fixed-effects model |
df |
Numeric vector of degrees of freedom |
ok |
Logical vector indicating successful fits |
Low-level function that calls the C++ meta-analysis implementation and returns exact contrast statistics c' (X' W X)^(-1) c for provided contrasts.
fmri_meta_fit_contrasts( Y, V, X, Cmat, method = c("pm", "dl", "fe", "reml"), robust = c("none", "huber"), huber_c = 1.345, robust_iter = 2, n_threads = getOption("fmrireg.num_threads", 0) )fmri_meta_fit_contrasts( Y, V, X, Cmat, method = c("pm", "dl", "fe", "reml"), robust = c("none", "huber"), huber_c = 1.345, robust_iter = 2, n_threads = getOption("fmrireg.num_threads", 0) )
Y, V, X
|
See fmri_meta_fit |
Cmat |
Numeric matrix K x J (columns are contrasts over predictors) |
method |
Character scalar; meta-analysis method: "pm" (Paule-Mandel), "dl" (DerSimonian-Laird),
"fe" (fixed-effects). "reml" is not implemented: it is accepted for
backward compatibility, runs Paule-Mandel with a message, and the result
records |
robust |
Character scalar; robust estimation method: "none" or "huber" |
huber_c |
Numeric scalar; tuning constant for Huber M-estimator (default: 1.345). Smaller values provide more robust estimates but may reduce efficiency. |
robust_iter |
Integer scalar; number of IRLS iterations for robust estimation (default: 2) |
n_threads |
Integer scalar; number of OpenMP threads (0 = use all available) |
List with base meta outputs plus c_beta, c_se, c_z
Fit Meta-Analysis and return packed covariance per voxel
fmri_meta_fit_cov( Y, V, X, method = c("pm", "dl", "fe", "reml"), robust = c("none", "huber"), huber_c = 1.345, robust_iter = 2, n_threads = getOption("fmrireg.num_threads", 0) )fmri_meta_fit_cov( Y, V, X, method = c("pm", "dl", "fe", "reml"), robust = c("none", "huber"), huber_c = 1.345, robust_iter = 2, n_threads = getOption("fmrireg.num_threads", 0) )
Y |
Numeric matrix of effect sizes (subjects x features) |
V |
Numeric matrix of variances (subjects x features) |
X |
Numeric matrix; design matrix (subjects x predictors), including intercept |
method |
Character scalar; meta-analysis method: "pm" (Paule-Mandel), "dl" (DerSimonian-Laird),
"fe" (fixed-effects). "reml" is not implemented: it is accepted for
backward compatibility, runs Paule-Mandel with a message, and the result
records |
robust |
Character scalar; robust estimation method: "none" or "huber" |
huber_c |
Numeric scalar; tuning constant for Huber M-estimator (default: 1.345). Smaller values provide more robust estimates but may reduce efficiency. |
robust_iter |
Integer scalar; number of IRLS iterations for robust estimation (default: 2) |
n_threads |
Integer scalar; number of OpenMP threads (0 = use all available) |
List with base outputs and cov_tri (tsize x P) where tsize = K*(K+1)/2
Wrapper for meta-analysis that supports an optional voxelwise covariate. This extends the basic fmri_meta_fit to handle per-voxel covariates.
fmri_meta_fit_extended( Y, V, X, method = c("pm", "dl", "fe", "reml"), robust = c("none", "huber"), huber_c = 1.345, robust_iter = 2, voxelwise = NULL, center_voxelwise = TRUE, voxel_name = "voxel_cov", n_threads = getOption("fmrireg.num_threads", 0) )fmri_meta_fit_extended( Y, V, X, method = c("pm", "dl", "fe", "reml"), robust = c("none", "huber"), huber_c = 1.345, robust_iter = 2, voxelwise = NULL, center_voxelwise = TRUE, voxel_name = "voxel_cov", n_threads = getOption("fmrireg.num_threads", 0) )
Y |
Matrix of effect sizes (S x P) |
V |
Matrix of variances (S x P) |
X |
Design matrix (S x K) |
method |
Meta-analysis method |
robust |
Robust estimation method |
huber_c |
Huber tuning constant |
robust_iter |
Number of IRLS iterations |
voxelwise |
Optional voxelwise covariate matrix (S x P) |
center_voxelwise |
Logical; center voxelwise covariate per feature |
voxel_name |
Name for voxelwise coefficient |
n_threads |
Number of threads |
List with meta-analysis results
This function constructs an fMRI regression model consisting of an event model and a baseline model. The resulting model can be used for the analysis of fMRI data.
fmri_model(event_model, baseline_model, dataset)fmri_model(event_model, baseline_model, dataset)
event_model |
An object of class "event_model" representing the event-related part of the fMRI regression model. |
baseline_model |
An object of class "baseline_model" representing the baseline-related part of the fMRI regression model. |
dataset |
An |
An object of class fmri_model containing the event and baseline models along with the dataset.
event_model, baseline_model
Wrapper for OLS t-tests that supports an optional voxelwise covariate.
fmri_ols_fit( Y, X, voxelwise = NULL, center_voxelwise = TRUE, voxel_name = "voxel_cov" )fmri_ols_fit( Y, X, voxelwise = NULL, center_voxelwise = TRUE, voxel_name = "voxel_cov" )
Y |
Outcome matrix (S x P) |
X |
Design matrix (S x K) |
voxelwise |
Optional voxelwise covariate matrix (S x P) |
center_voxelwise |
Logical; center voxelwise covariate per feature |
voxel_name |
Name for voxelwise coefficient |
List with OLS results
This function fits a robust linear regression model for fMRI data analysis using the specified model formula, block structure, and dataset. The model can be fit using either a runwise or chunkwise data splitting strategy.
fmri_rlm( formula, block, baseline_model = NULL, dataset, durations = 0, drop_empty = TRUE, strategy = c("runwise", "chunkwise"), nchunks = 10, cor_struct = c("iid", "ar1", "ar2", "arp"), cor_iter = 1L, cor_global = FALSE, ar_p = NULL, ar1_exact_first = FALSE, robust_psi = c("huber", "bisquare"), robust_k_huber = 1.345, robust_c_tukey = 4.685, robust_max_iter = 2L, robust_scale_scope = c("run", "global", "voxel"), ... )fmri_rlm( formula, block, baseline_model = NULL, dataset, durations = 0, drop_empty = TRUE, strategy = c("runwise", "chunkwise"), nchunks = 10, cor_struct = c("iid", "ar1", "ar2", "arp"), cor_iter = 1L, cor_global = FALSE, ar_p = NULL, ar1_exact_first = FALSE, robust_psi = c("huber", "bisquare"), robust_k_huber = 1.345, robust_c_tukey = 4.685, robust_max_iter = 2L, robust_scale_scope = c("run", "global", "voxel"), ... )
formula |
A model formula describing the event structure or an
|
block |
Formula describing run/block structure. |
baseline_model |
Optional baseline/nuisance model. |
dataset |
An |
durations |
Event durations passed to model construction. |
drop_empty |
Remove empty factor levels during model construction. |
strategy |
Fit runs separately and combine them, or fit a joint model partitioned into voxel chunks. |
nchunks |
Number of data chunks when strategy is "chunkwise". Default is 10. |
cor_struct |
Correlation structure: "iid", "ar1", "ar2", or "arp". Default is "iid". |
cor_iter |
Number of iterations for AR parameter estimation. Default is 1. |
cor_global |
Whether to use global AR parameters. Default is FALSE. |
ar_p |
Order of autoregressive model for "arp" structure. Default is NULL. |
ar1_exact_first |
Whether to use exact AR(1) for first iteration. Default is FALSE. |
robust_psi |
Robust psi function: "huber" or "bisquare". Default is "huber". |
robust_k_huber |
Tuning constant for Huber's psi. Default is 1.345. |
robust_c_tukey |
Tuning constant for Tukey's bisquare. Default is 4.685. |
robust_max_iter |
Maximum iterations for robust fitting. Default is 2. |
robust_scale_scope |
Scope for robust scale estimation: "run", "global", or "voxel". Default is "run". |
... |
Additional arguments passed to fmri_lm |
A fitted robust linear regression model for fMRI data analysis.
etab <- data.frame(onset=c(1,30,15,25), fac=factor(c("A", "B", "A", "B")), run=c(1,1,2,2)) etab2 <- data.frame(onset=c(1,30,65,75), fac=factor(c("A", "B", "A", "B")), run=c(1,1,1,1)) mat <- matrix(rnorm(100*100), 100,100) dset <- matrix_frame(mat, TR=1, run_length=c(50,50),event_table=etab) dset2 <- matrix_frame(mat, TR=1, run_length=c(100),event_table=etab2) lm.1 <- fmri_rlm(onset ~ hrf(fac), block= ~ run, dataset=dset) lm.2 <- fmri_rlm(onset ~ hrf(fac), block= ~ run, dataset=dset2)etab <- data.frame(onset=c(1,30,15,25), fac=factor(c("A", "B", "A", "B")), run=c(1,1,2,2)) etab2 <- data.frame(onset=c(1,30,65,75), fac=factor(c("A", "B", "A", "B")), run=c(1,1,1,1)) mat <- matrix(rnorm(100*100), 100,100) dset <- matrix_frame(mat, TR=1, run_length=c(50,50),event_table=etab) dset2 <- matrix_frame(mat, TR=1, run_length=c(100),event_table=etab2) lm.1 <- fmri_rlm(onset ~ hrf(fac), block= ~ run, dataset=dset) lm.2 <- fmri_rlm(onset ~ hrf(fac), block= ~ run, dataset=dset2)
Captures a complete model specification that is constant across subjects.
Per-subject data (BOLD paths, event table, run lengths, confounds) is bound
later with instantiate() to produce serializable fmri_job recipes that
can be executed locally, in parallel via future, or on a cluster.
fmri_template( formula, block, baseline = baseline_spec(), durations = 0, contrasts = NULL, control = fmri_lm_control(estimation = estimation_spec("runwise_meta")), compute = compute_spec(voxel_chunks = 1L), engine = NULL, engine_args = list(), reducer = NULL )fmri_template( formula, block, baseline = baseline_spec(), durations = 0, contrasts = NULL, control = fmri_lm_control(estimation = estimation_spec("runwise_meta")), compute = compute_spec(voxel_chunks = 1L), engine = NULL, engine_args = list(), reducer = NULL )
formula |
Event model formula, e.g. |
block |
Block / run-structure formula, e.g. |
baseline |
A |
durations |
Event durations passed through to the event model. |
contrasts |
Optional contrast specification (e.g. from
|
control |
An |
compute |
A |
engine |
Optional fitting engine name (see |
engine_args |
Optional list of engine arguments. |
reducer |
Optional function |
The template is pure, serializable data: it holds no voxel data and (by
design) no captured execution environment. The reducer is the one
field that can accidentally capture state; see Reducer
serializability below.
An object of class fmri_template.
A reducer should be a top-level / package function (or a closure that
captures nothing), so it survives serialization to a worker node. A closure
that captures local variables will drag those bindings along when the job is
written with saveRDS(); fmri_template() warns in that case.
baseline_spec(), instantiate(), fmri_lm_control()
tmpl <- fmri_template(onset ~ hrf(condition), ~ run, baseline = baseline_spec(degree = 3)) tmpltmpl <- fmri_template(onset ~ hrf(condition), ~ run, baseline = baseline_spec(degree = 3)) tmpl
Performs group-level t-tests on fMRI data with support for one-sample, two-sample, paired, and ANCOVA designs. Provides both meta-analysis and classical t-test engines, mirroring AFNI 3dttest++ functionality within the fmrireg framework.
fmri_ttest( gd, formula = ~1, engine = c("auto", "meta", "classic", "welch"), paired = FALSE, mu0 = 0, contrast = NULL, mc = NULL, alpha = 0.05, sign = c("AminusB", "BminusA"), mask = NULL, voxelwise_cov = NULL, center_voxelwise = TRUE, voxel_name = "voxel_cov", weights = c("ivw", "equal", "custom"), weights_custom = NULL, combine = NULL )fmri_ttest( gd, formula = ~1, engine = c("auto", "meta", "classic", "welch"), paired = FALSE, mu0 = 0, contrast = NULL, mc = NULL, alpha = 0.05, sign = c("AminusB", "BminusA"), mask = NULL, voxelwise_cov = NULL, center_voxelwise = TRUE, voxel_name = "voxel_cov", weights = c("ivw", "equal", "custom"), weights_custom = NULL, combine = NULL )
gd |
A group_data object or data frame with subject data |
formula |
R formula for between-subjects design. Examples:
|
engine |
Character string; analysis engine to use:
|
paired |
Logical; if TRUE, perform paired t-test on differences (default: FALSE) |
mu0 |
Numeric scalar; constant to test against for one-sample tests (default: 0) |
contrast |
Optional named numeric vector for linear combination of coefficients |
mc |
Character string or NULL; multiple comparisons correction:
|
alpha |
Numeric scalar; FDR level if mc is not NULL (default: 0.05) |
sign |
Character string; sign convention for group differences:
|
mask |
Optional mask object or path |
voxelwise_cov |
Optional S x P matrix of voxelwise covariates |
center_voxelwise |
Logical; center voxelwise covariate per feature (default: TRUE) |
voxel_name |
Character string; name for voxelwise coefficient (default: "voxel_cov") |
weights |
Character string for meta-engine weighting: "ivw" (default,
uses provided SE), "equal" (equal weights), or "custom" (supply
|
weights_custom |
Numeric vector (length S) or matrix (S x P) of custom
weights when |
combine |
Optional. When using the meta engine with t-only inputs
(i.e., per-subject t-statistics and df), specify the t-combination
method: "stouffer", "fisher", or "lancaster". Passed through to
|
An fmri_ttest_fit object containing:
beta: Matrix of coefficients
se: Matrix of standard errors (if available)
t: Matrix of t-statistics (classic engine)
z: Matrix of z-scores
p: Matrix of p-values
df: Matrix of degrees of freedom
q: Matrix of FDR-adjusted p-values (if mc is used)
z_contrast: Vector of contrast z-scores (if contrast is used)
p_contrast: Vector of contrast p-values (if contrast is used)
## Not run: # One-sample t-test fit <- fmri_ttest(gd, formula = ~ 1) # Two-sample t-test fit <- fmri_ttest(gd, formula = ~ 1 + group) # Paired t-test fit <- fmri_ttest(gd_diff, formula = ~ 1, paired = TRUE) # ANCOVA with age covariate fit <- fmri_ttest(gd, formula = ~ 1 + group + age) # With spatial FDR correction fit <- fmri_ttest(gd, formula = ~ 1 + group, mc = "spatial_fdr", alpha = 0.05) ## End(Not run)## Not run: # One-sample t-test fit <- fmri_ttest(gd, formula = ~ 1) # Two-sample t-test fit <- fmri_ttest(gd, formula = ~ 1 + group) # Paired t-test fit <- fmri_ttest(gd_diff, formula = ~ 1, paired = TRUE) # ANCOVA with age covariate fit <- fmri_ttest(gd, formula = ~ 1 + group + age) # With spatial FDR correction fit <- fmri_ttest(gd, formula = ~ 1 + group, mc = "spatial_fdr", alpha = 0.05) ## End(Not run)
Runs the package-owned command line interface for fmrireg. The current
public command surface focuses on bundled benchmark dataset inspection through
the benchmark subcommands.
fmrireg_cli(args = commandArgs(trailingOnly = TRUE))fmrireg_cli(args = commandArgs(trailingOnly = TRUE))
args |
Character vector of command line arguments. |
Integer exit status. Returns 0 on success, 1 for
command-level validation failures, and 2 for usage or runtime errors.
fmrireg_cli(c("benchmark", "list")) fmrireg_cli(c( "benchmark", "summary", "--dataset", "BM_Canonical_HighSNR", "--json" ))fmrireg_cli(c("benchmark", "list")) fmrireg_cli(c( "benchmark", "summary", "--dataset", "BM_Canonical_HighSNR", "--json" ))
Uses bidser to discover preprocessed BOLD scans, events, confounds, TR,
and run lengths for each subject of a BIDS dataset, returning an
fmri_manifest (a list of per-subject bindings) suitable for
instantiate().
from_bids( proj, task, space = NULL, confounds = NULL, mask = NULL, desc = "preproc", subjects = NULL, ... )from_bids( proj, task, space = NULL, confounds = NULL, mask = NULL, desc = "preproc", subjects = NULL, ... )
proj |
A |
task |
BIDS task label. |
space |
BIDS space for preprocessed scans (e.g.
|
confounds |
Confound selection passed to |
mask |
A brain-mask path (length-1 character) applied to every subject's file-backed dataset. Required for file-backed BIDS data, which needs a mask (auto-resolution of the fmriprep brain mask is a planned enhancement). |
desc |
BIDS |
subjects |
Optional character vector of subject labels (bare ids, e.g.
|
... |
Forwarded to |
The design formula's variables must match the columns of the BIDS
events.tsv (e.g. trial_type); a run column is added for
the block structure. Confound columns are selected by confounds (a
character vector or bidser::confound_set() result) passed straight to
bidser::read_confounds().
An fmri_manifest: a list of bindings, one per subject.
as_manifest(), instantiate(), preflight()
## Not run: proj <- bidser::bids_project("study/", fmriprep = TRUE) mani <- from_bids(proj, task = "stroop", space = "MNI152NLin2009cAsym", confounds = bidser::confound_set("motion6"), mask = "study/derivatives/.../space-MNI..._desc-brain_mask.nii.gz") jobs <- instantiate(template, mani) ## End(Not run)## Not run: proj <- bidser::bids_project("study/", fmriprep = TRUE) mani <- from_bids(proj, task = "stroop", space = "MNI152NLin2009cAsym", confounds = bidser::confound_set("motion6"), mask = "study/derivatives/.../space-MNI..._desc-brain_mask.nii.gz") jobs <- instantiate(template, mani) ## End(Not run)
Returns a matrix N_cells x N_contrasts - row r is the design cell in
row r of des (for expand.grid() input, the first factor varies fastest,
which is also the column order of fmridesign crossed terms such as
hrf(A, B)),
columns are independent contrasts (difference-coded for the factors you ask
for, grand-mean for the rest). Suitable for tcrossprod(dm, C) or
lm.fit(design, y) followed by %*% coef in the usual way.
generate_interaction_contrast(des, factors) generate_main_effect_contrast(des, factor)generate_interaction_contrast(des, factors) generate_main_effect_contrast(des, factor)
des |
data.frame with one column per factor (must be |
factors |
character vector: which factor(s) get difference coding.
- |
factor |
Single factor name for the main effect. |
numeric matrix nrow = prod levels(f) , ncol = prod (Li - 1) for the chosen factors.
des <- expand.grid(Time = factor(1:4), Cond = factor(c("face","scene"))) # Main effect of Time (4-1 = 3 contrasts) M <- generate_main_effect_contrast(des, "Time") # Full TimexCond interaction ( (4-1)*(2-1) = 3 contrasts ) I <- generate_interaction_contrast(des, c("Time","Cond")) dim(I) # 8 rows (cells) x 3 columns (contrasts)des <- expand.grid(Time = factor(1:4), Cond = factor(c("face","scene"))) # Main effect of Time (4-1 = 3 contrasts) M <- generate_main_effect_contrast(des, "Time") # Full TimexCond interaction ( (4-1)*(2-1) = 3 contrasts ) I <- generate_interaction_contrast(des, c("Time","Cond")) dim(I) # 8 rows (cells) x 3 columns (contrasts)
Provides a detailed summary of a specific benchmark dataset including dimensions, experimental design, and ground truth information.
get_benchmark_summary(dataset_name)get_benchmark_summary(dataset_name)
dataset_name |
Character string specifying which dataset to summarize |
A list with summary information about the dataset
# Get summary of a specific dataset summary_info <- get_benchmark_summary("BM_Canonical_HighSNR") print(summary_info)# Get summary of a specific dataset summary_info <- get_benchmark_summary("BM_Canonical_HighSNR") print(summary_info)
Get Available Contrasts
get_contrasts(gd)get_contrasts(gd)
gd |
A group_data_csv object |
Character vector of contrast names
gd <- fmrireg:::.demo_group_data_csv() get_contrasts(gd)gd <- fmrireg:::.demo_group_data_csv() get_contrasts(gd)
Get Covariates
get_covariates(x)get_covariates(x)
x |
A group_data object |
Data frame of covariates or NULL
gd <- fmrireg:::.demo_group_data_csv() get_covariates(gd)gd <- fmrireg:::.demo_group_data_csv() get_covariates(gd)
Get Available ROIs
get_rois(gd)get_rois(gd)
gd |
A group_data_csv object |
Character vector of ROI names
gd <- fmrireg:::.demo_group_data_csv() get_rois(gd)gd <- fmrireg:::.demo_group_data_csv() get_rois(gd)
Get Subject IDs
get_subjects(x)get_subjects(x)
x |
A group_data object |
Character vector of unique subject IDs
gd <- fmrireg:::.demo_group_data_csv() get_subjects(gd)gd <- fmrireg:::.demo_group_data_csv() get_subjects(gd)
A convenience wrapper around estimate_betas for least squares separate (LSS) estimation.
This is primarily designed for single trial estimation, where each individual trial/event
gets its own separate beta estimate rather than averaging across trials of the same condition.
glm_lss( dataset, model_obj, basis_obj, basemod = NULL, block = ~1, use_cpp = FALSE, progress = TRUE, ... )glm_lss( dataset, model_obj, basis_obj, basemod = NULL, block = ~1, use_cpp = FALSE, progress = TRUE, ... )
dataset |
An |
model_obj |
An |
basis_obj |
An HRF basis object (e.g., from |
basemod |
A |
block |
A formula specifying the block factor (default: ~ 1 for single block) |
use_cpp |
Deprecated. The C++ implementation has been retired. This parameter is ignored; fmrilss is always used. |
progress |
Logical; show progress bar (default: TRUE) |
... |
Additional arguments passed to |
Primary Use Case - Single Trial Estimation:
Trial-wise beta estimation: Each trial gets its own beta coefficient
Single trial analysis: Useful for decoding, representational similarity analysis (RSA)
Trial-by-trial variability: Captures individual trial responses rather than condition averages
Avoiding trial averaging: Preserves trial-specific information that would be lost in standard GLM
Method Details: LSS (Least Squares Separate) fits a separate model for each trial, where the trial of interest gets its own regressor while all other trials of the same condition are modeled together. This approach avoids the collinearity issues that would arise from including separate regressors for every trial simultaneously.
For standard condition-level estimation (averaging trials within conditions), use glm_ols() instead.
A list of class "fmri_betas" containing the estimated trial-wise coefficients
estimate_betas for the underlying estimation function,
glm_ols for condition-level estimation
## Not run: # Create event model and data event_data <- data.frame( onset = c(10, 30, 50, 70), condition = factor(c("A", "B", "A", "B")), run = rep(1, 4) ) sframe <- fmrihrf::sampling_frame(blocklens = 100, TR = 2) model_obj <- event_model(onset ~ hrf(condition), data = event_data, block = ~ run, sampling_frame = sframe) # Create data matrix (100 timepoints, 10 voxels) Y <- matrix(rnorm(1000), 100, 10) # Create an fmri_frame with event table dset <- matrix_frame(Y, TR = 2, run_length = 100, event_table = event_data) # Fit with LSS - estimates separate beta for each individual trial fit <- glm_lss(dset, model_obj, fmrihrf::HRF_SPMG1) dim(fit$betas_ran) # 4 trials x 10 voxels (NOT averaged by condition) # This is useful for: # - Decoding analysis (predicting condition from single trial patterns) # - RSA (representational similarity analysis) # - Studying trial-by-trial variability ## End(Not run)## Not run: # Create event model and data event_data <- data.frame( onset = c(10, 30, 50, 70), condition = factor(c("A", "B", "A", "B")), run = rep(1, 4) ) sframe <- fmrihrf::sampling_frame(blocklens = 100, TR = 2) model_obj <- event_model(onset ~ hrf(condition), data = event_data, block = ~ run, sampling_frame = sframe) # Create data matrix (100 timepoints, 10 voxels) Y <- matrix(rnorm(1000), 100, 10) # Create an fmri_frame with event table dset <- matrix_frame(Y, TR = 2, run_length = 100, event_table = event_data) # Fit with LSS - estimates separate beta for each individual trial fit <- glm_lss(dset, model_obj, fmrihrf::HRF_SPMG1) dim(fit$betas_ran) # 4 trials x 10 voxels (NOT averaged by condition) # This is useful for: # - Decoding analysis (predicting condition from single trial patterns) # - RSA (representational similarity analysis) # - Studying trial-by-trial variability ## End(Not run)
A convenience wrapper around estimate_betas for ordinary least squares (OLS) estimation.
This function provides a simplified interface for fitting GLMs using OLS on matrix datasets.
glm_ols( dataset, model_obj, basis_obj, basemod = NULL, block = ~1, progress = TRUE, ... )glm_ols( dataset, model_obj, basis_obj, basemod = NULL, block = ~1, progress = TRUE, ... )
dataset |
An |
model_obj |
An |
basis_obj |
An HRF basis object (e.g., from |
basemod |
A |
block |
A formula specifying the block factor (default: ~ 1 for single block) |
progress |
Logical; show progress bar (default: TRUE) |
... |
Additional arguments passed to |
Use Cases:
Condition-level estimation: Estimates average responses for each experimental condition
General linear modeling: Standard GLM approach for group-level or condition-level effects
Multi-trial averaging: Combines trials of the same condition to estimate mean responses
For single-trial estimation where each trial gets its own beta estimate, use glm_lss() instead.
A list of class "fmri_betas" containing the estimated coefficients
estimate_betas for the underlying estimation function,
glm_lss for single trial estimation
## Not run: # Create event model and data event_data <- data.frame( onset = c(10, 30, 50, 70), condition = factor(c("A", "B", "A", "B")), run = rep(1, 4) ) sframe <- fmrihrf::sampling_frame(blocklens = 100, TR = 2) model_obj <- event_model(onset ~ hrf(condition), data = event_data, block = ~ run, sampling_frame = sframe) # Create data matrix (100 timepoints, 10 voxels) Y <- matrix(rnorm(1000), 100, 10) # Create an fmri_frame with event table dset <- matrix_frame(Y, TR = 2, run_length = 100, event_table = event_data) # Fit with OLS - estimates average response for each condition fit <- glm_ols(dset, model_obj, fmrihrf::HRF_SPMG1) dim(fit$betas_ran) # 2 conditions x 10 voxels ## End(Not run)## Not run: # Create event model and data event_data <- data.frame( onset = c(10, 30, 50, 70), condition = factor(c("A", "B", "A", "B")), run = rep(1, 4) ) sframe <- fmrihrf::sampling_frame(blocklens = 100, TR = 2) model_obj <- event_model(onset ~ hrf(condition), data = event_data, block = ~ run, sampling_frame = sframe) # Create data matrix (100 timepoints, 10 voxels) Y <- matrix(rnorm(1000), 100, 10) # Create an fmri_frame with event table dset <- matrix_frame(Y, TR = 2, run_length = 100, event_table = event_data) # Fit with OLS - estimates average response for each condition fit <- glm_ols(dset, model_obj, fmrihrf::HRF_SPMG1) dim(fit$betas_ran) # 2 conditions x 10 voxels ## End(Not run)
Generic constructor that creates a group dataset from various input formats
for use in group-level meta-analysis with fmri_meta.
group_data(data, format = c("auto", "h5", "nifti", "csv", "fmrilm"), ...)group_data(data, format = c("auto", "h5", "nifti", "csv", "fmrilm"), ...)
data |
Input data. Format depends on the
|
format |
Character string specifying the input format. One of "auto" (default), "h5", "nifti", "csv", or "fmrilm". If "auto", attempts to detect format. |
... |
Additional arguments passed to format-specific constructors |
A group_data object suitable for meta-analysis
## Not run: # From HDF5 files created by write_results.fmri_lm gd <- group_data( h5_paths, format = "h5", subjects = subject_ids, covariates = covariates_df, contrast = "FaceVsPlace" ) # From NIfTI files gd <- group_data( list(beta = beta_paths, se = se_paths), format = "nifti", subjects = subject_ids, mask = "group_mask.nii.gz" ) ## End(Not run)## Not run: # From HDF5 files created by write_results.fmri_lm gd <- group_data( h5_paths, format = "h5", subjects = subject_ids, covariates = covariates_df, contrast = "FaceVsPlace" ) # From NIfTI files gd <- group_data( list(beta = beta_paths, se = se_paths), format = "nifti", subjects = subject_ids, mask = "group_mask.nii.gz" ) ## End(Not run)
Creates a group dataset from tabular data containing pre-extracted statistics such as ROI means, effect sizes, and standard errors. This format is useful for ROI-based analyses or when working with summary statistics.
group_data_from_csv( data, effect_cols, subject_col = "subject", roi_col = NULL, contrast_col = NULL, covariate_cols = NULL, wide_format = FALSE )group_data_from_csv( data, effect_cols, subject_col = "subject", roi_col = NULL, contrast_col = NULL, covariate_cols = NULL, wide_format = FALSE )
data |
Either a path to a CSV file or a data frame containing the data |
effect_cols |
Named vector or list specifying column names for effect statistics. E.g., c(beta = "mean_activation", se = "std_error") or c(t = "t_stat", df = "df") |
subject_col |
Character string specifying the column containing subject IDs |
roi_col |
Character string specifying the column containing ROI names (optional) |
contrast_col |
Character string specifying the column containing contrast names (optional) |
covariate_cols |
Character vector of column names to use as covariates (optional) |
wide_format |
Logical. If TRUE, expects wide format with ROIs as columns (default: FALSE) |
A group_data_csv object
## Not run: # Long format: one row per subject-ROI combination gd <- group_data_from_csv( "roi_statistics.csv", effect_cols = c(beta = "mean_beta", se = "se_beta"), subject_col = "participant_id", roi_col = "roi_name", covariate_cols = c("age", "sex", "group") ) # Wide format: one row per subject, ROIs as columns gd <- group_data_from_csv( "subject_summary.csv", effect_cols = c(beta = "roi_"), # Prefix for ROI columns subject_col = "subject", wide_format = TRUE ) # From data frame with multiple contrasts df <- read.csv("contrast_results.csv") gd <- group_data_from_csv( df, effect_cols = c(beta = "estimate", se = "std_error", t = "t_value"), subject_col = "subject_id", contrast_col = "contrast_name", roi_col = "region" ) ## End(Not run)## Not run: # Long format: one row per subject-ROI combination gd <- group_data_from_csv( "roi_statistics.csv", effect_cols = c(beta = "mean_beta", se = "se_beta"), subject_col = "participant_id", roi_col = "roi_name", covariate_cols = c("age", "sex", "group") ) # Wide format: one row per subject, ROIs as columns gd <- group_data_from_csv( "subject_summary.csv", effect_cols = c(beta = "roi_"), # Prefix for ROI columns subject_col = "subject", wide_format = TRUE ) # From data frame with multiple contrasts df <- read.csv("contrast_results.csv") gd <- group_data_from_csv( df, effect_cols = c(beta = "estimate", se = "std_error", t = "t_value"), subject_col = "subject_id", contrast_col = "contrast_name", roi_col = "region" ) ## End(Not run)
Creates a group dataset directly from a list of fitted fmri_lm objects
group_data_from_fmrilm( lm_list, contrast = NULL, stat = c("beta", "se"), subjects = NULL, covariates = NULL )group_data_from_fmrilm( lm_list, contrast = NULL, stat = c("beta", "se"), subjects = NULL, covariates = NULL )
lm_list |
List of fmri_lm objects |
contrast |
Character string specifying which contrast to extract |
stat |
Character vector of statistics to extract |
subjects |
Character vector of subject identifiers |
covariates |
Data frame of subject-level covariates |
A group_data_h5 object (in-memory variant)
Creates a group dataset from HDF5 files produced by write_results.fmri_lm.
These files use the fmristore LabeledVolumeSet format for efficient storage of
multiple statistical maps.
group_data_from_h5( paths, subjects = NULL, covariates = NULL, mask = NULL, contrast = NULL, stat = c("beta", "se"), validate = TRUE )group_data_from_h5( paths, subjects = NULL, covariates = NULL, mask = NULL, contrast = NULL, stat = c("beta", "se"), validate = TRUE )
paths |
Character vector of HDF5 file paths, one per subject |
subjects |
Character vector of subject identifiers. If NULL, extracted from file paths. |
covariates |
Data frame of subject-level covariates (optional) |
mask |
Path to mask file or mask object (optional) |
contrast |
Character string specifying which contrast to extract (for multi-contrast files) |
stat |
Character vector of statistics to extract (e.g., c("beta", "se", "tstat")) |
validate |
Logical. Validate that all files exist and contain expected data (default: TRUE) |
A group_data_h5 object
## Not run: # Read HDF5 files from write_results.fmri_lm subjects <- data.frame( subject = sprintf("sub-%02d", 1:20), group = rep(c("young", "old"), each = 10), age = c(rnorm(10, 25, 3), rnorm(10, 70, 5)) ) h5_paths <- sprintf("derivatives/sub-%02d_task-nback_desc-GLMstatmap_bold.h5", 1:20) gd <- group_data_from_h5( h5_paths, subjects = subjects$subject, covariates = subjects[c("group", "age")], contrast = "FaceVsPlace", stat = c("beta", "se") ) ## End(Not run)## Not run: # Read HDF5 files from write_results.fmri_lm subjects <- data.frame( subject = sprintf("sub-%02d", 1:20), group = rep(c("young", "old"), each = 10), age = c(rnorm(10, 25, 3), rnorm(10, 70, 5)) ) h5_paths <- sprintf("derivatives/sub-%02d_task-nback_desc-GLMstatmap_bold.h5", 1:20) gd <- group_data_from_h5( h5_paths, subjects = subjects$subject, covariates = subjects[c("group", "age")], contrast = "FaceVsPlace", stat = c("beta", "se") ) ## End(Not run)
Creates a group dataset from NIfTI files containing effect sizes and their standard errors or variances. Supports various input configurations including beta/SE pairs, beta/variance pairs, or t-statistics with degrees of freedom.
group_data_from_nifti( beta_paths = NULL, se_paths = NULL, var_paths = NULL, t_paths = NULL, df = NULL, subjects = NULL, covariates = NULL, mask = NULL, target_space = NULL, validate = TRUE )group_data_from_nifti( beta_paths = NULL, se_paths = NULL, var_paths = NULL, t_paths = NULL, df = NULL, subjects = NULL, covariates = NULL, mask = NULL, target_space = NULL, validate = TRUE )
beta_paths |
Character vector of paths to beta/effect size NIfTI files |
se_paths |
Character vector of paths to standard error NIfTI files |
var_paths |
Character vector of paths to variance NIfTI files (alternative to se_paths) |
t_paths |
Character vector of paths to t-statistic NIfTI files |
df |
Degrees of freedom (scalar or vector). Required if using t_paths. |
subjects |
Character vector of subject identifiers. If NULL, extracted from file paths. |
covariates |
Data frame of subject-level covariates (optional) |
mask |
Path to mask NIfTI file or mask object (optional but recommended) |
target_space |
Path to template NIfTI for spatial alignment checking |
validate |
Logical. Validate that all files exist and have matching dimensions (default: TRUE) |
A group_data_nifti object
## Not run: # From FSL FEAT output (COPE and VARCOPE files) gd <- group_data_from_nifti( beta_paths = Sys.glob("feat_output/sub-*/cope1.nii.gz"), var_paths = Sys.glob("feat_output/sub-*/varcope1.nii.gz"), subjects = sprintf("sub-%02d", 1:20), mask = "group_mask.nii.gz" ) # From SPM contrast images gd <- group_data_from_nifti( beta_paths = Sys.glob("SPM/sub*/con_0001.nii"), se_paths = Sys.glob("SPM/sub*/se_0001.nii"), mask = "SPM/mask.nii" ) # From t-statistics only (for Stouffer's Z combination) gd <- group_data_from_nifti( t_paths = Sys.glob("stats/sub-*/tstat1.nii.gz"), df = 100, # Or vector of per-subject df mask = "group_mask.nii.gz" ) ## End(Not run)## Not run: # From FSL FEAT output (COPE and VARCOPE files) gd <- group_data_from_nifti( beta_paths = Sys.glob("feat_output/sub-*/cope1.nii.gz"), var_paths = Sys.glob("feat_output/sub-*/varcope1.nii.gz"), subjects = sprintf("sub-%02d", 1:20), mask = "group_mask.nii.gz" ) # From SPM contrast images gd <- group_data_from_nifti( beta_paths = Sys.glob("SPM/sub*/con_0001.nii"), se_paths = Sys.glob("SPM/sub*/se_0001.nii"), mask = "SPM/mask.nii" ) # From t-statistics only (for Stouffer's Z combination) gd <- group_data_from_nifti( t_paths = Sys.glob("stats/sub-*/tstat1.nii.gz"), df = 100, # Or vector of per-subject df mask = "group_mask.nii.gz" ) ## End(Not run)
This function computes a temporal similarity matrix from a series of hemodynamic response functions.
hrf_smoothing_kernel( len, TR = 2, form = onset ~ trialwise(), buffer_scans = 3L, normalise = TRUE, method = c("gram", "cosine") )hrf_smoothing_kernel( len, TR = 2, form = onset ~ trialwise(), buffer_scans = 3L, normalise = TRUE, method = c("gram", "cosine") )
len |
The number of scans. |
TR |
The repetition time (default is 2 seconds). |
form |
the |
buffer_scans |
The number of scans to buffer before and after the event. |
normalise |
Whether to normalise the kernel. |
method |
"gram" (inner products of the design rows) or "cosine"
(the same, scaled to unit diagonal). With |
a smoothing matrix
form <- onset ~ trialwise(basis="gaussian") sk <- hrf_smoothing_kernel(100, TR=1.5, form)form <- onset ~ trialwise(basis="gaussian") sk <- hrf_smoothing_kernel(100, TR=1.5, form)
Copies the package CLI wrapper from the installed exec/ directory into a
user-selected destination directory and, on Unix-like systems, marks the copied
file as executable.
install_cli(dest_dir, overwrite = FALSE, commands = NULL)install_cli(dest_dir, overwrite = FALSE, commands = NULL)
dest_dir |
Destination directory for installed command wrappers.
Required: there is no default, so nothing is written to the home
directory unless asked (e.g. |
overwrite |
Logical; overwrite an existing installed wrapper. |
commands |
Character vector of command names to install. Defaults to all available commands. |
Character vector of installed paths, returned invisibly.
## Not run: install_cli("~/.local/bin", overwrite = TRUE) ## End(Not run)## Not run: install_cli("~/.local/bin", overwrite = TRUE) ## End(Not run)
Binds an fmri_template to one or more per-subject data bindings, producing
serializable fmri_job recipes. No data is loaded: file-backed scans are kept
as paths and realized lazily by run_job().
instantiate(template, x, ...)instantiate(template, x, ...)
template |
An fmri_template. |
x |
A single binding (a named list with at least |
... |
Unused. |
A single fmri_job if x is one binding, otherwise a list of
fmri_jobs.
fmri_job(), run_job(), from_bids()
tmpl <- fmri_template(onset ~ hrf(condition), ~ run) b <- list(id = "sub-01", scans = matrix(rnorm(80 * 3), 80, 3), TR = 2, run_length = c(40, 40), events = data.frame(onset = c(5, 25, 45, 65), condition = factor(c("A", "B", "A", "B")), run = c(1, 1, 2, 2))) job <- instantiate(tmpl, b)tmpl <- fmri_template(onset ~ hrf(condition), ~ run) b <- list(id = "sub-01", scans = matrix(rnorm(80 * 3), 80, 3), TR = 2, run_length = c(40, 40), events = data.frame(onset = c(5, 25, 45, 65), condition = factor(c("A", "B", "A", "B")), run = c(1, 1, 2, 2))) job <- instantiate(tmpl, b)
A fmristore::LatentNeuroVec stores temporal scores (time by component)
and spatial loadings (voxel by component). This helper places the scores in
the assay and the loadings in a fmridataset::basis_space() whose parent
is the mask's volume_space, so fmrireg's latent engines can fit in
component space and reconstruct voxel maps through
fmridataset::basis_synthesis().
latent_frame( x, TR, run_length = NULL, event_table = NULL, censor = NULL, assay = "scores" )latent_frame( x, TR, run_length = NULL, event_table = NULL, censor = NULL, assay = "scores" )
x |
A |
TR |
Repetition time in seconds; one value, or one per run. |
run_length |
Volumes per run; defaults to a single run spanning all score rows. |
event_table, censor
|
As in |
assay |
Name of the assay. |
An fmri_frame over a synthesis-only basis_space.
matrix_frame(), fmridataset::basis_space_from_decoder()
Returns a summary of all available benchmark datasets with their descriptions.
list_benchmark_datasets()list_benchmark_datasets()
A data.frame with dataset names and descriptions
# See what benchmark datasets are available list_benchmark_datasets()# See what benchmark datasets are available list_benchmark_datasets()
This function provides easy access to the benchmark datasets included with the fmrireg package. These datasets are designed for testing HRF fitting, beta estimation, and other fMRI analysis methods.
load_benchmark_dataset(dataset_name = "BM_Canonical_HighSNR")load_benchmark_dataset(dataset_name = "BM_Canonical_HighSNR")
dataset_name |
Character string specifying which dataset to load. Options include:
|
A list containing the specified benchmark dataset(s) with the following components:
description: Text description of the dataset
Y_noisy: Matrix of noisy BOLD time series (time x voxels)
Y_clean: Matrix of clean BOLD time series (when available)
X_list_true_hrf: List of design matrices convolved with true HRF
true_hrf_parameters: Reconstructed true HRF object(s) plus
complete, versioned generation recipes
oracle_contract: Scope and tolerance of the valid stored
numerical oracle, including limitations for partial-truth datasets
event_onsets: Vector of event onset times
condition_labels: Vector of condition labels for each event
true_betas_condition: Matrix of true condition-level beta values
true_amplitudes_trial: Matrix of true trial-level amplitudes
TR: Repetition time
total_time: Total scan duration
run_length: Number of acquired time points
noise_parameters: Information about noise generation
simulation_seed: Random seed used for generation
target_snr: Target signal-to-noise ratio
# Load a specific dataset high_snr_data <- load_benchmark_dataset("BM_Canonical_HighSNR") # Get information about all available datasets metadata <- load_benchmark_dataset("metadata") # Load all datasets all_data <- load_benchmark_dataset("all") # Access the BOLD data Y <- high_snr_data$Y_noisy # Get event information onsets <- high_snr_data$event_onsets conditions <- high_snr_data$condition_labels# Load a specific dataset high_snr_data <- load_benchmark_dataset("BM_Canonical_HighSNR") # Get information about all available datasets metadata <- load_benchmark_dataset("metadata") # Load all datasets all_data <- load_benchmark_dataset("all") # Access the BOLD data Y <- high_snr_data$Y_noisy # Get event information onsets <- high_snr_data$event_onsets conditions <- high_snr_data$condition_labels
Control object to enable the optional sketched GLM engine.
lowrank_control( parcels = NULL, landmarks = NULL, k_neighbors = 16L, time_sketch = list(method = "gaussian", m = NULL), ncomp = NULL, noise_pcs = 0L )lowrank_control( parcels = NULL, landmarks = NULL, k_neighbors = 16L, time_sketch = list(method = "gaussian", m = NULL), ncomp = NULL, noise_pcs = 0L )
parcels |
Optional parceling, e.g., a neuroim2::ClusteredNeuroVol or integer vector (length = number of voxels in mask). |
landmarks |
Optional integer; number of landmark voxels for optional Nyström extension (NULL = off). |
k_neighbors |
Integer; k for k-NN in Nyström extension. |
time_sketch |
List controlling the temporal sketch, with elements
The former All sketches are normalised ( For speed, prefer For sketch-and-solve fits, Estimability is judged on the full (whitened) design with the same
pivoted-QR rule that exact fits use. Aliased coefficients are reported
as |
ncomp |
Optional integer; number of latent components within parcels (PCA). |
noise_pcs |
Integer; optional GLMdenoise-style PCs from low-R2 parcels. |
A list with class "lowrank_control".
Retrieves the estimated moving-average parameters from a fitted fMRI linear model that used ARMA or MA temporal-noise modeling.
ma_parameters(object, ...) ## S3 method for class 'fmri_lm' ma_parameters(object, scope = c("average", "per_run", "raw"), ...)ma_parameters(object, ...) ## S3 method for class 'fmri_lm' ma_parameters(object, scope = c("average", "per_run", "raw"), ...)
object |
An object of class |
... |
Additional arguments (currently unused). |
scope |
Character; |
Depending on scope, either a numeric vector of averaged MA
coefficients, a list of per-run coefficient vectors, or the raw stored
structure. Returns NULL when no MA modeling was performed.
fmrireg fits models to fmridataset::fmri_frame() objects. This helper
builds one from a time-by-feature matrix plus the run structure and event
table that the design machinery needs. Observation IDs are derived
deterministically as run-<r>-vol-<index>; feature IDs default to the
matrix column names when they are unique and to feature-<j> otherwise.
matrix_frame( datamat, TR, run_length, event_table = NULL, censor = NULL, feature_ids = NULL, assay = "signal" )matrix_frame( datamat, TR, run_length, event_table = NULL, censor = NULL, feature_ids = NULL, assay = "signal" )
datamat |
A numeric matrix with one row per acquired volume and one column per feature (voxel, vertex, parcel, or component). |
TR |
Repetition time in seconds; one value, or one per run. |
run_length |
Integer vector giving the number of volumes in each run.
Must sum to |
event_table |
Optional data frame of events (onsets, conditions, and
the block variable). It is stored on the frame as a keyed
|
censor |
Optional censoring indicator with one entry per volume: a
logical vector, a 0/1 vector, or a vector of 1-based volume indices.
Stored as the logical |
feature_ids |
Optional stable feature IDs (one per column). |
assay |
Name of the assay holding |
An fmri_frame whose feature space is an
fmridataset::index_space().
neurovec_frame(), nifti_frame(), latent_frame(),
fmridataset::fmri_frame()
Y <- matrix(rnorm(80 * 3), 80, 3) events <- data.frame(onset = c(5, 25, 45, 65), condition = factor(c("A", "B", "A", "B")), run = c(1, 1, 2, 2)) frame <- matrix_frame(Y, TR = 2, run_length = c(40, 40), event_table = events) fmridataset::temporal_schema(frame)$run_lengthsY <- matrix(rnorm(80 * 3), 80, 3) events <- data.frame(onset = c(5, 25, 45, 65), condition = factor(c("A", "B", "A", "B")), run = c(1, 1, 2, 2)) frame <- matrix_frame(Y, TR = 2, run_length = c(40, 40), event_table = events) fmridataset::temporal_schema(frame)$run_lengths
Computes the effective sample size based on the heterogeneity estimate. This is useful for understanding the impact of between-study heterogeneity.
meta_effective_n(v, tau2)meta_effective_n(v, tau2)
v |
Numeric vector of within-study variances |
tau2 |
Numeric scalar; between-study variance (tau-squared) |
Numeric scalar; effective sample size
meta_effective_n(v = rep(0.05, 3), tau2 = 0.01)meta_effective_n(v = rep(0.05, 3), tau2 = 0.01)
Meta-regression with ONE voxelwise covariate
meta_fit_vcov_cpp( Y, V, X, C, method, robust, huber_c = 1.345, robust_iter = 2L, n_threads = 0L )meta_fit_vcov_cpp( Y, V, X, C, method, robust, huber_c = 1.345, robust_iter = 2L, n_threads = 0L )
Y |
S x P matrix of effect sizes |
V |
S x P matrix of variances |
X |
S x K design matrix |
C |
S x P matrix of voxelwise covariates |
method |
Meta-analysis method |
robust |
Robust estimation method |
huber_c |
Huber tuning constant |
robust_iter |
Number of IRLS iterations |
n_threads |
Number of OpenMP threads |
List with results
This function solves a mixed model using Rcpp and roptim for optimization. It estimates variance components in a mixed model, potentially speeding up computations compared to the pure R implementation.
mixed_solve_cpp( y, Z = NULL, K = NULL, X = NULL, method = "REML", bounds = c(1e-09, 1e+09), SE = FALSE, return_Hinv = FALSE )mixed_solve_cpp( y, Z = NULL, K = NULL, X = NULL, method = "REML", bounds = c(1e-09, 1e+09), SE = FALSE, return_Hinv = FALSE )
y |
Response vector. |
Z |
Design matrix for random effects (default: identity matrix of size n). |
K |
Kinship matrix (default: NULL). |
X |
Design matrix for fixed effects (default: vector of ones). |
method |
Optimization method, either "REML" or "ML" (default: "REML"). |
bounds |
Bounds for the optimizer (default: c(1e-9, 1e9)). |
SE |
Logical, whether to return standard errors (default: FALSE). |
return_Hinv |
Logical, whether to return the inverse of H (default: FALSE). |
A list containing:
Vu |
Estimated variance component for random effects. |
Ve |
Estimated variance component for residuals. |
beta |
Estimated fixed effects coefficients. |
u |
Estimated random effects coefficients. |
LL |
Log-likelihood of the model. |
beta.SE |
Standard errors of fixed effects coefficients (if SE = TRUE). |
u.SE |
Standard errors of random effects coefficients (if SE = TRUE). |
Hinv |
Inverse of H (if return_Hinv = TRUE). |
## Not run: # Example usage with random data set.seed(123) n <- 100 y <- rnorm(n) Z <- matrix(rnorm(n * 5), n, 5) K <- diag(5) X <- matrix(1, n, 1) result <- mixed_solve_cpp(y, Z, K, X) ## End(Not run)## Not run: # Example usage with random data set.seed(123) n <- 100 y <- rnorm(n) Z <- matrix(rnorm(n * 5), n, 5) K <- diag(5) X <- matrix(1, n, 1) result <- mixed_solve_cpp(y, Z, K, X) ## End(Not run)
Extract Number of Subjects
n_subjects(x)n_subjects(x)
x |
A group_data object |
Integer number of unique subjects
gd <- fmrireg:::.demo_group_data_csv() n_subjects(gd)gd <- fmrireg:::.demo_group_data_csv() n_subjects(gd)
Extracts the masked time series of one or more neuroim2::NeuroVec runs
and builds an fmridataset::fmri_frame() whose feature space is the
fmridataset::volume_space() described by mask. Voxel time series are
read with neuroim2::series(), so a SparseNeuroVec does not need to be
densified.
neurovec_frame( scans, mask, TR, run_length = NULL, event_table = NULL, censor = NULL, assay = "bold" )neurovec_frame( scans, mask, TR, run_length = NULL, event_table = NULL, censor = NULL, assay = "bold" )
scans |
A |
mask |
A |
TR |
Repetition time in seconds; one value, or one per run. |
run_length |
Volumes per run; defaults to the fourth dimension of each scan. |
event_table, censor
|
As in |
assay |
Name of the assay. |
An fmri_frame over a volume_space.
Wraps fmridataset::nifti_array_source(), so headers and the mask are
read at construction but no BOLD volumes are read until the fit needs them.
For fMRIPrep derivatives see fmridataset::read_bids_bold().
nifti_frame( scans, mask, TR, run_length = NULL, event_table = NULL, censor = NULL, base_path = ".", assay = "bold" )nifti_frame( scans, mask, TR, run_length = NULL, event_table = NULL, censor = NULL, base_path = ".", assay = "bold" )
scans |
Character vector of NIfTI paths, one per run, all on the same grid. |
mask |
A NIfTI mask path, a |
TR |
Repetition time in seconds; one value, or one per run. |
run_length |
Volumes per run; defaults to the fourth header dimension of each file. |
event_table, censor
|
As in |
base_path |
Directory that relative |
assay |
Name of the assay. |
An fmri_frame over the NIfTI source's volume_space.
matrix_frame(), neurovec_frame(),
fmridataset::nifti_array_source()
Temporal-noise specification
noise_spec( struct = c("iid", "ar1", "ar2", "arp"), p = NULL, q = 0L, iter_gls = 1L, pooling = c("run", "global", "parcel"), shared_estimator = c("pooled_acvf", "mean_series"), parcels = NULL, voxelwise = FALSE, exact_first = FALSE, censor = NULL, shrink_c0 = 100L )noise_spec( struct = c("iid", "ar1", "ar2", "arp"), p = NULL, q = 0L, iter_gls = 1L, pooling = c("run", "global", "parcel"), shared_estimator = c("pooled_acvf", "mean_series"), parcels = NULL, voxelwise = FALSE, exact_first = FALSE, censor = NULL, shrink_c0 = 100L )
struct |
One of iid, ar1, ar2, or arp. |
p |
AR order for |
q |
Nonnegative moving-average order. Positive values request an
ARMA(p, q) model, with the AR order determined by |
iter_gls |
Maximum number of GLS refinement iterations. Standard
design-corrected AR fitting uses the initial OLS residuals once and holds
that estimate fixed for the GLS solve. Iterative refinement remains
available to ARMA models ( |
pooling |
Temporal covariance pooling scope. With built-in shared AR
estimation, |
shared_estimator |
Spatial estimator for a shared temporal covariance.
|
parcels |
Optional parcel labels for parcel pooling. |
voxelwise |
Estimate temporal covariance separately by voxel. In the
built-in fitter this currently requires AR-only runwise meta-estimation,
|
exact_first |
Use exact first-observation AR scaling. |
censor |
Optional censor indices, logical mask, |
shrink_c0 |
Parcel shrinkage constant used by supporting engines. |
Subjects with a non-finite outcome for a feature are dropped from that feature only (complete-case analysis per feature); a feature with no more usable subjects than predictors, or whose remaining design is rank deficient, returns NaN.
ols_t_cpp(Y, X)ols_t_cpp(Y, X)
Y |
S x P matrix (subjects x features) |
X |
S x K design matrix with intercept if desired |
List with beta (K x P), se (K x P), t (K x P), df (the residual degrees of freedom: a scalar when every feature uses all S subjects, otherwise a length-P vector), n_obs (length P, subjects used per feature), ok (P)
OLS with ONE voxelwise covariate
ols_t_vcov_cpp(Y, X, C)ols_t_vcov_cpp(Y, X, C)
Y |
S x P matrix of outcomes |
X |
S x K design matrix |
C |
S x P matrix of voxelwise covariates |
List with beta ((K+1) x P), se, t, df
Extract p-values associated with parameter estimates or test statistics from a fitted model object. This is part of a family of functions for extracting statistical measures.
p_values(x, ...) ## S3 method for class 'fmri_lm' p_values(x, type = c("betas", "contrasts"), ...)p_values(x, ...) ## S3 method for class 'fmri_lm' p_values(x, type = c("betas", "contrasts"), ...)
x |
The fitted model object |
... |
Additional arguments passed to methods |
type |
For |
A tibble or matrix containing p-values
Other statistical_measures:
coef_names(),
standard_error(),
stats()
# Create example data event_data <- data.frame( condition = factor(c("A", "B", "A", "B")), onsets = c(1, 10, 20, 30), run = c(1, 1, 1, 1) ) # Create sampling frame and dataset sframe <- sampling_frame(blocklens = 50, TR = 2) dset <- matrix_frame( matrix(rnorm(50 * 2), 50, 2), TR = 2, run_length = 50, event_table = event_data ) # Fit model fit <- fmri_lm( onsets ~ hrf(condition), block = ~run, dataset = dset ) # Extract p-values pvals <- p_values(fit)# Create example data event_data <- data.frame( condition = factor(c("A", "B", "A", "B")), onsets = c(1, 10, 20, 30), run = c(1, 1, 1, 1) ) # Create sampling frame and dataset sframe <- sampling_frame(blocklens = 50, TR = 2) dset <- matrix_frame( matrix(rnorm(50 * 2), 50, 2), TR = 2, run_length = 50, event_table = event_data ) # Fit model fit <- fmri_lm( onsets ~ hrf(condition), block = ~run, dataset = dset ) # Extract p-values pvals <- p_values(fit)
Support functions for paired differences, sign flipping, and wrapper functions for OLS and meta-analysis with voxelwise covariates. Compute Paired Within-Subject Differences
paired_diff_block(blkA, blkB, rho = 0)paired_diff_block(blkA, blkB, rho = 0)
blkA |
First group_data block |
blkB |
Second group_data block |
rho |
Optional within-subject correlation between A and B. Can be:
Default is 0 (independence). |
Creates within-subject differences (A - B) for paired t-tests from two blocks with identical subjects and features.
A new block with Y = Y_A - Y_B and propagated variance if available
Draws the design-matrix columns as time series in stacked small multiples that share one time axis. Each event regressor gets its own row, coloured by condition (designs with more than 12 event columns get one row per event term instead). Baseline regressors are drawn in grey, one row per baseline term (drift, nuisance). Run-specific columns are drawn only within their own run, and no line is joined across a run boundary. Run intercepts are omitted, since a constant carries no temporal information.
## S3 method for class 'fmri_model' plot(x, baseline = TRUE, ...)## S3 method for class 'fmri_model' plot(x, baseline = TRUE, ...)
x |
An |
baseline |
Logical; include baseline rows (default |
... |
Unused. |
A ggplot2 object.
Validates one or more fmri_jobs against the structural contract their
template implies, collecting issues per job. Intended to run on the driver
before run_jobs() / export_jobs().
preflight(x, check_files = FALSE, on_issue = c("warn", "error", "collect"))preflight(x, check_files = FALSE, on_issue = c("warn", "error", "collect"))
x |
A single fmri_job or a list of them. |
check_files |
Logical; for file-backed jobs, verify scan paths exist. |
on_issue |
One of |
Checks performed per job: template validity; every variable referenced by the
design formula and block is a column of that job's event table;
TR is positive; run lengths are consistent with the data
(matrix_frame: rows match sum(run_length);
nifti_frame: one run_length per scan file); nuisance regressor
rows match the total number of scans; and, when check_files = TRUE,
that file-backed scans exist on disk.
The design-column check uses all.vars() on the formula, so it is
deliberately conservative: a formula with a variable-valued HRF argument
(e.g. hrf(x, basis = my_basis)) may flag my_basis as a missing
column. Arguments of feature() are excluded, because those series live
in the formula environment rather than the event table. Event tables supplied
as file paths are not yet validated here.
Invisibly, an object of class fmri_preflight with an
$issues data frame (job_id, message) and $ok.
validate_template(), run_jobs(), export_jobs()
tmpl <- fmri_template(onset ~ hrf(condition), ~ run) job <- instantiate(tmpl, list(id = "sub-01", scans = matrix(rnorm(80 * 2), 80, 2), TR = 2, run_length = c(40, 40), events = data.frame(onset = c(5, 45), condition = factor(c("A", "B")), run = c(1, 2)))) preflight(job)tmpl <- fmri_template(onset ~ hrf(condition), ~ run) job <- instantiate(tmpl, list(id = "sub-01", scans = matrix(rnorm(80 * 2), 80, 2), TR = 2, run_length = c(40, 40), events = data.frame(onset = c(5, 45), condition = factor(c("A", "B")), run = c(1, 2)))) preflight(job)
Print Meta-Analysis Results
## S3 method for class 'fmri_meta' print(x, ...)## S3 method for class 'fmri_meta' print(x, ...)
x |
An fmri_meta object |
... |
Additional print arguments |
Invisibly returns the input object x
toy_meta <- structure( list( coefficients = matrix(0.2, nrow = 1), se = matrix(0.05, nrow = 1), method = "DL", robust = "none", formula = ~ condition, n_subjects = 12, n_rois = 1 ), class = "fmri_meta" ) print(toy_meta)toy_meta <- structure( list( coefficients = matrix(0.2, nrow = 1), se = matrix(0.05, nrow = 1), method = "DL", robust = "none", formula = ~ condition, n_subjects = 12, n_rois = 1 ), class = "fmri_meta" ) print(toy_meta)
Provides a colorful and informative printout.
## S3 method for class 'fmri_model' print(x, ...) ## S3 method for class 'fmri_lm' print(x, ...) ## S3 method for class 'fmri_rlm' print(x, ...)## S3 method for class 'fmri_model' print(x, ...) ## S3 method for class 'fmri_lm' print(x, ...) ## S3 method for class 'fmri_rlm' print(x, ...)
x |
An fmri_rlm object |
... |
Additional arguments passed to print.fmri_lm |
Invisibly returns the input object x
Print method for fmri_ttest_fit
## S3 method for class 'fmri_ttest_fit' print(x, ...)## S3 method for class 'fmri_ttest_fit' print(x, ...)
x |
An fmri_ttest_fit object |
... |
Additional print arguments |
Invisibly returns the input object x
Print Group Data Object
## S3 method for class 'group_data' print(x, ...)## S3 method for class 'group_data' print(x, ...)
x |
A group_data object |
... |
Additional print arguments |
Invisibly returns the input object x
Print Spatial FDR Results
## S3 method for class 'spatial_fdr_result' print(x, ...)## S3 method for class 'spatial_fdr_result' print(x, ...)
x |
A spatial_fdr_result object |
... |
Additional print arguments |
Invisibly returns the input object x
Nuisance-projection specification
projection_spec( method = c("none", "soft_subspace"), nuisance_mask = NULL, nuisance_matrix = NULL, lambda = "auto", warn_redundant = TRUE )projection_spec( method = c("none", "soft_subspace"), nuisance_mask = NULL, nuisance_matrix = NULL, lambda = "auto", warn_redundant = TRUE )
method |
Projection rule. |
nuisance_mask |
Optional logical or numeric selector for nuisance features. |
nuisance_matrix |
Optional explicit nuisance-feature matrix. |
lambda |
Nonnegative regularization strength, |
warn_redundant |
Warn when projection overlaps modeled nuisance terms. |
Compute P-values from Meta-Analysis
pvalues(object, two_tailed = TRUE)pvalues(object, two_tailed = TRUE)
object |
An fmri_meta object |
two_tailed |
Logical. Use two-tailed test (default: TRUE) |
Matrix of p-values
toy_meta <- structure( list( coefficients = matrix(c(0.2, -0.1), nrow = 1), se = matrix(c(0.05, 0.08), nrow = 1) ), class = "fmri_meta" ) pvalues(toy_meta)toy_meta <- structure( list( coefficients = matrix(c(0.2, -0.1), nrow = 1), se = matrix(c(0.05, 0.08), nrow = 1) ), class = "fmri_meta" ) pvalues(toy_meta)
Transforms correlations to Fisher's Z scale for meta-analysis.
r_to_z(r, n)r_to_z(r, n)
r |
Numeric vector or matrix of correlations |
n |
Integer scalar; sample size |
List with components:
z |
Numeric vector or matrix; Fisher's Z transformed correlations |
v |
Numeric vector or matrix; sampling variances |
r_to_z(r = 0.4, n = 30)r_to_z(r = 0.4, n = 30)
Reads complete data from all subjects' HDF5 files. Warning: This can use a lot of memory for whole-brain data.
read_h5_full(gd, stat = NULL)read_h5_full(gd, stat = NULL)
gd |
A group_data_h5 object |
stat |
Character vector of statistics to extract |
Array with dimensions (voxels, subjects, statistics)
Read jobs exported by export_jobs()
read_jobs(dir)read_jobs(dir)
dir |
Directory containing |
The list of fmri_jobs.
Reads complete data from all subjects' NIfTI files using memory mapping when possible.
read_nifti_full(gd, use_mask = NULL)read_nifti_full(gd, use_mask = NULL)
gd |
A group_data_nifti object |
use_mask |
Logical. Apply mask to data (default: TRUE if mask exists) |
List with data matrices
Reconstructs the fmri_frame from the job's dataset_spec. For
file-backed specs this is where data first becomes addressable (still lazily,
through the frame's NIfTI array source).
realize_dataset(job)realize_dataset(job)
job |
An fmri_job. |
An fmri_frame.
Extracts the fitted condition estimates as a long data frame (one row per
term x voxel), with standard error and a Wald statistic when available.
Built on coef() / standard_error() for robustness.
reduce_betas(add_id = TRUE, include_baseline = FALSE)reduce_betas(add_id = TRUE, include_baseline = FALSE)
add_id |
Logical; prepend a |
include_baseline |
Logical; include baseline/nuisance terms. Default
|
A reducer function suitable for fmri_template().
Other reducers:
reduce_contrasts(),
reduce_identity(),
reduce_write_results()
Extracts the fitted contrasts as a long data frame (one row per
contrast x voxel), optionally prefixed with the job id so results stack
cleanly across subjects. Surfaces template-level contrasts (those passed to
fmri_template(), applied by run_job()) when present, otherwise the
model's formula-embedded contrasts via coef(). Returns a zero-row frame
when no contrasts are defined.
reduce_contrasts(add_id = TRUE)reduce_contrasts(add_id = TRUE)
add_id |
Logical; prepend a |
A reducer function suitable for fmri_template().
Other reducers:
reduce_betas(),
reduce_identity(),
reduce_write_results()
The escape hatch when you want the entire fmri_lm object per subject.
Note this is the largest possible output.
reduce_identity()reduce_identity()
A reducer function suitable for fmri_template().
Other reducers:
reduce_betas(),
reduce_contrasts(),
reduce_write_results()
Calls write_results() on the fitted model, keying the output filenames by
the job's meta (subject, task, space). Returns the
vector of created file paths. This is the default reducer for large studies:
workers write compact statistical maps that group_data() / fmri_meta()
can pick up for the group level.
reduce_write_results( format = c("h5", "nifti", "gds"), stats = c("beta", "tstat"), contrasts = NULL, path = NULL, desc = "GLM", overwrite = FALSE )reduce_write_results( format = c("h5", "nifti", "gds"), stats = c("beta", "tstat"), contrasts = NULL, path = NULL, desc = "GLM", overwrite = FALSE )
format |
Output format passed to |
stats |
Which contrast statistics to write (mapped to
|
contrasts |
Optional restriction to specific contrasts. |
path |
Output directory (passed to |
desc |
BIDS |
overwrite |
Overwrite existing files. |
A reducer function suitable for fmri_template().
write_results() requires BIDS entities: each job's meta should
carry subject (defaults to the job id), task, and
space. A missing task/space surfaces as an error when
the reducer runs.
Other reducers:
reduce_betas(),
reduce_contrasts(),
reduce_identity()
Register an HRF/basis constructor
register_basis(name, constructor)register_basis(name, constructor)
name |
Character scalar used in formulas (e.g. |
constructor |
Function returning an object understood by |
Invisibly, TRUE.
fmri_lm
Register a plugin engine for fmri_lm
register_engine(name, fit, preflight = NULL, capabilities = list())register_engine(name, fit, preflight = NULL, capabilities = list())
name |
Character scalar identifier advertised to users (e.g. "friman"). |
fit |
Function invoked as |
preflight |
Optional function invoked before fitting; receives the same
arguments as |
capabilities |
Optional named list describing engine support for global
|
Invisibly, TRUE.
Adds a named backend that run_jobs() can dispatch to. This is how to plug in
a custom scheduler driver without modifying fmrireg. The builtins
"sequential" and "future" are always available.
register_run_backend(name, fn, overwrite = FALSE)register_run_backend(name, fn, overwrite = FALSE)
name |
Backend name (a string). |
fn |
A function |
overwrite |
Allow replacing an existing backend of the same name. |
The backend name, invisibly.
register_run_backend("first_only", function(jobs, run_one, ...) { list(run_one(jobs[[1]])) })register_run_backend("first_only", function(jobs, run_one, ...) { list(run_one(jobs[[1]])) })
This function is used internally by the formula processing system to lazily resolve basis names (e.g., "cca2", "cca3", "spmg1") into actual basis objects. It's part of the plugin API that allows extension packages to register custom basis functions.
resolve_basis(name, ...)resolve_basis(name, ...)
name |
Character string naming a registered basis function |
... |
Additional arguments passed to the basis constructor |
An HRF basis object
## Not run: # Resolve the cca3 basis with specific parameters basis <- resolve_basis("cca3", span = 30, TR = 2) ## End(Not run)## Not run: # Resolve the cca3 basis with specific parameters basis <- resolve_basis("cca3", span = 30, TR = 2) ## End(Not run)
Robust-estimation specification
robust_spec( type = c("none", "huber", "bisquare"), k_huber = 1.345, c_tukey = 4.685, max_iter = 2L, scale_scope = c("run", "global", "voxel"), reestimate_phi = FALSE )robust_spec( type = c("none", "huber", "bisquare"), k_huber = 1.345, c_tukey = 4.685, max_iter = 2L, scale_scope = c("run", "global", "voxel"), reestimate_phi = FALSE )
type |
Robust loss, or |
k_huber |
Tuning constant for the Huber loss. |
c_tukey |
Tuning constant for Tukey's bisquare loss. |
max_iter |
Maximum number of robust reweighting iterations. |
scale_scope |
Scope at which robust scale is estimated. |
reestimate_phi |
Re-estimate temporal-noise parameters after robust reweighting. |
List registered execution backends
run_backends()run_backends()
A character vector of backend names.
register_run_backend(), run_jobs()
Realizes the dataset, assembles the model, fits it with the template's control options (and engine, if any), and applies the template's reducer.
run_job(job, progress = FALSE)run_job(job, progress = FALSE)
job |
An fmri_job. |
progress |
Logical; show a fitting progress bar. |
The reducer's output, or the fitted fmri_lm if the template has
no reducer.
Executes a list of fmri_jobs, returning an fmri_batch_result. By
default each job is isolated: a failure is captured rather than aborting the
batch. Set parallel = TRUE to dispatch through the active
future::plan() (e.g. multisession, cluster, or a
future.batchtools cluster backend).
run_jobs( jobs, parallel = FALSE, progress = FALSE, on_error = c("isolate", "stop"), backend = NULL )run_jobs( jobs, parallel = FALSE, progress = FALSE, on_error = c("isolate", "stop"), backend = NULL )
jobs |
A single fmri_job or a list of them. |
parallel |
Logical shorthand; |
progress |
Logical; per-job fitting progress. |
on_error |
Either |
backend |
Optional execution backend name (see |
An object of class fmri_batch_result.
run_job(), batch_values(), batch_errors(), register_run_backend()
## Not run: jobs <- instantiate(tmpl, manifest) future::plan(future::multisession, workers = 4) res <- run_jobs(jobs, parallel = TRUE) values <- batch_values(res) ## End(Not run)## Not run: jobs <- instantiate(tmpl, manifest) future::plan(future::multisession, workers = 4) res <- run_jobs(jobs, parallel = TRUE) values <- batch_values(res) ## End(Not run)
Extract Standard Errors from Meta-Analysis
se(object)se(object)
object |
An fmri_meta object |
Matrix of standard errors
toy_meta <- structure( list(se = matrix(c(0.05, 0.06), nrow = 1)), class = "fmri_meta" ) se(toy_meta)toy_meta <- structure( list(se = matrix(c(0.05, 0.06), nrow = 1)), class = "fmri_meta" ) se(toy_meta)
This function simulates an fMRI time series for multiple experimental conditions with specified parameters. It generates a controlled event-related teaching design with randomized inter-stimulus intervals and condition orders.
simulate_bold_signal( ncond, hrf = fmrihrf::HRF_SPMG1, nreps = 12, amps = rep(1, ncond), isi = c(3, 6), ampsd = 0, TR = 1.5 )simulate_bold_signal( ncond, hrf = fmrihrf::HRF_SPMG1, nreps = 12, amps = rep(1, ncond), isi = c(3, 6), ampsd = 0, TR = 1.5 )
ncond |
The number of conditions to simulate. |
hrf |
The hemodynamic response function to use (default is fmrihrf::HRF_SPMG1). |
nreps |
The number of repetitions per condition (default is 12). |
amps |
A vector of amplitudes for each condition (default is a vector of 1s with length ncond). |
isi |
A vector of length 2 specifying the range of inter-stimulus intervals to sample from (default is c(3, 6) seconds). |
ampsd |
The standard deviation of the amplitudes (default is 0). |
TR |
The repetition time of the fMRI acquisition (default is 1.5 seconds). |
A list with the following components:
onset: A vector of the onset times for each trial
condition: A vector of condition labels for each trial
mat: A matrix containing the simulated fMRI time series:
Column 1: Time points (in seconds)
Columns 2:(ncond+1): Simulated BOLD responses for each condition
# Simulate 3 conditions with different amplitudes sim <- simulate_bold_signal(ncond = 3, amps = c(1, 1.5, 2), TR = 2) # Plot the simulated time series matplot(sim$mat[,1], sim$mat[,-1], type = "l", xlab = "Time (s)", ylab = "BOLD Response")# Simulate 3 conditions with different amplitudes sim <- simulate_bold_signal(ncond = 3, amps = c(1, 1.5, 2), TR = 2) # Plot the simulated time series matplot(sim$mat[,1], sim$mat[,-1], type = "l", xlab = "Time (s)", ylab = "BOLD Response")
Generates time-series (columns) with a single set of onsets, but
resampled amplitudes/durations for each column if amplitude_sd>0
or duration_sd>0. Each column also gets independent noise. The result
is a list containing:
time_series: a fmri_frame with .
The event_table uses the first column's amplitude/duration draws.
ampmat: an matrix of per-column amplitudes.
durmat: an matrix of per-column durations.
hrf_info: info about the HRF.
noise_params: info about noise generation (type + AR coefficients + SD).
simulate_fmri_matrix( n = 1, total_time = 240, TR = 2, hrf = fmrihrf::HRF_SPMG1, n_events = 10, onsets = NULL, isi_dist = c("even", "uniform", "exponential"), isi_min = 2, isi_max = 6, isi_rate = 0.25, durations = 0, duration_sd = 0, duration_dist = c("lognormal", "gamma"), amplitudes = 1, amplitude_sd = 0, amplitude_dist = c("lognormal", "gamma", "gaussian"), single_trial = FALSE, noise_type = c("none", "white", "ar1", "ar2"), noise_ar = NULL, noise_sd = 1, random_seed = NULL, verbose = FALSE, buffer = 16 )simulate_fmri_matrix( n = 1, total_time = 240, TR = 2, hrf = fmrihrf::HRF_SPMG1, n_events = 10, onsets = NULL, isi_dist = c("even", "uniform", "exponential"), isi_min = 2, isi_max = 6, isi_rate = 0.25, durations = 0, duration_sd = 0, duration_dist = c("lognormal", "gamma"), amplitudes = 1, amplitude_sd = 0, amplitude_dist = c("lognormal", "gamma", "gaussian"), single_trial = FALSE, noise_type = c("none", "white", "ar1", "ar2"), noise_ar = NULL, noise_sd = 1, random_seed = NULL, verbose = FALSE, buffer = 16 )
n |
Number of time-series (columns). |
total_time |
Numeric. Total scan length (seconds). |
TR |
Numeric. Repetition time (seconds). |
hrf |
Hemodynamic response function, e.g. |
n_events |
Number of events (ignored if |
onsets |
Optional non-empty numeric vector of finite event onsets in
|
isi_dist |
One of |
isi_min, isi_max
|
For |
isi_rate |
For |
durations |
Numeric, scalar or length- |
duration_sd |
Numeric. If >0, random variation in durations. |
duration_dist |
|
amplitudes |
Numeric, scalar or length- |
amplitude_sd |
Numeric. If >0, random variation in amplitudes. |
amplitude_dist |
|
single_trial |
If TRUE, each event is a separate single-trial regressor that gets summed. |
noise_type |
|
noise_ar |
Numeric vector for AR(1) or AR(2). If missing or insufficient, defaults are used (0.3 for AR(1); c(0.3,0.2) for AR(2)). |
noise_sd |
Std dev of the noise. |
random_seed |
Optional integer for reproducibility. |
verbose |
If TRUE, prints messages. |
buffer |
Numeric seconds appended to the end of the time grid to avoid edge truncation (default: 16). |
If noise_type="ar1" and you do not provide noise_ar, we
default to c(0.3).
If noise_type="ar2" and you do not provide a 2-element noise_ar,
we default to c(0.3, 0.2).
Onsets are either provided or generated once for all columns.
Amplitudes/durations are re-sampled inside the loop so each
column can differ randomly. The final arrays ampmat and durmat
each have one column per time-series.
The fmri_frame's event_table records the first column's
amplitudes/durations. If you need each column's, see ampmat and
durmat.
A list containing:
time_seriesA fmri_frame with data
and event_table for the first column's random draws.
ampmatAn numeric matrix of amplitudes.
durmatAn numeric matrix of durations.
hrf_infoA list with HRF metadata.
noise_paramsA list describing noise generation.
This function simulates structured fMRI-like teaching noise by combining:
Temporal autocorrelation using an ARMA model
Low-frequency drift
Physiological noise (cardiac and respiratory)
simulate_noise_vector( n, TR = 1.5, ar = c(0.3), ma = c(0.5), sd = 1, drift_freq = 1/128, drift_amplitude = 2, physio = TRUE, seed = NULL )simulate_noise_vector( n, TR = 1.5, ar = c(0.3), ma = c(0.5), sd = 1, drift_freq = 1/128, drift_amplitude = 2, physio = TRUE, seed = NULL )
n |
The number of time points in the fMRI time series |
TR |
The repetition time in seconds (default is 1.5) |
ar |
A numeric vector containing autoregressive (AR) coefficients (default is c(0.3)) |
ma |
A numeric vector containing moving average (MA) coefficients (default is c(0.5)) |
sd |
The standard deviation of the white noise component (default is 1) |
drift_freq |
Frequency of the low-frequency drift in Hz (default is 1/128) |
drift_amplitude |
Amplitude of the low-frequency drift (default is 2) |
physio |
Logical; whether to add simulated physiological noise (default is TRUE) |
seed |
An optional seed for reproducibility (default is NULL) |
A numeric vector containing the simulated fMRI noise
# Simulate noise for a 5-minute scan with TR=2s n_timepoints <- 150 # 5 minutes * 60 seconds / 2s TR noise <- simulate_noise_vector(n_timepoints, TR = 2) plot(noise, type = "l", xlab = "Time Point", ylab = "Signal")# Simulate noise for a 5-minute scan with TR=2s n_timepoints <- 150 # 5 minutes * 60 seconds / 2s TR noise <- simulate_noise_vector(n_timepoints, TR = 2) plot(noise, type = "l", xlab = "Time Point", ylab = "Signal")
This function simulates a complete fMRI dataset by combining task-related signals with a declared structured-noise model. It returns both the clean signals and the noisy data.
simulate_simple_dataset( ncond, nreps = 12, TR = 1.5, snr = 0.5, hrf = fmrihrf::HRF_SPMG1, seed = NULL )simulate_simple_dataset( ncond, nreps = 12, TR = 1.5, snr = 0.5, hrf = fmrihrf::HRF_SPMG1, seed = NULL )
ncond |
Number of conditions to simulate |
nreps |
Number of repetitions per condition (default is 12) |
TR |
Repetition time in seconds (default is 1.5) |
snr |
Finite positive signal-to-noise ratio, defined as the empirical standard deviation of the clean condition signals divided by the empirical standard deviation of the complete structured-noise matrix (default is 0.5) |
hrf |
Hemodynamic response function to use (default is fmrihrf::HRF_SPMG1) |
seed |
Optional seed for reproducibility (default is NULL). Reusing a
seed with a different |
A list containing:
clean: The simulated signals without noise (from simulate_bold_signal)
noisy: The signals with added noise
noise: The simulated noise component
onsets: Trial onset times
conditions: Condition labels for each trial
# Simulate a dataset with 3 conditions data <- simulate_simple_dataset(ncond = 3, TR = 2, snr = 0.5) # Plot clean and noisy data old_par <- par(no.readonly = TRUE) par(mfrow = c(2,1)) matplot(data$clean$mat[,1], data$clean$mat[,-1], type = "l", main = "Clean Signal", xlab = "Time (s)", ylab = "BOLD") matplot(data$noisy[,1], data$noisy[,-1], type = "l", main = "Noisy Signal", xlab = "Time (s)", ylab = "BOLD") par(old_par)# Simulate a dataset with 3 conditions data <- simulate_simple_dataset(ncond = 3, TR = 2, snr = 0.5) # Plot clean and noisy data old_par <- par(no.readonly = TRUE) par(mfrow = c(2,1)) matplot(data$clean$mat[,1], data$clean$mat[,-1], type = "l", main = "Clean Signal", xlab = "Time (s)", ylab = "BOLD") matplot(data$noisy[,1], data$noisy[,-1], type = "l", main = "Noisy Signal", xlab = "Time (s)", ylab = "BOLD") par(old_par)
Computes the soft (ridge-regularized) projection matrix that removes variance explainable by the nuisance subspace while avoiding overfitting.
soft_projection(N, lambda = "auto", Y = NULL)soft_projection(N, lambda = "auto", Y = NULL)
N |
Nuisance matrix (time x nuisance_voxels). Typically extracted from white matter and CSF voxels. Can have many columns (thousands); computational cost depends on min(nrow, ncol) due to SVD. |
lambda |
Ridge penalty controlling shrinkage strength:
|
Y |
Optional data matrix for GCV-based lambda selection. Required if
|
A list with class "soft_projection" containing:
P_lambda |
Function that applies the projection to a matrix |
lambda |
The selected/specified lambda value |
method |
Method used: "singular_value_heuristic", "gcv", or "user_specified" |
effective_df |
Effective degrees of freedom removed (sum of shrinkage factors) |
n_nuisance |
Number of nuisance columns in N |
n_timepoints |
Number of timepoints |
# Create nuisance matrix (e.g., from WM/CSF voxels) set.seed(123) N <- matrix(rnorm(100 * 20), nrow = 100, ncol = 20) proj <- soft_projection(N, lambda = "auto") # Apply to data and design Y <- matrix(rnorm(100 * 50), nrow = 100, ncol = 50) Y_clean <- proj$P_lambda(Y) # Full workflow: project both data and design X <- cbind(1, rnorm(100), rnorm(100)) # intercept + 2 predictors cleaned <- apply_soft_projection(proj, Y, X) # Now fit GLM with cleaned$Y and cleaned$X # Using GCV for lambda selection (data-driven) proj_gcv <- soft_projection(N, lambda = "gcv", Y = Y) print(proj_gcv) # Shows selected lambda and effective df# Create nuisance matrix (e.g., from WM/CSF voxels) set.seed(123) N <- matrix(rnorm(100 * 20), nrow = 100, ncol = 20) proj <- soft_projection(N, lambda = "auto") # Apply to data and design Y <- matrix(rnorm(100 * 50), nrow = 100, ncol = 50) Y_clean <- proj$P_lambda(Y) # Full workflow: project both data and design X <- cbind(1, rnorm(100), rnorm(100)) # intercept + 2 predictors cleaned <- apply_soft_projection(proj, Y, X) # Now fit GLM with cleaned$Y and cleaned$X # Using GCV for lambda selection (data-driven) proj_gcv <- soft_projection(N, lambda = "gcv", Y = Y) print(proj_gcv) # Shows selected lambda and effective df
Creates a configuration object for soft subspace projection.
soft_subspace_options( enabled = FALSE, nuisance_mask = NULL, nuisance_matrix = NULL, lambda = "auto", warn_redundant = TRUE )soft_subspace_options( enabled = FALSE, nuisance_mask = NULL, nuisance_matrix = NULL, lambda = "auto", warn_redundant = TRUE )
enabled |
Logical. Whether to apply soft subspace projection. |
nuisance_mask |
Path to NIfTI mask or logical vector indicating nuisance voxels. |
nuisance_matrix |
Pre-computed nuisance timeseries matrix (alternative to mask). |
lambda |
Ridge penalty: numeric, "auto", or "gcv". |
warn_redundant |
Logical. Warn if baseline model contains nuisance terms. |
A list of class "soft_subspace_options".
# Using a mask file opts <- soft_subspace_options( enabled = TRUE, nuisance_mask = "path/to/wm_csf_mask.nii.gz", lambda = "auto" ) # Using pre-computed nuisance matrix N <- matrix(rnorm(100 * 20), 100, 20) opts <- soft_subspace_options( enabled = TRUE, nuisance_matrix = N, lambda = 0.5 )# Using a mask file opts <- soft_subspace_options( enabled = TRUE, nuisance_mask = "path/to/wm_csf_mask.nii.gz", lambda = "auto" ) # Using pre-computed nuisance matrix N <- matrix(rnorm(100 * 20), 100, 20) opts <- soft_subspace_options( enabled = TRUE, nuisance_matrix = N, lambda = 0.5 )
Ridge-regularized projection to remove nuisance variance without explicitly choosing components or adding confound regressors. This is "CompCor without choosing components."
Traditional CompCor extracts K principal components from white matter/CSF and regresses them out. This requires choosing K, which is arbitrary.
Soft subspace projection instead treats the entire WM/CSF timeseries as a nuisance basis and removes it with shrinkage. Each nuisance direction is partially removed proportional to its variance, avoiding the hard keep/delete decision of PCA truncation.
Given nuisance matrix N (time x nuisance_voxels), the soft projection is:
Applied to data Y and design X:
Y_clean = P_lambda %*% Y
X_clean = P_lambda %*% X (important to avoid bias)
The shrinkage parameter lambda controls aggressiveness:
Small lambda: aggressive removal (risk removing signal)
Large lambda: gentle removal (risk leaving nuisance)
Three selection methods are available:
Singular value heuristic: lambda = median(d^2) where d
are singular values of N. Fast, stable, no tuning. Components with variance
below median are heavily shrunk; above median lightly shrunk.
Generalized cross-validation minimizes RSS/(1-df/n)^2
for ridge regression of Y on N. Finds lambda giving best leave-one-out
prediction without actually doing LOO. Requires Y.
Computational cost is O(k) per evaluation where k = min(n, p).
User-specified lambda value.
The soft subspace projection workflow has three steps:
Extract nuisance timeseries from WM/CSF mask (or provide pre-computed)
Create the soft projection operator
Apply to both data Y and design matrix X before GLM fitting
For use within fmri_lm(), see the convenience parameter
nuisance_projection or the more detailed soft_subspace_options.
Soft subspace projection is most beneficial when:
You have physiological noise from WM/CSF that isn't captured by motion parameters
Traditional CompCor requires arbitrary component selection
You want to avoid adding many confound regressors to the design
Consider alternatives when:
Motion parameters alone sufficiently control artifacts
You prefer explicit confound regressors for interpretability
Data has very few timepoints relative to nuisance dimensions
Performs spatially-aware FDR control using structure-adaptive weighted Benjamini-Hochberg (SABHA-style) procedure. This method leverages spatial structure in the data to increase power while controlling the false discovery rate.
spatial_fdr( z = NULL, p = NULL, group = NULL, alpha = 0.05, tau = 0.5, lambda = 1, neighbors = NULL, min_pi0 = 0.05, empirical_null = TRUE, verbose = FALSE )spatial_fdr( z = NULL, p = NULL, group = NULL, alpha = 0.05, tau = 0.5, lambda = 1, neighbors = NULL, min_pi0 = 0.05, empirical_null = TRUE, verbose = FALSE )
z |
Numeric vector of Z-values (one per feature). Provide either z or p, not both. |
p |
Numeric vector of p-values (two-sided). Provide either z or p, not both. |
group |
Integer, character or factor vector of group IDs for each feature
(e.g., parcel IDs, block IDs). Must have same length as z or p. Groups are
numbered 1..G by sorted unique value, or by factor level (unused levels
included) when |
alpha |
Numeric scalar; FDR level to control (default: 0.05) |
tau |
Numeric scalar; Storey threshold in (0,1) for |
lambda |
Numeric scalar; smoothing strength for groupwise |
neighbors |
Optional list of length G (number of groups) where element g
is an integer vector of 1-based neighbor IDs of group g, both in the group
numbering described under |
min_pi0 |
Numeric scalar; lower bound for |
empirical_null |
Logical; if TRUE, estimate null distribution parameters
( |
verbose |
Logical; print progress messages (default: FALSE) |
This function implements a spatially-aware multiple testing procedure that:
Estimates the proportion of null hypotheses (pi0) within spatial groups
Optionally smooths these estimates across neighboring groups
Uses the pi0 estimates to weight the Benjamini-Hochberg procedure
Provides more power in regions with true signal while maintaining FDR control
Weights are of each feature's group, normalised to sum to the
number of finite p-values, and the weighted BH procedure of Genovese,
Roeder & Wasserman (2006) is applied: reject .
Empirical null. When z is supplied with
empirical_null = TRUE, and are estimated by
Efron's (2004) central matching: a quadratic is fitted by Poisson regression
to the log histogram counts of z within robust SDs of the
median. Because a minority of non-null features only changes the intercept
of the log density near the centre and not its curvature, this estimate is
consistent under and stays close to it when, for example, 10\
of features are true signals. With fewer than 500 finite z-values (or if the
fit fails) the median and normal-consistent MAD of all finite z-values are
used instead; with fewer than 10 the theoretical null is used.
The method is particularly effective for:
Voxelwise analyses with spatial clustering of signal
Parcel-based analyses with anatomical or functional grouping
Any scenario where hypotheses can be grouped spatially or functionally
Object of class "spatial_fdr_result" containing:
reject |
Logical vector indicating rejected hypotheses (discoveries) |
q |
Numeric vector of FDR-adjusted p-values (q-values) |
p |
Numeric vector of two-sided p-values used for testing |
weights |
Numeric vector of normalized weights used in weighted BH |
pi0_raw |
Numeric vector of raw |
pi0_smooth |
Numeric vector of smoothed |
threshold |
Numeric scalar; BH threshold used for rejection |
k |
Integer scalar; number of rejections |
mu0 |
Numeric scalar; estimated null mean (if empirical_null = TRUE) |
sigma0 |
Numeric scalar; estimated null SD (if empirical_null = TRUE) |
groups |
Integer vector of group IDs (1..G) in the numbering described
under |
G |
Integer scalar; number of groups |
m_g |
Integer vector; number of features per group |
alpha |
Numeric scalar; FDR level used |
coef_name |
Character scalar; name of coefficient (for S3 method) |
Benjamini & Hochberg (1995). Controlling the false discovery rate.
Storey (2002). A direct approach to false discovery rates.
Hu et al. (2010). False discovery rate control with groups (SABHA).
Genovese, Roeder & Wasserman (2006). False discovery control with p-value weighting.
Efron (2004). Large-scale simultaneous hypothesis testing: the choice of a null hypothesis.
# Simple synthetic example set.seed(123) n <- 1000 z_scores <- c(rnorm(800), rnorm(200, mean = 2)) # 200 true signals group_ids <- rep(1:10, each = 100) # 10 groups of 100 features each # Basic usage without spatial smoothing result <- spatial_fdr(z = z_scores, group = group_ids, alpha = 0.05) summary(result) # Create simple neighbor structure (each group neighbors with adjacent groups) neighbors <- lapply(1:10, function(i) { c(if(i > 1) i-1, if(i < 10) i+1) }) # With spatial smoothing result_smooth <- spatial_fdr(z = z_scores, group = group_ids, neighbors = neighbors, lambda = 1.0) print(result_smooth)# Simple synthetic example set.seed(123) n <- 1000 z_scores <- c(rnorm(800), rnorm(200, mean = 2)) # 200 true signals group_ids <- rep(1:10, each = 100) # 10 groups of 100 features each # Basic usage without spatial smoothing result <- spatial_fdr(z = z_scores, group = group_ids, alpha = 0.05) summary(result) # Create simple neighbor structure (each group neighbors with adjacent groups) neighbors <- lapply(1:10, function(i) { c(if(i > 1) i-1, if(i < 10) i+1) }) # With spatial smoothing result_smooth <- spatial_fdr(z = z_scores, group = group_ids, neighbors = neighbors, lambda = 1.0) print(result_smooth)
Extract standard errors of parameter estimates from a fitted model object. This is part of a family of functions for extracting statistical measures.
standard_error(x, ...) ## S3 method for class 'fmri_latent_lm' standard_error(x, type = c("betas", "contrasts"), recon = FALSE, ...) ## S3 method for class 'fmri_lm' standard_error(x, type = c("betas", "contrasts"), ...)standard_error(x, ...) ## S3 method for class 'fmri_latent_lm' standard_error(x, type = c("betas", "contrasts"), recon = FALSE, ...) ## S3 method for class 'fmri_lm' standard_error(x, type = c("betas", "contrasts"), ...)
x |
The fitted model object |
... |
Additional arguments passed to methods |
type |
For |
recon |
Logical; whether to reconstruct the full matrix representation (default: FALSE) |
A tibble or matrix containing standard errors of parameter
estimates. For fmri_lm fits, one row per voxel and one column per
term, in the same orientation as coef().
Other statistical_measures:
coef_names(),
p_values(),
stats()
# Create example data event_data <- data.frame( condition = factor(c("A", "B", "A", "B")), onsets = c(1, 10, 20, 30), run = c(1, 1, 1, 1) ) # Create sampling frame and dataset sframe <- sampling_frame(blocklens = 50, TR = 2) dset <- matrix_frame( matrix(rnorm(50 * 2), 50, 2), TR = 2, run_length = 50, event_table = event_data ) # Fit model fit <- fmri_lm( onsets ~ hrf(condition), block = ~run, dataset = dset ) # Extract standard errors se <- standard_error(fit)# Create example data event_data <- data.frame( condition = factor(c("A", "B", "A", "B")), onsets = c(1, 10, 20, 30), run = c(1, 1, 1, 1) ) # Create sampling frame and dataset sframe <- sampling_frame(blocklens = 50, TR = 2) dset <- matrix_frame( matrix(rnorm(50 * 2), 50, 2), TR = 2, run_length = 50, event_table = event_data ) # Fit model fit <- fmri_lm( onsets ~ hrf(condition), block = ~run, dataset = dset ) # Extract standard errors se <- standard_error(fit)
Extract test statistics (e.g., t-statistics, F-statistics) from a fitted model object. This is part of a family of functions for extracting statistical measures.
stats(x, ...) ## S3 method for class 'fmri_lm' stats(x, type = c("betas", "contrasts", "F"), ...)stats(x, ...) ## S3 method for class 'fmri_lm' stats(x, type = c("betas", "contrasts", "F"), ...)
x |
The fitted model object |
... |
Additional arguments passed to methods |
type |
For |
A tibble or matrix containing test statistics. For fmri_lm
fits, a tibble with one row per voxel and one column per term, in the
same orientation as coef().
coef() for the estimates.
Other statistical_measures:
coef_names(),
p_values(),
standard_error()
# Create example data event_data <- data.frame( condition = factor(c("A", "B", "A", "B")), onsets = c(1, 10, 20, 30), run = c(1, 1, 1, 1) ) # Create sampling frame and dataset sframe <- sampling_frame(blocklens = 50, TR = 2) dset <- matrix_frame( matrix(rnorm(50 * 2), 50, 2), TR = 2, run_length = 50, event_table = event_data ) # Fit model fit <- fmri_lm( onsets ~ hrf(condition), block = ~run, dataset = dset ) # Extract test statistics tstats <- stats(fit)# Create example data event_data <- data.frame( condition = factor(c("A", "B", "A", "B")), onsets = c(1, 10, 20, 30), run = c(1, 1, 1, 1) ) # Create sampling frame and dataset sframe <- sampling_frame(blocklens = 50, TR = 2) dset <- matrix_frame( matrix(rnorm(50 * 2), 50, 2), TR = 2, run_length = 50, event_table = event_data ) # Fit model fit <- fmri_lm( onsets ~ hrf(condition), block = ~run, dataset = dset ) # Extract test statistics tstats <- stats(fit)
Summary of Meta-Analysis Results
## S3 method for class 'fmri_meta' summary(object, threshold = 0.05, ...)## S3 method for class 'fmri_meta' summary(object, threshold = 0.05, ...)
object |
An fmri_meta object |
threshold |
P-value threshold for significance (default: 0.05) |
... |
Additional summary arguments |
A list containing summary statistics invisibly
toy_meta <- structure( list( coefficients = matrix(c(0.3, -0.1), nrow = 1, dimnames = list(NULL, c("A", "B"))), se = matrix(c(0.05, 0.07), nrow = 1), method = "DL", robust = "none", formula = ~ condition, n_subjects = 12, n_rois = 1 ), class = "fmri_meta" ) summary(toy_meta, threshold = 0.1)toy_meta <- structure( list( coefficients = matrix(c(0.3, -0.1), nrow = 1, dimnames = list(NULL, c("A", "B"))), se = matrix(c(0.05, 0.07), nrow = 1), method = "DL", robust = "none", formula = ~ condition, n_subjects = 12, n_rois = 1 ), class = "fmri_meta" ) summary(toy_meta, threshold = 0.1)
Summary method for fmri_ttest_fit
## S3 method for class 'fmri_ttest_fit' summary(object, ...)## S3 method for class 'fmri_ttest_fit' summary(object, ...)
object |
An fmri_ttest_fit object |
... |
Additional summary arguments |
Invisibly returns the input object
Summary of Group Data Object
## S3 method for class 'group_data' summary(object, ...)## S3 method for class 'group_data' summary(object, ...)
object |
A group_data object |
... |
Additional summary arguments |
Invisibly returns the input object
Summary of Spatial FDR Results
## S3 method for class 'spatial_fdr_result' summary(object, ...)## S3 method for class 'spatial_fdr_result' summary(object, ...)
object |
A spatial_fdr_result object |
... |
Additional summary arguments |
Invisibly returns the input object
Helper function to convert t-statistics and df to beta and SE estimates
t_to_beta_se(t, df, n = NULL)t_to_beta_se(t, df, n = NULL)
t |
T-statistic values |
df |
Degrees of freedom |
n |
Sample size (optional, improves SE estimation) |
List with beta and se estimates
Converts t-statistics and degrees of freedom to standardized mean differences (Cohen's d) and their sampling variances for meta-analysis.
t_to_d(t, df, n = NULL)t_to_d(t, df, n = NULL)
t |
Numeric vector or matrix of t-statistics |
df |
Numeric scalar or vector; degrees of freedom (matching t) |
n |
Numeric scalar; sample size per group (for two-sample t-tests) |
List with components:
d |
Numeric vector or matrix; standardized mean differences |
v |
Numeric vector or matrix; sampling variances |
t_to_d(t = 2, df = 18) t_to_d(t = 2, df = 18, n = 20)t_to_d(t = 2, df = 18) t_to_d(t = 2, df = 18, n = 20)
Extract design matrices for individual terms from an fmri_model object.
## S3 method for class 'fmri_model' term_matrices(x, blocknum = NULL, ...)## S3 method for class 'fmri_model' term_matrices(x, blocknum = NULL, ...)
x |
An fmri_model object |
blocknum |
Optional vector of block numbers to extract matrices for |
... |
Additional arguments (currently unused) |
A list of matrices, one for each term in the model
fm <- fmrireg:::.demo_fmri_model() term_matrices(fm)fm <- fmrireg:::.demo_fmri_model() term_matrices(fm)
A restrained theme used by all fmrireg plots: hairline major grid only, left-aligned titles, and muted axis ink so the data carry the contrast. Add it to any ggplot, or override it by adding your own theme afterwards.
theme_fmrireg( base_size = 11, base_family = "", grid = c("y", "x", "xy", "none") )theme_fmrireg( base_size = 11, base_family = "", grid = c("y", "x", "xy", "none") )
base_size |
Base font size in points. |
base_family |
Base font family. |
grid |
Which major gridlines to keep: |
fmrireg plots use this theme by default. To use a different base theme for
every fmrireg plot, set options(fmrireg.plot_theme = <theme>), for example
options(fmrireg.plot_theme = ggplot2::theme_bw()).
A ggplot2::theme() object.
library(ggplot2) ggplot(mtcars, aes(wt, mpg)) + geom_point() + theme_fmrireg()library(ggplot2) ggplot(mtcars, aes(wt, mpg)) + geom_point() + theme_fmrireg()
Minimal tidy generic to support fmri_meta tidy() without requiring broom.
tidy(x, ...)tidy(x, ...)
x |
object |
... |
passed to methods |
A tidy data frame with model coefficients and statistics
Convert fitted_hrf() output into a long tibble suitable for plotting.
tidy_fitted_hrf( x, sample_at = seq(0, 24, by = 1), term = NULL, term_match = c("contains", "exact", "regex"), voxel = 1L, average_voxels = FALSE, ... )tidy_fitted_hrf( x, sample_at = seq(0, 24, by = 1), term = NULL, term_match = c("contains", "exact", "regex"), voxel = 1L, average_voxels = FALSE, ... )
x |
An |
sample_at |
Numeric vector of time points passed to |
term |
Optional term selector. When |
term_match |
Matching mode for |
voxel |
Integer voxel index to extract from each term's prediction matrix. |
average_voxels |
Logical; if |
... |
Additional arguments passed to |
A tibble with columns term, time, condition, estimate, and
voxel.
Tidy a smooth HRF estimate
## S3 method for class 'fmri_hrf_estimate' tidy(x, curve = NULL, voxel = NULL, ...)## S3 method for class 'fmri_hrf_estimate' tidy(x, curve = NULL, voxel = NULL, ...)
x |
An |
curve |
Optional curve names to retain. |
voxel |
Optional voxel names or indices to retain. |
... |
Unused. |
A tibble with one row per time, curve, and voxel.
Tidy Meta-Analysis Results
## S3 method for class 'fmri_meta' tidy(x, conf.int = FALSE, conf.level = 0.95, ...)## S3 method for class 'fmri_meta' tidy(x, conf.int = FALSE, conf.level = 0.95, ...)
x |
An fmri_meta object |
conf.int |
Logical. Include confidence intervals (default: FALSE) |
conf.level |
Confidence level (default: 0.95) |
... |
Additional arguments |
A tibble with tidy results
Validate an fmri_job
validate_job(x)validate_job(x)
x |
An fmri_job. |
TRUE invisibly if valid; otherwise an error is raised.
Re-checks the structural invariants of a template. Useful in preflight before fanning a model out over many subjects.
validate_template(x)validate_template(x)
x |
An fmri_template. |
TRUE invisibly if valid; otherwise an error is raised.
Extract the structured variance model
variance_model(x)variance_model(x)
x |
A fitted |
An fmri_lm_variance_model object.
Variance and reference-distribution specification
variance_spec( method = c("model", "hac", "sandwich"), max_lag = "auto", taper = c("none", "tukey", "parzen"), debias = TRUE, df = c("residual", "satterthwaite") )variance_spec( method = c("model", "hac", "sandwich"), max_lag = "auto", taper = c("none", "tukey", "parzen"), debias = TRUE, df = c("residual", "satterthwaite") )
method |
Covariance estimator: model-based, HAC, or sandwich. |
max_lag |
Maximum HAC lag, or |
taper |
HAC lag-window taper. |
debias |
Apply the estimator's finite-sample correction. |
df |
Reference degrees-of-freedom method. |
Functions for computing volume-level quality metrics (DVARS) and converting them to weights for weighted least squares fitting. This implements "soft scrubbing" - downweighting bad volumes rather than hard censoring them.
Traditional fMRI artifact removal uses hard censoring ("scrubbing"), where volumes exceeding a threshold are completely removed. This loses temporal information and can create discontinuities.
Soft scrubbing instead assigns each volume a weight between 0 and 1 based on its quality. High-quality volumes receive full weight; artifacts receive reduced weight. This preserves temporal continuity while still downweighting problematic data.
DVARS (Derivative of VARiance across voxels) measures the root mean square of the temporal derivative across all voxels:
Interpretation:
Low DVARS: signal changed smoothly from previous volume (good)
High DVARS: signal changed rapidly (possible artifact or motion)
Normalized DVARS: values ~1 are typical; >1.5 suggests artifacts
Three methods convert DVARS to weights, offering different trade-offs:
Formula: w = 1 / (1 + dvars^2)
Properties: Smooth, continuous decay. A volume with DVARS=1 gets weight 0.5. Most conservative choice - provides gentle downweighting even for moderately elevated DVARS. Use when: You want smooth, gradual artifact handling without sharp transitions.
Formula: Sigmoid decay above threshold
Properties: Volumes below threshold get full weight; above threshold, weights decay smoothly. Steepness parameter controls how rapidly weights drop. Use when: You have a clear idea of what "acceptable" DVARS looks like (e.g., 1.5x median) and want to preserve good volumes fully.
Formula: Tukey bisquare (1 - u^2)^2 for |u| <= 1
Properties: Complete downweighting for extreme values. Volumes beyond 2x threshold get zero weight. Most aggressive choice. Use when: You have clear artifacts that should be fully excluded, similar to hard scrubbing but with smooth transitions.
Both approaches handle artifacts, but through different mechanisms:
Volume weighting downweights entire timepoints uniformly across all voxels. Best when artifacts affect the whole brain (head motion, scanner spikes).
Robust fitting (Huber/Tukey bisquare) downweights outlier residuals voxel-by-voxel. Best when artifacts are spatially localized or when temporal structure matters.
Combined approach: Use volume weighting for global artifacts
plus robust fitting for residual voxel-level outliers. Set
robust = TRUE, volume_weights = TRUE in fmri_lm().
For use within fmri_lm(), see the convenience parameter
volume_weights or the more detailed volume_weights_options.
Convenience function that computes DVARS and converts to weights in one step. This is the main user-facing function for volume quality weighting.
volume_weights( Y, method = "inverse_squared", threshold = 1.5, return_dvars = FALSE )volume_weights( Y, method = "inverse_squared", threshold = 1.5, return_dvars = FALSE )
Y |
Numeric matrix of fMRI data (time x voxels). |
method |
Weighting method passed to |
threshold |
Threshold passed to |
return_dvars |
Logical. If TRUE, return both weights and DVARS. |
If return_dvars is FALSE, a numeric vector of weights. If TRUE, a list with components "weights" and "dvars".
set.seed(123) Y <- matrix(rnorm(100 * 50), nrow = 100, ncol = 50) Y[50, ] <- Y[50, ] + 5 # Add artifact # One-step computation of weights result <- volume_weights(Y, return_dvars = TRUE) cat("DVARS at artifact:", round(result$dvars[50], 2), "\n") cat("Weight at artifact:", round(result$weights[50], 3), "\n") # With Tukey method for more aggressive downweighting w_tukey <- volume_weights(Y, method = "tukey")set.seed(123) Y <- matrix(rnorm(100 * 50), nrow = 100, ncol = 50) Y[50, ] <- Y[50, ] + 5 # Add artifact # One-step computation of weights result <- volume_weights(Y, return_dvars = TRUE) cat("DVARS at artifact:", round(result$dvars[50], 2), "\n") cat("Weight at artifact:", round(result$weights[50], 3), "\n") # With Tukey method for more aggressive downweighting w_tukey <- volume_weights(Y, method = "tukey")
Volume-weight specification
weights_spec( method = c("none", "inverse_squared", "soft_threshold", "tukey"), threshold = 1.5, values = NULL )weights_spec( method = c("none", "inverse_squared", "soft_threshold", "tukey"), threshold = 1.5, values = NULL )
method |
Volume-weight construction rule. |
threshold |
Positive threshold or tuning constant for the selected rule. |
values |
Optional explicit finite, nonnegative volume weights. Supplying
|
Welch two-sample t-test across features
welch_t_cpp(Y, g_in)welch_t_cpp(Y, g_in)
Y |
S x P matrix |
g_in |
Length S vector of group indicators (1/2 or 0/1) |
List with muA, muB, t, df (Welch), nA, nB
Generic function to export statistical maps and analysis results from fitted fMRI models to standardized file formats with appropriate metadata.
write_results(x, ...)write_results(x, ...)
x |
A fitted fMRI model object |
... |
Additional arguments passed to methods |
Invisible list of created file paths
Exports statistical maps from an fmri_lm object with BIDS-compliant naming and JSON metadata sidecars. Supported image outputs are HDF5 LabeledVolumeSet files and NIfTI files; an fmrigds representation can also be requested.
## S3 method for class 'fmri_lm' write_results( x, path = NULL, subject = NULL, task = NULL, space = NULL, desc = "GLM", format = c("h5", "nifti", "gds"), strategy = c("by_stat", "by_contrast"), save_betas = TRUE, contrasts = NULL, contrast_match = c("auto", "exact", "regex"), contrast_stats = c("beta", "tstat", "pval", "se"), overwrite = FALSE, validate_inputs = TRUE, ... )## S3 method for class 'fmri_lm' write_results( x, path = NULL, subject = NULL, task = NULL, space = NULL, desc = "GLM", format = c("h5", "nifti", "gds"), strategy = c("by_stat", "by_contrast"), save_betas = TRUE, contrasts = NULL, contrast_match = c("auto", "exact", "regex"), contrast_stats = c("beta", "tstat", "pval", "se"), overwrite = FALSE, validate_inputs = TRUE, ... )
x |
An fmri_lm object containing fitted model results |
path |
Output directory path. If NULL, uses current working directory |
subject |
Subject identifier (e.g., "01", "1001"). Required. |
task |
Task identifier (e.g., "nback", "rest"). Required for BIDS compliance. |
space |
Spatial reference (e.g., "MNI152NLin2009cAsym"). Optional but recommended. |
desc |
BIDS |
format |
Output format(s). One or more of |
strategy |
Storage strategy: "by_stat" (group contrasts by statistic) or "by_contrast" (separate files) |
save_betas |
Logical. Save raw regressor betas (default: TRUE) |
contrasts |
Character vector of contrast names to save. NULL saves all contrasts |
contrast_match |
How to match |
contrast_stats |
Character vector of contrast statistics to save (default: c("beta", "tstat", "pval", "se")) |
overwrite |
Logical. Overwrite existing files (default: FALSE) |
validate_inputs |
Logical. Validate fmrilm object structure (default: TRUE) |
... |
Additional arguments passed to internal functions |
Invisible list of file paths created
## Not run: # Save all results using default settings write_results(fitted_model, subject = "01", task = "nback") # Save only specific contrasts and statistics write_results(fitted_model, subject = "01", task = "nback", space = "MNI152NLin2009cAsym", contrasts = c("FacesVsPlaces", "GoVsNoGo"), contrast_stats = c("beta", "tstat")) ## End(Not run)## Not run: # Save all results using default settings write_results(fitted_model, subject = "01", task = "nback") # Save only specific contrasts and statistics write_results(fitted_model, subject = "01", task = "nback", space = "MNI152NLin2009cAsym", contrasts = c("FacesVsPlaces", "GoVsNoGo"), contrast_stats = c("beta", "tstat")) ## End(Not run)
Exports spatial meta-analysis maps with the same core export contract used by
write_results.fmri_lm: BIDS-style filenames, HDF5 and NIfTI image
outputs, JSON metadata sidecars, vectorized format, explicit overwrite
handling, and atomic finalization.
## S3 method for class 'fmri_meta' write_results( x, path = NULL, subject = NULL, task = NULL, space = NULL, desc = "Meta", format = c("h5"), strategy = c("by_stat", "by_coefficient"), coefficients = NULL, coefficient_match = c("auto", "exact", "regex"), coefficient_stats = c("beta", "se", "z", "pval"), heterogeneity = TRUE, overwrite = FALSE, validate_inputs = TRUE, ... )## S3 method for class 'fmri_meta' write_results( x, path = NULL, subject = NULL, task = NULL, space = NULL, desc = "Meta", format = c("h5"), strategy = c("by_stat", "by_coefficient"), coefficients = NULL, coefficient_match = c("auto", "exact", "regex"), coefficient_stats = c("beta", "se", "z", "pval"), heterogeneity = TRUE, overwrite = FALSE, validate_inputs = TRUE, ... )
x |
An fmri_meta object |
path |
Output directory path. If NULL, uses current working directory. |
subject |
Optional subject or group identifier. Unlike subject-level GLM exports, meta-analysis outputs may omit this for group-level maps. |
task |
Optional task identifier. |
space |
Optional spatial reference label. |
desc |
Description of the analysis (default: "Meta"). |
format |
Output format(s). Use |
strategy |
Storage strategy. |
coefficients |
Character vector of coefficient names to save. NULL saves all coefficients. |
coefficient_match |
How to match |
coefficient_stats |
Character vector of coefficient statistics to save.
Supported values are |
heterogeneity |
Logical. Save heterogeneity maps ( |
overwrite |
Logical. Overwrite existing files (default: FALSE). |
validate_inputs |
Logical. Validate fmri_meta object structure (default: TRUE). |
... |
Additional arguments passed to internal functions. |
Invisible list of created files
Back-transform Fisher's Z to Correlation
z_to_r(z)z_to_r(z)
z |
Numeric vector or matrix; Fisher's Z values |
Numeric vector or matrix; correlations
z_to_r(0.2)z_to_r(0.2)
Compute Z-scores from Meta-Analysis
zscores(object)zscores(object)
object |
An fmri_meta object |
Matrix of z-scores
toy_meta <- structure( list( coefficients = matrix(c(0.2, -0.1), nrow = 1), se = matrix(c(0.05, 0.08), nrow = 1) ), class = "fmri_meta" ) zscores(toy_meta)toy_meta <- structure( list( coefficients = matrix(c(0.2, -0.1), nrow = 1), se = matrix(c(0.05, 0.08), nrow = 1) ), class = "fmri_meta" ) zscores(toy_meta)