Performs spatial Bayesian GLM on the cortical surface for fMRI task activation
BayesGLM_cifti(
cifti_fname,
surfL_fname = NULL,
surfR_fname = NULL,
brainstructures = c("left", "right"),
design = NULL,
onsets = NULL,
TR = NULL,
nuisance = NULL,
dHRF = c(0, 1, 2),
dHRF_as = c("auto", "nuisance", "task"),
hpf = NULL,
DCT = if (is.null(hpf)) {
4
} else {
NULL
},
resamp_res = 10000,
task_names = NULL,
session_names = NULL,
combine_sessions = TRUE,
scale_BOLD = c("auto", "mean", "sd", "none"),
scale_design = TRUE,
Bayes = TRUE,
ar_order = 6,
ar_smooth = 5,
aic = FALSE,
num.threads = 4,
return_INLA = c("trimmed", "full", "minimal"),
verbose = 1,
meanTol = 1e-06,
varTol = 1e-06
)An object of class "BayesGLM_cifti": a list with elements
The task coefficients for the Bayesian model.
The task 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.
The names of the sessions.
The number of sessions (before averaging, if applicable).
The task part of the design matrix, after centering and scaling, but before any nuisance regression or prewhitening.
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.
Left cortex surface geometry in GIFTI format
("*.surf.gii"). This can be a file path to a GIFTI file or a "surf"
object from the ciftiTools package. This argument is only used if
brainstructures includes "left" and Bayes==TRUE. If
it's not provided, the HCP group-average inflated surface included in the
ciftiTools package will be used.
Right cortex surface geometry in GIFTI format
("*.surf.gii"). This can be a file path to a GIFTI file or a "surf"
object from the ciftiTools package. This argument is only used if
brainstructures includes "right" and Bayes==TRUE. If
it's not provided, the HCP group-average inflated surface included in the
ciftiTools package will be used.
Character vector indicating which brain structure(s)
to analyze: "left" (left cortical surface) and/or "right"
(right cortical surface). Default: c("left","right") (both
hemispheres). Note that the subcortical models have not yet been implemented.
Either provide design directly, or provide
both onsets and TR from which the design matrix or matrices
will be constructed.
design is a \(T \times K\) task design matrix. Each column
represents the expected BOLD response due to each task, a convolution of
the hemodynamic response function (HRF) and the task stimulus. Note that
the scale of the regressors will affect the scale and interpretation of the
beta coefficients, so imposing a proper scale is recommended; see the
scale_design argument, which by default is TRUE. Task names
should be the column names, if not provided by the task_names
argument. For multi-session modeling, this argument should be a list of
such matrices. To model HRF derivatives, calculate the derivatives of the
task columns beforehand (see the helper function cderiv which
computes the discrete central derivative) and either add them to
design to model them as tasks, or nuisance to model them as
nuisance signals; it's recommended to then drop the first and last
timepoints because the discrete central derivative doesn't exist at the
time series boundaries. Do note that INLA computation times increase
greatly if the design matrix has more than five columns, so it might be
required to add these derivatives to nuisance rather than
design.
onsets is an \(L\)-length list in which the name of each element is
the name of the corresponding task, and the value of each element is a
matrix of onsets (first column) and durations (second column) for each
stimuli (each row) of the corresponding task. The units of both columns
is seconds. For multi-session modeling, this argument should be a list of
such lists. To model HRF derivatives, use the arguments dHRF and
dHRF_as. If dHRF==0 or dHRF_as=="nuisance", the total
number of columns in the design matrix, \(K\), will equal \(L\).
If dHRF_as=="task", \(K\) will equal \(L\) times dHRF+1.
TR is the temporal resolution of the data, in seconds.
(Optional) A \(T \times J\) matrix of nuisance signals. These are regressed from the fMRI data and the design matrix prior to the GLM computation. For multi-session modeling, this argument should be a list of such matrices.
Only applicable if onsets and TR are
provided. These arguments enable the modeling of HRF derivatives.
Set dHRF to 1 to model the temporal derivatives of each task,
2 to add the second derivatives too, or 0 to not model the
derivatives. Default: 1.
If dHRF > 0, dHRF_as controls whether the derivatives are
modeled as "nuisance" signals to regress out, "tasks", or
"auto" (default) to treat them as tasks unless the total number of
columns in the design matrix would exceed five.
Add DCT bases to nuisance to apply a temporal
high-pass filter to the data? Only one of these arguments should be provided.
hpf should be the filter frequency; if it is provided, TR
must be provided too. The number of DCT bases to include will be computed
to yield a filter with as close a frequency to hpf as possible.
Alternatively, DCT can be provided to directly specify the number
of DCT bases to include.
Default: DCT=4. For typical TR, four DCT bases amounts to a
lower frequency cutoff than the approximately .01 Hz used in most studies.
We selected this default to err on the side of retaining more low-frequency
information, but we recommend setting these arguments to values most
appropriate for the data analysis at hand.
Using at least two DCT bases is as sufficient as using linear and quadratic
drift terms in the design matrix. So if DCT detrending is being used, there
is no need to add linear and quadratic drift terms to nuisance.
The number of vertices to which each cortical surface
should be resampled, or NULL to not resample. For computational
feasibility, a value of 10000 or lower is recommended.
(Optional) Names of tasks represented in design matrix.
(Optional, and only relevant for multi-session modeling)
Names of each session. Default: NULL. In BayesGLM this
argument will overwrite the names of the list entries in data, if
both exist.
If multiple sessions are provided, should their data be combined and analyzed as a single session?
If TRUE (default), the multiple sessions will be concatenated along
time after scaling and nuisance regression, but before prewhitening. If
FALSE, each session will be analyzed separately, except that a single
estimate of the AR model coefficients for prewhitening is used, estimated
across all sessions.
Option for scaling the BOLD response.
"auto" (default) will use "mean" scaling except if demeaned
data is detected (if any mean is less than one), in which case "sd"
scaling will be used instead.
"mean" scaling will scale the data to percent local signal change.
"sd" scaling will scale the data by local standard deviation.
"none" will only center the data, not scale it.
Scale the design matrix by dividing each column by its
maximum and then subtracting the mean? Default: TRUE. If
FALSE, the design matrix is centered but not scaled.
If TRUE (default), will fit a spatial Bayesian GLM in
addition to the classical GLM. (The classical GLM is always returned.)
(numeric) Controls prewhitening. If greater than zero, this
should be a number indicating the order of the autoregressive model to use
for prewhitening. If zero, do not prewhiten. Default: 6. For
multi-session models, note that a single AR model is used; the parameters
are estimated by averaging the estimates from each session.
(numeric) FWHM parameter for smoothing the AR model
coefficient estimates for prewhitening. Remember that
\(\sigma = \frac{FWHM}{2*sqrt(2*log(2)}\). Set to 0 or NULL
to not do any smoothing. Default: 5.
Use the AIC to select AR model order 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) to return only the more relevant results, which
is enough for both id_activations and BayesGLM2,
"minimal" to return just enough for BayesGLM2 but not
id_activations, or "full" to return the full output of
inla.
Should updates be printed? Use 1 (default) for
occasional updates, 2 for occasional updates as well as running INLA
in verbose mode (if applicable), or 0 for no 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.
INLA computation times increase greatly when the number of columns in the
design matrix exceeds five. So if there are more than five tasks, or
three or more tasks each with its temporal derivative being modeled as a
task, BayesGLM will raise a warning. In cases like the latter, we
recommend modeling the temporal derivatives as nuisance signals using the
nuisance argument, rather than modeling them as tasks.
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.