Performs spatial Bayesian GLM for task fMRI activation with CIFTI-format data. The cortex is modeled as a surface mesh, and subcortical structures are modeled as distinct volumetric regions. Includes the pre-processing steps of nuisance regression, prewhitening, scaling, and variance normalization. Supports both single- and multi-session analysis. Can also compute just the classical (spatially-independent)
BayesGLM(
BOLD,
brainstructures = c("left", "right"),
subROI = c("Amygdala-L", "Amygdala-R", "Caudate-L", "Caudate-R", "Hippocampus-L",
"Hippocampus-R", "Thalamus-L", "Thalamus-R"),
design,
nuisance = NULL,
scrub = NULL,
hpf = NULL,
TR = NULL,
surfL = NULL,
surfR = NULL,
resamp_res = 10000,
nbhd_order = 1,
buffer = c(1, 1, 3, 4, 4),
session_names = NULL,
scale_BOLD = c("mean", "sd", "none"),
Bayes = TRUE,
hyperpriors = c("informative", "default"),
ar_order = 6,
ar_smooth = 5,
aic = FALSE,
n_threads = 4,
return_INLA = c("trimmed", "full", "minimal"),
verbose = 1,
meanTol = 1e-06,
varTol = 1e-06
)An object of class "BayesGLM": a list with elements
The field coefficients for the Bayesian model.
The field coefficients for the classical model.
The entire list of GLM results, except for parameters estimated for the classical model.
Parameters estimated for the classical model from the GLM.
data.frame summarizing the spatial features of each brain structure modeled.
data.frame with the name and nTime of each BOLD session.
data.frame with the name, related task, and HRF_order of each field.
fMRI timeseries data in CIFTI format ("*.dtseries.nii").
For single-session analysis this can be a file path to a CIFTI file or a
"xifti" object from the ciftiTools package. For multi-session
analysis this can be a vector of file paths or a list of "xifti"
objects.
If BOLD is a "xifti" object(s), the surfaces, if any, will be
used for the spatial model. However, if surfL and surfR are
provided, they will override any surfaces in BOLD.
Character vector indicating which brain structure(s)
of BOLD to analyze: "left" cortex; "right" cortex;
and/or "subcortical" structures. Or "all" to model all three.
Default: c("left","right") (cortex only).
Which subcortical ROIs should be analyzed? Can be "all"
to analyze all subcortex ROIs. See the ciftiTools_Name column of
ciftiTools:::substructure_table() for a list of possible
subcortical ROIs.
A numeric matrix or data.frame, or a
"BayesfMRI_design" object from make_design. Can also
be an array where the third dimension is the same length as the number of
data locations, to model each location with its own design.
(Optional) A \(T \times N_{nuis}\) matrix of nuisance signals, where \(T\) is the number of timepoints and \(N\) is the number of nuisance signals, or a list of these for multi-session analysis. Nuisance signals are regressed from the fMRI data and design matrix prior to GLM computation. Nuisance signals can include motion regressors, HRF derivatives not being modeled as tasks, and other sources of noise.
Detrending/high-pass filtering is accomplished by adding DCT bases to the
nuisance matrix; see the parameters hpf and DCT.
Do not add spike regressors for scrubbing to the nuisance matrix.
Rather, provide these in scrub so that their corresponding timepoints
are also removed from the BOLD data after nuisance regression.
(Optional) A \(T \times N_{scrub}\) matrix of spike regressors
(one 1 value at the timepoint to scrub, and 0 for all other values), or a
logical vector indicating the timepoints to scrub (TRUE to scrub, and
FALSE to keep). For multi-session data, a session-length list of
such matrices or logical vectors.
The spike regressors will be included in the nuisance
regression, and afterwards the timepoints indicated in scrub will be
removed from the BOLD data and design matrix.
Add DCT bases to nuisance to apply a temporal high-pass
filter to the data, for detrending? hpf is the filter frequency.
Use NULL to skip detrending. Detrending is strongly recommended for
fMRI data, to help reduce the autocorrelation in the residuals, so
NULL will induce a warning. Use "already" to disable the
warning while skipping highpass filtering.
Using at least two DCT bases is as sufficient for detrending as using linear
and quadratic drift terms in the nuisance matrix. So if DCT detrending is
being used here, there is no need to add linear and quadratic drift terms to
nuisance.
Temporal resolution of the data, in seconds.
For cortex spatial model. Left and right cortex surface
geometry in GIFTI format ("*.surf.gii"). These can be a file path to
a GIFTI file or a "surf" object from ciftiTools.
Surfaces can alternatively be provided through the $surf metadata in
BOLD if it is "xifti" data. If neither are provided, by default the
HCP group-average fs_LR inflated surfaces included in ciftiTools will be
used for the cortex spatial model.
For cortex spatial model. The number of vertices to which
each cortical surface should be resampled, or NULL to not resample.
For computational feasibility, a value of 10000 (default) or lower is
recommended for Bayesian spatial modeling. If Bayes=FALSE,
resamp_res can be set to NULL for full-resolution classical
modeling.
For volumetric model. What order neighborhood around data
locations to keep? 0 for no neighbors, 1 for 1st-order
neighbors, 2 for 1st- and 2nd-order neighbors, etc. Smaller values
will provide greater computational efficiency at the cost of higher variance
around the edge of the data.
For volumetric model. The number of extra voxel layers around
the bounding box. Set to NULL for no buffer. (We recommend not
changing buffer unless you know what you're doing. Instead, to reduce
the number of boundary voxels, adjust nbhd_order).
The names of the task-fMRI BOLD sessions, for
multi-session analysis. If not provided here, will be inferred from
names(BOLD), inferred from names(design), or generated
automatically, in that order.
Controls scaling the BOLD response at each location.
Scale the data to percent local signal change.
Scale the data by local standard deviation.
Center the data but do not scale it.
Perform spatial Bayesian modeling? Default: TRUE. If
FALSE, only perform classical (massive univariate) modeling. (The classical GLM
result is always returned, whether Bayes is TRUE or FALSE.)
Should informative or default non-informative hyperpriors be assumed on SPDE hyperparameters?
(For prewhitening) The order of the autoregressive (AR) model
to use for prewhitening. If 0, do not prewhiten. Default: 6.
For multi-session modeling, note that a single AR model is used; its coefficients will be the average estimate from each session.
(For prewhitening) The FWHM parameter for spatially
smoothing the coefficient estimates for the AR model to use for
prewhitening. Recall that
\(\sigma = \frac{FWHM}{2*sqrt(2*log(2)}\). Set to 0 to not smooth
the estimates. Default: 5.
(For prewhitening) Use the Akaike information criterion (AIC) to
select AR model orders between 0 and ar_order? Default:
FALSE.
The maximum number of threads to use for parallel
computations: prewhitening parameter estimation, and the inla-program model
estimation. Default: 4. Note that parallel prewhitening requires the
parallel package.
Return the INLA model object? (It can be large.) Use
"trimmed" (default) returns the results sufficient for
activations and BayesGLM2; "minimal"
returns enough for BayesGLM2 but not
activations; "full" returns the full inla
output.
1 (default) to print occasional updates during model
computation; 2 for occasional updates as well as running INLA in
verbose mode (if Bayes), or 0 for no printed updates.
Tolerance for mean and variance of each data location.
Locations which do not meet these thresholds are masked out of the analysis.
Default: 1e-6 for both.
This function uses a system wrapper for the 'wb_command' executable. The user must first download and install the Connectome Workbench, available from https://www.humanconnectome.org/software/get-connectome-workbench .
This function requires the INLA package, which is not a CRAN package.
See https://www.r-inla.org/download-install for easy installation instructions.
INLA computation times increase greatly when the number of columns in the
design matrix exceeds five: when there are more than five tasks, or more
than three tasks each with a temporal derivative modeled as a field. In
cases like the latter, we recommend modeling the temporal derivatives as
nuisance signals using the option dHRF_as="nuisance", rather than
modeling the temporal derivatives as fields.
To use BayesGLM, the design matrix must first be constructed
with make_design.