Type: Package
Title: Spatial Bayesian Methods for Task Functional MRI Studies
Version: 0.10.1
Maintainer: Amanda Mejia <mandy.mejia@gmail.com>
Description: Performs a spatial Bayesian general linear model (GLM) for task functional magnetic resonance imaging (fMRI) data on the cortical surface. Additional models include group analysis and inference to detect thresholded areas of activation. Includes direct support for the 'CIFTI' neuroimaging file format. For more information see A. F. Mejia, Y. R. Yue, D. Bolin, F. Lindgren, M. A. Lindquist (2020) <doi:10.1080/01621459.2019.1611582> and D. Spencer, Y. R. Yue, D. Bolin, S. Ryan, A. F. Mejia (2022) <doi:10.1016/j.neuroimage.2022.118908>.
License: GPL-3
URL: https://github.com/mandymejia/BayesfMRI
BugReports: https://github.com/mandymejia/BayesfMRI/issues
Additional_repositories: https://inla.r-inla-download.org/R/testing
Depends: R (≥ 3.6.0)
Imports: car, ciftiTools (≥ 0.17.2), excursions, foreach, fMRItools (≥ 0.5.3), MASS, Matrix, matrixStats, methods, Rcpp, stats, sp, utils, viridisLite
Suggests: covr, abind, grDevices, hrf, INLA (≥ 0.0-1468840039), knitr, MatrixModels, parallel, purrr, rmarkdown, SQUAREM, testthat (≥ 3.0.0), spelling
Encoding: UTF-8
RoxygenNote: 7.3.2
LinkingTo: RcppEigen, Rcpp
Language: en-US
NeedsCompilation: yes
Packaged: 2025-03-07 23:13:18 UTC; ddpham
Author: Amanda Mejia [aut, cre], Damon Pham ORCID iD [ctb], David Bolin [ctb], Yu (Ryan) Yue [ctb], Daniel Spencer ORCID iD [aut], Sarah Ryan [ctb]
Repository: CRAN
Date/Publication: 2025-03-07 23:50:02 UTC

BayesfMRI: Spatial Bayesian Methods for Task Functional MRI Studies

Description

logo

Performs a spatial Bayesian general linear model (GLM) for task functional magnetic resonance imaging (fMRI) data on the cortical surface. Additional models include group analysis and inference to detect thresholded areas of activation. Includes direct support for the 'CIFTI' neuroimaging file format. For more information see A. F. Mejia, Y. R. Yue, D. Bolin, F. Lindgren, M. A. Lindquist (2020) doi:10.1080/01621459.2019.1611582 and D. Spencer, Y. R. Yue, D. Bolin, S. Ryan, A. F. Mejia (2022) doi:10.1016/j.neuroimage.2022.118908.

Author(s)

Maintainer: Amanda Mejia mandy.mejia@gmail.com

Authors:

Other contributors:

See Also

Useful links:


Perform the EM algorithm of the Bayesian GLM fitting

Description

Perform the EM algorithm of the Bayesian GLM fitting

Usage

.findTheta(theta, spde, y, X, QK, Psi, A, Ns, tol, verbose = FALSE)

Arguments

theta

the vector of initial values for theta

spde

a list containing the sparse matrix elements Cmat, Gmat, and GtCinvG

y

the vector of response values

X

the sparse matrix of the data values

QK

a sparse matrix of the prior precision found using the initial values of the hyperparameters

Psi

a sparse matrix representation of the basis function mapping the data locations to the mesh vertices

A

a precomputed matrix crossprod(X%*%Psi)

Ns

the number of columns for the random matrix used in the Hutchinson estimator

tol

a value for the tolerance used for a stopping rule (compared to the squared norm of the differences between theta(s) and theta(s-1))

verbose

(logical) Should intermediate output be displayed?


Get the prewhitening matrix for a single data location

Description

Get the prewhitening matrix for a single data location

Usage

.getSqrtInvCpp(AR_coefs, nTime, avg_var)

Arguments

AR_coefs

a length-p vector where p is the AR order

nTime

(integer) the length of the time series that is being prewhitened

avg_var

a scalar value of the residual variances of the AR model


Find the initial values of kappa2 and phi

Description

Find the initial values of kappa2 and phi

Usage

.initialKP(theta, spde, w, n_sess, tol, verbose)

Arguments

theta

a vector of length two containing the range and scale parameters kappa2 and phi, in that order

spde

a list containing the sparse matrix elements Cmat, Gmat, and GtCinvG

w

the beta_hat estimates for a single task

n_sess

the number of sessions

tol

the stopping rule tolerance

verbose

(logical) Should intermediate output be displayed?


Find the log of the determinant of Q_tilde

Description

Find the log of the determinant of Q_tilde

Usage

.logDetQt(kappa2, in_list, n_sess)

Arguments

kappa2

a scalar

in_list

a list with elements Cmat, Gmat, and GtCinvG

n_sess

the integer number of sessions


Corrected AIC

Description

Computes corrected AIC (AICc).

Usage

AICc(y, demean = FALSE, order.max = 10)

Arguments

y

The autocorrelated data

demean

Demean y? Default: FALSE.

order.max

The model order limit. Default: 10.

Value

The cAIC


BOLD

Description

BOLD

Arguments

BOLD

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.


BayesGLM for CIFTI

Description

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)

Usage

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
)

Arguments

BOLD

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.

brainstructures

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).

subROI

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.

design

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.

nuisance

(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.

scrub

(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.

hpf

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.

TR

Temporal resolution of the data, in seconds.

surfL, surfR

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.

resamp_res

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.

nbhd_order

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.

buffer

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).

session_names

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.

scale_BOLD

Controls scaling the BOLD response at each location.

"mean":

Scale the data to percent local signal change.

"sd":

Scale the data by local standard deviation.

"none":

Center the data but do not scale it.

Bayes

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.)

hyperpriors

Should informative or default non-informative hyperpriors be assumed on SPDE hyperparameters?

ar_order

(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.

ar_smooth

(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.

aic

(For prewhitening) Use the Akaike information criterion (AIC) to select AR model orders between 0 and ar_order? Default: FALSE.

n_threads

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_INLA

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.

verbose

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.

meanTol, varTol

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.

Details

To use BayesGLM, the design matrix must first be constructed with make_design.

Value

An object of class "BayesGLM": a list with elements

betas_Bayesian

The field coefficients for the Bayesian model.

betas_classical

The field coefficients for the classical model.

GLMs_Bayesian

The entire list of GLM results, except for parameters estimated for the classical model.

GLMs_classical

Parameters estimated for the classical model from the GLM.

brainstructures

data.frame summarizing the spatial features of each brain structure modeled.

sessions

data.frame with the name and nTime of each BOLD session.

fields

data.frame with the name, related task, and HRF_order of each field.

Connectome Workbench Requirement

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 .

INLA Requirement

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 Latent Fields Limit

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.


Group-level Bayesian GLM

Description

Performs group-level Bayesian GLM estimation and inference using the joint approach described in Mejia et al. (2020).

Usage

BayesGLM2(
  results,
  contrasts = NULL,
  quantiles = NULL,
  excursion_type = NULL,
  contrast_names = NULL,
  gamma = 0,
  alpha = 0.05,
  nsamp_theta = 50,
  nsamp_beta = 100,
  num_cores = NULL,
  verbose = 1
)

Arguments

results

Either (1) a length N list of "BGLM" objects, or (2) a length N character vector of files storing "BGLM" objects saved with saveRDS. "fit_bglm" objects also are accepted.

contrasts

(Optional) A list of contrast vectors that specify the group-level summaries of interest. If NULL (DEFAULT), use contrasts that compute the average of each field (field HRF) across all subjects/sessions.

Each contrast vector is length KSN specifying a group-level summary of interest, where K is the number of fields in the first-level design matrices, S is the number of sessions, and N is the number of subjects. The vector is grouped by fields, then sessions, then subjects.

For a single session/subject, the contrast vector for the first field would be:

c0 <- c(1, rep(0, K-1)) #indexes the first field for a single session

so the full contrast vector for the group average over all sessions/subjects for the first field would be:

contrasts = rep(c0, S*N) /(S*N).

To obtain the group average for the first field, for just the first session, input zeros for the remaining sessions:

c2 <- c(c0, rep(0, K*(S-1))) contrasts = rep(c2, N) /N.

To obtain the group mean difference between two sessions (S=2) for the first field:

c3 <- c(c0, -c0) contrasts = rep(c3, N) / N.

To obtain the mean over sessions of the first field, just for the first subject:

c4 <- rep(c0, S) c(c4, rep(0, K*S*(N-1))) / S.

quantiles

(Optional) Vector of posterior quantiles to return in addition to the posterior mean.

excursion_type

(For inference only) The type of excursion function for the contrast (">", "<", "!="), or a vector thereof (each element corresponding to one contrast). If NULL, no inference performed.

contrast_names

(Optional) Names of contrasts.

gamma

(For inference only) Activation threshold for the excursion set, or a vector thereof (each element corresponding to one contrast). Default: 0.

alpha

(For inference only) Significance level for activation for the excursion set, or a vector thereof (each element corresponding to one contrast). Default: .05.

nsamp_theta

Number of theta values to sample from posterior. Default: 50.

nsamp_beta

Number of beta vectors to sample conditional on each theta value sampled. Default: 100.

num_cores

The number of cores to use for sampling betas in parallel. If NULL (default), do not run in parallel.

verbose

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.

Value

A list containing the estimates, PPMs and areas of activation for each contrast.

INLA Requirement

This function requires the INLA package, which is not a CRAN package. See https://www.r-inla.org/download-install for easy installation instructions.


Bayes GLM arg checks

Description

Checks arguments for BayesGLM and fit_bayesglm

Usage

BayesGLM_argChecks(
  scale_BOLD = c("mean", "sd", "none"),
  Bayes = TRUE,
  EM = FALSE,
  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,
  emTol = 0.001
)

Arguments

scale_BOLD

See BayesGLM.

Bayes, EM

See BayesGLM.

ar_order, ar_smooth, aic

See BayesGLM.

n_threads

See BayesGLM.

return_INLA

See BayesGLM.

verbose

See BayesGLM.

meanTol, varTol, emTol

See BayesGLM.

Details

Avoids duplicated code between BayesGLM and fit_bayesglm

Value

The arguments that may have changed, in a list: scale_BOLD, do_Bayesian, do_EM, and do_pw.


Format fit_bayesglm results into "xifti" objects

Description

Format beta estimates or RSS of the fit_bayesglm results into list of "xifti" objects for each session.

Usage

BayesGLM_format_cifti(
  BGLMs,
  do,
  spatial,
  submeta,
  session_names,
  field_names,
  method = c("classical", "Bayesian")
)

Arguments

BGLMs

The list of "BGLM" objects, from fit_bayesglm.

do, spatial, submeta, session_names, field_names

See BayesGLM.

method

"classical" or "Bayesian"

Value

The list of "xifti" objects


Format design

Description

Format design for BayesGLM, fit_bayesglm, multiGLM, and multiGLM_fun.

Usage

BayesGLM_format_design(
  design,
  scale_design = TRUE,
  nS_expect = NULL,
  nT_expect = NULL,
  nD_expect = NULL,
  per_location_design = NULL
)

Arguments

design

The design argument input. Will be formatted to a nS-length list.

scale_design

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.

nS_expect

The expected number of sessions, if known.

nT_expect

The expected number of timepoints, if known. For multi-session data this is a session-length vector.

nD_expect

The expected number of designs, if known. For per-location modeling this is equal to nVd0, the initial number of data locations. For multi-session data this is a session-length vector.

per_location_design

FALSE if per-location modeling is not being performed (i.e. for multiGLM); TRUE if it is; or, NULL to infer based on the dimensions of design (TRUE if the design has three dimensions.)

Value

design


Format nuisance

Description

Format nuisance for BayesGLM, fit_bayesglm, multiGLM, and multiGLM_fun.

Usage

BayesGLM_format_nuisance(nuisance, nS_expect = NULL, nT_expect = NULL)

Arguments

nuisance

The nuisance argument input. Will be formatted to a nS-length list.

nS_expect

The expected number of sessions, if known.

nT_expect

The expected number of timepoints, if known. For multi-session data this is a session-length vector.

Value

nuisance


Format scrub

Description

Format scrub for BayesGLM, fit_bayesglm, multiGLM, and multiGLM_fun.

Usage

BayesGLM_format_scrub(scrub, nS_expect = NULL, nT_expect = NULL)

Arguments

scrub

The scrub argument input. Will be formatted to a nS-length list.

nS_expect

The expected number of sessions, if known.

nT_expect

The expected number of timepoints, if known. For multi-session data this is a session-length vector.

Value

scrub


Is a valid design?

Description

Is a valid design?

Usage

BayesGLM_is_valid_one_design(design)

Arguments

design

The design matrix/array


Is a valid nuisance?

Description

Is a valid nuisance?

Usage

BayesGLM_is_valid_one_nuisance(nuisance)

Arguments

nuisance

The nuisance matrix


Is a valid scrub?

Description

Is a valid scrub?

Usage

BayesGLM_is_valid_one_scrub(scrub)

Arguments

scrub

The scrub matrix


Get session_names for GLM

Description

Get session_names for GLM

Usage

BayesGLM_session_names(nS, session_names, BOLD_names, design_names)

Bayes

Description

Bayes

Arguments

Bayes

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.)


Connectome Workbench

Description

Connectome Workbench

Connectome Workbench Requirement

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 .


Expected log-likelihood function

Description

Expected log-likelihood function

Usage

ELL(Q, sigma2, model_data, Psi, mu, Sigma, A)

Arguments

Q

precision matrix

sigma2

noise variance

model_data

list with X and y

Psi

data locations to mesh locations matrix

mu

posterior mean of w

Sigma

posterior covariance of w

A

crossprod(X%*%Psi)

Value

scalar expected log-likelihood


EM

Description

EM

Arguments

EM

(logical) Should the EM implementation of the Bayesian GLM be used? Default: FALSE. This method is still in development.


F logwt

Description

Internal function used in joint approach to group-analysis for combining across models

Usage

F.logwt(theta, spde, mu_theta, Q_theta, nN)

Arguments

theta

A vector of hyperparameter values at which to compute the posterior log density

spde

A SPDE object from inla.spde2.matern() function, determines prior precision matrix

mu_theta

Posterior mean from combined subject-level models.

Q_theta

Posterior precision matrix from combined subject-level models.

nN

Number of subjects

Value

The prior density

INLA Requirement

This function requires the INLA package, which is not a CRAN package. See https://www.r-inla.org/download-install for easy installation instructions.


Fixed point function for the joint BayesGLM EM update algorithm

Description

Fixed point function for the joint BayesGLM EM update algorithm

Usage

GLMEM_fixptseparate(theta, spde, model_data, Psi, K, A, cl, Ns = 50)

Arguments

theta

a list containing kappa2, phi, and sigma2, in that order

spde

the spde object

model_data

the model_data object containing y and X

Psi

a conversion matrix (N by V) (or N by n)

K

number of covariates

A

The value for Matrix::crossprod(X%*%Psi) (saves time on computation)

cl

parallelization cluster

Ns

The number of samples used to approximate traces using the Hutchinson estimator. If set to 0, the exact trace is found.

Value

a vector with the same length as theta, the EM updates


Objective function for the BayesGLM EM algorithm

Description

This returns the negative of the expected log-likelihood function.

Usage

GLMEM_objfn(theta, spde, model_data, Psi, K, A, n_threads = NULL, Ns = NULL)

Arguments

theta

a vector containing kappa2, phi, and sigma2, in that order

spde

the spde object

model_data

the model_data object containing y and X

Psi

a conversion matrix (N by V) (or N by n)

K

number of covariates

A

The value for Matrix::crossprod(X%*%Psi) (saves time on computation)

n_threads

Needed for SQUAREM (it is an argument to the fixed-point functions)

Ns

The number of samples used to approximate traces using the Hutchinson estimator. If set to 0, the exact trace is found.

Value

A scalar value for the negative expected log-likelihood


Classical GLM

Description

Classical GLM for fit_bayesglm (internal function)

Usage

GLM_classical(
  BOLD,
  design,
  nK2,
  nV_input,
  field_names,
  design_type,
  valid_cols,
  nT,
  do_pw,
  compute_SE = TRUE
)

Arguments

BOLD

BOLD timeseries in vector form (TVx1), result of sparse_and_PW

design

List of large sparse design matrices (each TVxV), one per regressor, result of sparse_and_PW

nK2, nV_input, field_names, design_type

See fit_bayesglm.

valid_cols, nT

See fit_bayesglm.

do_pw

Has prewhitening been performed on the data and design?

compute_SE

Compute SE of model coefficients?

Value

A list of results


Classical GLM for multiple models

Description

Classical GLM for multiple models

Usage

GLM_compare(...)

Standardize data variance, and prewhiten if applicable

Description

Standardize data variance and prewhiten if applicable, for the GLM.

Usage

GLM_est_resid_var_pw(
  BOLD,
  design,
  spatial,
  session_names,
  field_names,
  design_type,
  valid_cols,
  nT,
  ar_order,
  ar_smooth,
  aic,
  n_threads,
  do_pw
)

Arguments

BOLD, design, spatial

See fit_bayesglm.

session_names, field_names, design_type

See fit_bayesglm.

valid_cols, nT, do_pw

See fit_bayesglm.

Value

List of results


INLA

Description

INLA

INLA Requirement

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 Latent Fields

Description

INLA Latent Fields

INLA Latent Fields Limit

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.


Import INLA dependencies

Description

Roxygen entry for additional INLA dependencies


Q prime

Description

Gives the portion of the Q matrix independent of phi

Usage

Q_prime(kappa2, spde)

Arguments

kappa2

scalar

spde

spde object

Value

a dgCMatrix


SPDE from mesh model

Description

SPDE from mesh model

Usage

SPDE_from_vertex(spatial, qc_mask, logkappa = NULL, logtau = NULL)

Arguments

spatial

See BayesGLM.

qc_mask

QC mask.

logkappa, logtau

vector of min, max and initial value for prior on log kappa and log tau. Min and max are extreme quantiles, not hard constraints.

Value

List: spde and spatial.


SPDE from voxel model

Description

SPDE from voxel model

Usage

SPDE_from_voxel(spatial, qc_mask = NULL, logkappa = NULL, logtau = NULL)

Arguments

spatial

See BayesGLM.

qc_mask

QC mask.

logkappa, logtau

vector of min, max and initial value for prior on log kappa and log tau. Min and max are extreme quantiles, not hard constraints.

Value

List: spde and spatial.


TR

Description

TR

Arguments

TR

Temporal resolution of the data, in seconds.


Trace approximation function

Description

Trace approximation function

Usage

TrQEww(kappa2, spde, P, mu, Vh)

Arguments

kappa2

a scalar

spde

spde object

P

Matrix of dimension nk by N_s found by solve(Sig_inv,Vh)

mu

posterior mean

Vh

matrix of random variables with nrow(Sig_inv) rows and Ns columns

Value

a scalar


Trace of Q beta' beta

Description

Trace of Q beta' beta

Usage

TrQbb(kappa2, beta_hat, spde)

Arguments

kappa2

scalar

beta_hat

a vector

spde

an spde object

Value

a scalar


Hutchinson estimator of the trace

Description

Hutchinson estimator of the trace

Usage

TrSigB(P, B, vh)

Arguments

P

Matrix of dimension nk by N_s found by solve(Sig_inv,Vh)

B

Matrix of dimension nk by nk inside the desired trace product

vh

Matrix of dimension nk by N_s in which elements are -1 or 1 with equal probability.

Value

a scalar estimate of the trace of Sigma %*% B


Identify field activations

Description

Identify areas of activation for each field from the result of BayesGLM or fit_bayesglm.

Usage

activations(
  x,
  Bayes = TRUE,
  gamma = NULL,
  alpha = 0.05,
  correction = c("FWER", "FDR", "none"),
  fields = NULL,
  sessions = NULL,
  verbose = 1
)

id_activations(
  x,
  Bayes = TRUE,
  gamma = NULL,
  alpha = 0.05,
  correction = c("FWER", "FDR", "none"),
  fields = NULL,
  sessions = NULL,
  verbose = 1
)

Arguments

x

Result of BayesGLM or fit_bayesglm model call, of class "BGLM" or "fit_bglm".

Bayes

Use spatial Bayesian modeling to identify activations based on the joint posterior distribution? Default: TRUE. If FALSE, activations will be based on classical (massive univariate) GLM model, with multiple comparisons correction (see correction). Note that TRUE is only applicable if x includes Bayesian results (i.e. x <- BayesGLM(..., Bayes = TRUE) was run.)

gamma

Activation threshold, for example 1 for 1 percent signal change if scale_BOLD=="mean" during model estimation. Setting a gamma is required for the Bayesian method; NULL (default) will use a gamma of zero for the classical method.

alpha

Significance level for inference. Default: 0.05.

correction

For the classical method only: Type of multiple comparisons correction: "FWER" (Bonferroni correction, the default), "FDR" (Benjamini Hochberg), or "none".

fields

The field(s) to identify activations for. Give either the name(s) as a character vector, or the numerical indices. If NULL (default), analyze all fields.

sessions

The session(s) to identify activations for. Give either the name(s) as a character vector, or the numerical indices. If NULL (default), analyze the first session.

verbose

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.

Value

An "act_BGLM" or "act_fit_bglm" object, a list which indicates the activated locations along with related information.


Identification of areas of activation in a General Linear Model using classical methods

Description

Identification of areas of activation in a General Linear Model using classical methods

Usage

activations.classical(
  x,
  gamma = 0,
  alpha = 0.05,
  correction = c("FWER", "FDR", "none"),
  fields,
  session,
  mesh = NULL
)

Arguments

x

A BayesGLM object

gamma, alpha, correction, fields, session

See activations.

mesh

(Optional) An "inla.mesh" object (see make_mesh for surface data). Only necessary for computing surface areas of identified activations.

Value

A matrix corresponding to the 0-1 activation status for the model coefficients.


Identify activations using joint posterior probabilities

Description

Identifies areas of activation given an activation threshold and significance level using joint posterior probabilities

Usage

activations.posterior(x, fields, session, alpha = 0.05, gamma)

Arguments

x

Result of BayesGLM, of class "BGLM".

fields, session, alpha, gamma

See activations.

Details

For a given latent field, identifies locations that exceed a certain activation threshold (e.g. 1 percent signal change) at a given significance level, based on the joint posterior distribution of the latent field.

Value

A list with two elements: active, which gives a matrix of zeros and ones of the same dimension as x$field_estimates${session}, and excur_result, an object of class "excurobj" (see INLA's excursions.inla for more information).


aic

Description

aic

Arguments

aic

(For prewhitening) Use the Akaike information criterion (AIC) to select AR model orders between 0 and ar_order? Default: FALSE.


ar_order

Description

ar_order

Arguments

ar_order

(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.


ar_smooth

Description

ar_smooth

Arguments

ar_smooth

(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.


Beta posterior theta sampling

Description

Internal function used in joint approach to group-analysis

Usage

beta.posterior.thetasamp(
  theta,
  spde,
  Xcros,
  Xycros,
  contrasts,
  quantiles,
  excursion_type,
  gamma,
  alpha,
  nsamp_beta = 100
)

Arguments

theta

A single sample of theta (hyperparameters) from q(theta|y)

spde

A SPDE object from inla.spde2.matern() function.

Xcros

A crossproduct of design matrix.

Xycros

A crossproduct of design matrix and BOLD y.

contrasts

A list of vectors of length M*K specifying the contrasts of interest.

quantiles

Vector of posterior quantiles to return in addition to the posterior mean

excursion_type

Vector of excursion function type (">", "<", "!=") for each contrast

gamma

Vector of activation thresholds for each contrast

alpha

Significance level for activation for the excursion sets

nsamp_beta

The number of samples to draw from full conditional of beta given the current value of theta (p(beta|theta,y))

Value

A list containing mu, quantiles, and F

INLA Requirement

This function requires the INLA package, which is not a CRAN package. See https://www.r-inla.org/download-install for easy installation instructions.


brainstructures

Description

brainstructures

Arguments

brainstructures

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).


buffer

Description

buffer

Arguments

buffer

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).


cbind if first argument might be NULL

Description

cbind, but return the second argument if the first is NULL

Usage

cbind2(mat_or_NULL, to_add)

Arguments

mat_or_NULL

NULL or a numeric matrix

to_add

A numeric matrix with the same number of rows as mat_or_NULL

Value

cbind(mat_or_NULL, to_add), or just to_add if the first argument is NULL.


Check INLA and PARDISO

Description

Check INLA and PARDISO

Usage

check_INLA(require_PARDISO = FALSE)

Arguments

require_PARDISO

Is PARDISO required? Default: FALSE.

Value

NULL, invisibly


Sample from the multivariate normal distribution with Cholesky(Q)

Description

Sample from the multivariate normal distribution with Cholesky(Q)

Usage

cholQsample(n, mu, cholQ)

Arguments

n

number of samples

mu

mean vector

cholQ

Cholesky decomposition of the precision (found via Matrix::Cholesky(Q))

Value

An n \times p matrix of samples from the MVN distribution, where p is the length of mu.


contrasts

Description

contrasts

Arguments

contrasts

List of contrast vectors to be passed to inla::inla.


Function to prepare objects for use in Rcpp functions

Description

Function to prepare objects for use in Rcpp functions

Usage

create_listRcpp(spde)

Arguments

spde

an spde object

Value

The SPDE matrices with the correct data formats


design

Description

design

Arguments

design

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.


Set column values to zero for sparse matrix

Description

Set column values to zero for sparse matrix

Usage

dgCMatrix_cols_to_zero(mat, cols)

Arguments

mat

dgCMatrix

cols

The column indices to set to zero

Value

The modified sparse matrix


Mask out invalid data

Description

Mask out data locations that are invalid (missing data, low mean, or low variance) for any session.

Usage

do_QC(BOLD, meanTol = 1e-06, varTol = 1e-06, verbose = TRUE)

Arguments

BOLD

A session-length list of T \times V numeric BOLD data.

meanTol, varTol

Tolerance for mean and variance of each data location. Locations which do not meet these thresholds are masked out of the analysis. Defaults: 1e-6.

verbose

Print messages counting how many locations are removed? Default: TRUE.

Value

A logical vector indicating locations that are valid across all sessions.

Examples

nT <- 30
nV <- 400
BOLD1 <- matrix(rnorm(nT*nV), nrow=nT)
BOLD1[,seq(30,50)] <- NA
BOLD2 <- matrix(rnorm(nT*nV), nrow=nT)
BOLD2[,65] <- BOLD2[,65] / 1e10
BOLD <- list(sess1=BOLD1, sess2=BOLD2)
do_QC(BOLD)


emTol

Description

emTol

Arguments

emTol

The stopping tolerance for the EM algorithm. Default: 1e-3.


Extract Estimates of Activation

Description

Obtains the posterior mean or other summary statistic for each latent field

Usage

extract_estimates(INLA_model_obj, session_names, spatial, spde, stat = "mean")

Arguments

INLA_model_obj

An object of class "inla", a result of a call to inla.

session_names

Vector of fMRI session names

stat

A string representing the posterior summary statistic to be returned

Value

Estimates from inla model


faces

Description

faces

Arguments

faces

An F \times 3 matrix, where each row contains the vertex indices for a given triangular face in the mesh. F is the number of faces in the mesh.


field_names

Description

field_names

Arguments

field_names

(Optional) Names of fields represented in design matrix.


fit_bayesglm

Description

Performs spatial Bayesian GLM for task fMRI activation

Usage

fit_bayesglm(
  BOLD,
  design,
  nuisance = NULL,
  scrub = NULL,
  spatial,
  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
)

Arguments

BOLD, design, nuisance

Session-length list of numeric matrices/arrays, each with volumes along the first dimension.

scrub

Session-length list of spike regressors: numeric matrices, with volumes along the first dimension, valued at 1 for scrubbed volumes and 0 otherwise.

Scrubbing is performed by incorporating spike regressors in the nuisance matrix during nuisance regression (in a simultaneous framework), and then removing the scrubbed timepoints from the resulting BOLD and design.

spatial

Gives the spatial information:

surf

A list of two: vertices V \times 3 numeric matrix of vertex locations in XYZ coordinate space, and faces, F \times 3 matrix of positive integers defining the triangular faces.

mask

Mask of locations with valid data.

For voxel data, a list of six:

label

3D array of labeled locations to include in the model.

trans_mat

Projection matrix to convert voxel indices to XYZ position. Can be NULL.

trans_units

XYZ units. Can be NULL.

nbhd_order

See documentation for BayesGLM.

buffer

See documentation for BayesGLM.

scale_BOLD

Controls scaling the BOLD response at each location.

"mean":

Scale the data to percent local signal change.

"sd":

Scale the data by local standard deviation.

"none":

Center the data but do not scale it.

Bayes

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.)

hyperpriors

Should informative or default non-informative hyperpriors be assumed on SPDE hyperparameters?

ar_order

(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.

ar_smooth

(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.

aic

(For prewhitening) Use the Akaike information criterion (AIC) to select AR model orders between 0 and ar_order? Default: FALSE.

n_threads

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_INLA

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.

verbose

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.

meanTol, varTol

Tolerance for mean, variance and SNR of each data location. Locations which do not meet these thresholds are masked out of the analysis. Default: 1e-6 for mean and variance, 50 for SNR.

Value

A "BayesGLM" object: a list with elements

INLA_model_obj

The full result of the call to INLA::inla.

field_estimates

The estimated coefficients for the Bayesian model.

result_classical

Results from the classical model: field estimates, field standard error estimates, residuals, degrees of freedom, and the mask.

mesh

The model mesh.

mask

A mask of mesh indicating the locations inside mesh.

design

The design matrix, after centering and scaling, but before any nuisance regression or prewhitening.

field_names

The names of the fields.

session_names

The names of the sessions.

hyperpar_posteriors

Hyperparameter posterior densities.

theta_estimates

Theta estimates from the Bayesian model.

posterior_Sig_inv

For joint group modeling.

mu_theta

For joint group modeling.

Q_theta

For joint group modeling.

y

For joint group modeling: The BOLD data after any centering, scaling, nuisance regression, or prewhitening.

X

For joint group modeling: The design matrix after any centering, scaling, nuisance regression, or prewhitening.

prewhiten_info

Vectors of values across locations: phi (AR coefficients averaged across sessions), sigma_sq (residual variance averaged across sessions), and AIC (the maximum across sessions).

call

match.call() for this function call.

INLA Requirement

This function requires the INLA package, which is not a CRAN package. See https://www.r-inla.org/download-install for easy installation instructions.


Create FEM matrices

Description

Create FEM matrices

Usage

galerkin_db(FV, P, surface = FALSE)

Arguments

FV

Matrix of faces in triangularization

P

Matrix of vertex locations in triangularization

surface

(logical) Will this create the SPDE matrices for a surface or not?

Value

A list of matrices C and G appearing in sparse SPDE precision


Get number of locations for various masks

Description

Get number of locations for various masks: total, model, data.

Usage

get_nV(spatial)

Arguments

spatial

spatial

Value

A list of two: T for the total number of locations, and D for the number of data locations. If spatial is provided for voxel data, there is also DB for the number of data locations plus the number of boundary locations.


Extracts posterior density estimates for hyperparameters

Description

Extracts posterior density estimates for hyperparameters

Usage

get_posterior_densities(INLA_model_obj, spde, field_names)

Arguments

INLA_model_obj

An object of class "inla", a result of a call to inla()

spde

The model used for the latent fields in the inla() call, an object of class "inla.spde"

field_names

Descriptive names of model regressors (fields).

Value

Long-form data frame containing posterior densities for the hyperparameters associated with each latent field

INLA Requirement

This function requires the INLA package, which is not a CRAN package. See https://www.r-inla.org/download-install for easy installation instructions.


Extracts posterior density estimates for hyperparameters for volumetric SPDE

Description

Extracts posterior density estimates for hyperparameters for volumetric SPDE

Usage

get_posterior_densities2(INLA_model_obj, field_names)

Arguments

INLA_model_obj

An object of class "inla", a result of a call to inla()

field_names

Descriptive names of model regressors (fields).

Value

Long-form data frame containing posterior densities for the hyperparameters associated with each latent field

INLA Requirement

This function requires the INLA package, which is not a CRAN package. See https://www.r-inla.org/download-install for easy installation instructions.


hpf

Description

hpf

Arguments

hpf

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.


The fix point function for the initialization of kappa2 and phi

Description

The fix point function for the initialization of kappa2 and phi

Usage

init_fixpt(theta, spde, beta_hat)

Arguments

theta

a vector c(kappa2,phi)

spde

an spde object

beta_hat

vector

Value

scalar


Objective function for the initialization of kappa2 and phi

Description

Objective function for the initialization of kappa2 and phi

Usage

init_objfn(theta, spde, beta_hat)

Arguments

theta

a vector c(kappa2,phi)

spde

an spde object

beta_hat

vector

Value

scalar


Intersection mask for BayesGLM or activations result

Description

Intersection mask for BayesGLM or activations result

Usage

intersect_mask(x)

Arguments

x

The list of "fit_bglm", "BGLM", or "act_BGLM" objects.

Value

The intersections masks


Is a matrix or data.frame?

Description

Is this a matrix or data.frame?

Usage

is_matrix_or_df(x)

Arguments

x

The object

Value

Length-one logical.


Function to optimize over kappa2

Description

Function to optimize over kappa2

Usage

kappa_init_fn(kappa2, phi, spde, beta_hat)

Arguments

kappa2

scalar

phi

scalar

spde

an spde object

beta_hat

vector

Value

a scalar


Make log_kappa and log_tau

Description

Make log_kappa and log_tau

Usage

log_kappa_tau(spatial, hyperpriors, verbose)

Arguments

spatial, hyperpriors, verbose

See fit_bayesglm


Make A matrix

Description

Make A matrix

Usage

make_A_mat(spatial)

Arguments

spatial

See BayesGLM

Value

The A matrix


Make A matrix with resampling framework

Description

Make the A matrix for downsampling surface mesh data to a lower resolution.

Usage

make_A_mat_rs(surf, surf_rs)

Arguments

surf

The full-resolution "surf" object.

surf_rs

The downsampled "surf" object.

Value

The A matrix


Make the full SPDE precision based on theta, the spde, and the number of sessions

Description

Make the full SPDE precision based on theta, the spde, and the number of sessions

Usage

make_Q(theta, spde, n_sess)

Arguments

theta

the hyperparameter theta vector of length K * 2 + 1, where the first K elements are the kappas, the next K elements are the phis, and the last element (unused here) corresponds to sigma2

spde

a list containing three spde elements: Cmat, Gmat, and GtCinvG

n_sess

The integer number of sessions

Value

An SPDE prior matrix


Make data list for estimate_model

Description

Make data list to be passed to estimate_model

Usage

make_data_list(y, X, betas, repls)

Arguments

y

Vectorized BOLD data (all voxels, sessions, etc.)

X

List (length = number of sessions) of sparse design matrices size TVxVK from each session, each created using sparse_and_PW

betas

List (length = number of fields) of beta objects from make_replicates

repls

List (length = number of fields) of repl objects from make_replicates

Value

List


Make Mesh

Description

Make INLA triangular mesh from faces and vertices

Usage

make_mesh(vertices, faces)

Arguments

vertices

A V \times 3 matrix, where each row contains the Euclidean coordinates at which a given vertex in the mesh is located. V is the number of vertices in the mesh

faces

An F \times 3 matrix, where each row contains the vertex indices for a given triangular face in the mesh. F is the number of faces in the mesh.

Value

INLA triangular mesh

INLA Requirement

This function requires the INLA package, which is not a CRAN package. See https://www.r-inla.org/download-install for easy installation instructions.


Make replicates

Description

beta and repl vectors are of length nMesh \times nSess \times n_field. The ith repl vector is an indicator vector for the cells corresponding to the ith column of x. The ith beta vector contains data indices (e.g. 1,...,V) in the cells corresponding to the ith column of x.

Usage

make_replicates(nSess, field_names, spatial)

Arguments

nSess

The number of sessions sharing hyperparameters (can be different fields)

field_names

Vector of names for each field

spatial

Spatial info

Value

replicates vector and betas for sessions


Make sqrtInv_all

Description

Make sqrtInv_all for prewhitening

Usage

make_sqrtInv_all(nT, nV, do_pw, n_threads, ar_order, AR_coefs_avg, var_avg)

Arguments

nT, nV, do_pw, n_threads, ar_order, AR_coefs_avg, var_avg

See GLM_est_resid_var_pw.

Value

sqrtInv_all


mask: vertices

Description

mask: vertices

Arguments

mask

A length V logical vector indicating if each vertex is within the input mask.


max_threads

Description

max_threads

Arguments

max_threads

The maximum number of threads to use in the inla-program for model estimation. 0 (default) will use the maximum number of threads allowed by the system.


mean and variance tolerance

Description

mean and variance tolerance

Arguments

meanTol, varTol

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.


mesh: either

Description

mesh: either

Arguments

mesh

An "inla.mesh" object (see make_mesh for surface data)


mesh: INLA only

Description

mesh: INLA only

Arguments

mesh

An "inla.mesh" object (see make_mesh for surface data).


n_threads

Description

n_threads

Arguments

n_threads

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.


nbhd_order

Description

nbhd_order

Arguments

nbhd_order

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.


The negative of the objective function for kappa

Description

The negative of the objective function for kappa

Usage

neg_kappa_fn(kappa2, spde, phi, Sigma, mu)

Arguments

kappa2

scalar

spde

an spde object

phi

scalar

Sigma

dgCMatrix

mu

dgeMatrix

Value

a scalar


The negative of the objective function for kappa without Sig_inv

Description

The negative of the objective function for kappa without Sig_inv

Usage

neg_kappa_fn2(kappa2, spde, phi, P, mu, Vh)

Arguments

kappa2

scalar

spde

an spde object

phi

scalar

P

Matrix of dimension nk by N_s found by solve(Sig_inv,Vh)

mu

dgeMatrix

Vh

random matrix of -1 and 1 of dim dim(P)

Value

a scalar


Streamlined negative objective function for kappa2 using precompiled values

Description

Streamlined negative objective function for kappa2 using precompiled values

Usage

neg_kappa_fn3(kappa2, spde, a_star, b_star, n_sess)

Arguments

kappa2

scalar

spde

an spde object

a_star

precomputed coefficient (scalar)

b_star

precomputed coefficient (scalar)

n_sess

number of sessions (scalar)

Value

scalar output of the negative objective function


Streamlined negative objective function for kappa2 using precompiled values

Description

Streamlined negative objective function for kappa2 using precompiled values

Usage

neg_kappa_fn4(kappa2, spde, a_star, b_star, n_sess)

Arguments

kappa2

scalar

spde

a list of SPDE prior matrices from

a_star

precomputed coefficient (scalar)

b_star

precomputed coefficient (scalar)

n_sess

number of sessions (scalar)

Value

scalar output of the negative objective function


nuisance

Description

nuisance

Arguments

nuisance

(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.


S3 method: use view_xifti to plot a "BGLM" object

Description

S3 method: use view_xifti to plot a "BGLM" object

Usage

## S3 method for class 'BGLM'
plot(
  x,
  Bayes = NULL,
  idx = NULL,
  title = NULL,
  session = NULL,
  zlim = c(-1, 1),
  ...
)

Arguments

x

An object of class "BGLM"

Bayes

TRUE for plotting Bayesian results, FALSE for plotting classical GLM results. Default: NULL, which will use the Bayesian results if available and the classical results if not.

idx

Which field should be plotted? Give the numeric indices or the names. NULL (default) will show all fields. This argument overrides the idx argument to view_xifti.

title

If NULL, the field names associated with idx will be used.

session

Which session should be plotted? NULL (default) will use the first.

zlim

Overrides the zlim argument for view_xifti. Default: c(-1, 1).

...

Additional arguments to view_xifti

Value

Result of the call to ciftiTools::view_cifti.


S3 method: use view_xifti to plot a "BGLM2" object

Description

S3 method: use view_xifti to plot a "BGLM2" object

Usage

## S3 method for class 'BGLM2'
plot(x, idx = NULL, stat = c("contrasts", "activations"), zlim = c(-1, 1), ...)

Arguments

x

An object of class "BGLM2"

idx

Which contrast should be plotted? Give the numeric indices or the names. NULL (default) will show all contrasts. This argument overrides the idx argument to view_xifti.

stat

Estimates of the "contrasts" (default), or their thresholded "activations".

zlim

Overrides the zlim argument for view_xifti. Default: c(-1, 1).

...

Additional arguments to view_xifti

Value

Result of the call to ciftiTools::view_cifti.


S3 method: use view_xifti to plot a "act_BGLM" object

Description

S3 method: use view_xifti to plot a "act_BGLM" object

Usage

## S3 method for class 'act_BGLM'
plot(x, idx = NULL, title = NULL, session = NULL, ...)

Arguments

x

An object of class "act_BGLM"

idx

Which field should be plotted? Give the numeric indices or the names. NULL (default) will show all fields. This argument overrides the idx argument to view_xifti.

title

If NULL, the field names associated with idx will be used.

session

Which session should be plotted? NULL (default) will use the first.

...

Additional arguments to view_xifti

Value

Result of the call to ciftiTools::view_cifti_surface.


S3 method: use view_xifti to plot a "prev_BGLM" object

Description

S3 method: use view_xifti to plot a "prev_BGLM" object

Usage

## S3 method for class 'prev_BGLM'
plot(
  x,
  idx = NULL,
  session = NULL,
  drop_zeros = NULL,
  colors = "plasma",
  zlim = c(0, 1),
  ...
)

Arguments

x

An object of class "prev_BGLM"

idx

Which task should be plotted? Give the numeric indices or the names. NULL (default) will show all tasks. This argument overrides the idx argument to view_xifti.

session

Which session should be plotted? NULL (default) will use the first.

drop_zeros

Color locations without any activation across all results (zero prevalence) the same color as the medial wall? Default: NULL to drop the zeros if only one idx is being plotted.

colors, zlim

See view_xifti.

...

Additional arguments to view_xifti

Value

Result of the call to ciftiTools::view_cifti_surface.


Find values for coefficients used in objective function for kappa2

Description

Find values for coefficients used in objective function for kappa2

Usage

prep_kappa2_optim(spde, mu, phi, P, vh)

Arguments

spde

an spde object

mu

the mean

phi

scale parameter

P

Matrix of dimension nk by N_s found by solve(Sig_inv,Vh)

vh

random matrix of -1 and 1 of dim dim(P)

Value

a list with two elements, which are precomputed coefficients for the optimization function


Activations prevalence.

Description

Activations prevalence.

Usage

prevalence(
  act_list,
  gamma_idx = 1,
  p_test = NULL,
  alpha = 0.05,
  correction = c("FWER", "FDR", "none")
)

Arguments

act_list

List of activations from activations. All should have the same sessions, fields, and brainstructures.

gamma_idx

If activations at multiple thresholds were computed, which threshold should be used for prevalence? Default: the first (lowest).

p_test

For inference: the expected baseline rate of activation for all data locations, under the null hypothesis. Default: NULL (do not perform hypothesis testing).

alpha

Significance level for inference. Default: .05.

correction

For the classical method only: Type of multiple comparisons correction: "FWER" (Bonferroni correction, the default), "FDR" (Benjamini Hochberg), or "none".

Value

A list containing the prevalences of activation, as a proportion of the results from act_list.


Estimate residual autocorrelation for prewhitening

Description

Estimate residual autocorrelation for prewhitening

Usage

pw_estimate(resids, ar_order, aic = FALSE)

Arguments

resids

Estimated residuals in T \times V numeric matrix

ar_order, aic

Order of the AR model used to prewhiten the data at each location. If !aic (default), the order will be exactly ar_order. If aic, the order will be between zero and ar_order, as determined by the AIC.

Value

Estimated AR coefficients and residual variance at every vertex


Smooth AR coefficients and white noise variance

Description

Smooth AR coefficients and white noise variance

Usage

pw_smooth(spatial, AR, var, FWHM = 5)

Arguments

spatial

See fit_bayesglm internal code.

AR

A Vxp matrix of estimated AR coefficients, where V is the number of vertices and p is the AR model order

var

A vector length V containing the white noise variance estimates from the AR model

FWHM

FWHM parameter for smoothing. Remember that \sigma = \frac{FWHM}{2*sqrt(2*log(2)}. Set to 0 or NULL to not do any smoothing. Default: 5.

Value

Smoothed AR coefficients and residual variance at every vertex


Sample from a multivariate normal with mean and precision

Description

Sample from a multivariate normal with mean and precision

Usage

qsample(n, mu, Q)

Arguments

n

number of samples

mu

mean vector (length = p)

Q

sparse p x p positive definite precision matrix (class = dgCMatrix)

Value

An n x p matrix of samples


resamp_res

Description

resamp_res

Arguments

resamp_res

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.


Retroactively mask activations

Description

Retroactively mask activations

Usage

retro_mask_act(x, Masks)

Arguments

x

The activations object

Masks

The masks to be applied to each brain structure of x.

Value

The masked result


Retroactively mask locations from fit_bglm result.

Description

Retroactively mask locations from fit_bglm result.

Usage

retro_mask_fit_bglm(x, mask)

Arguments

x

The "fit_bglm" result

mask

The mask to be applied to x. It's relative to the full mesh/voxel array, not x.

Value

The masked result


Retroactively mask locations from mesh.

Description

Work in progress.

Usage

retro_mask_mesh(x, mask)

Arguments

x

The mesh

mask

The mask to be applied to x (on top of any masks already applied to it.)

Value

The masked result


return_INLA

Description

return_INLA

Arguments

return_INLA

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.


Sequential 2-means variable selection

Description

Sequential 2-means variable selection

Usage

s2m(x, b)

Arguments

x

A vector consisting of all variables of interest for a single draw from a posterior distribution

b

A scale parameter used to determine at what distance cluster centers are considered to be the same.

Value

The number of nonzero values detected within x


Sequential 2-means on array B

Description

Sequential 2-means on array B

Usage

s2m_B(B, sigma)

Arguments

B

An array of posterior samples (typically a matrix), in which the last margin corresponds to a single posterior sample

sigma

A scale parameter used to determine at what distance cluster centers are considered to be the same.

Value

An array of dimension head(dim(B),-1) with a point estimate of B based on the sequential 2-means method


Scale the BOLD timeseries

Description

Scale the BOLD timeseries

Usage

scale_BOLD(BOLD, scale = c("mean", "sd", "none"), v_means = NULL)

Arguments

BOLD

fMRI data as a locations by time (V \times T) numeric matrix.

scale

Option for scaling the BOLD response.

v_means

Original means of the BOLD data. ONLY provide if data has already been centered.

\code{"mean"} scaling will scale the data to percent local signal change.

\code{"sd"} scaling will scale the data by local standard deviation.

\code{"none"} will only center the data, not scale it.

Value

Scale to units of percent local signal change and centers


scale_BOLD

Description

scale_BOLD

Arguments

scale_BOLD

Controls scaling the BOLD response at each location.

"mean":

Scale the data to percent local signal change.

"sd":

Scale the data by local standard deviation.

"none":

Center the data but do not scale it.


Scale the design matrix

Description

Scale the design matrix

Usage

scale_design_mat(design_mat)

Arguments

design_mat

The original (unscaled) design matrix that is T x K, where T is the number of time points, and k is the number of field covariates

Value

A scaled design matrix


scrub

Description

scrub

Arguments

scrub

(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.


seed

Description

seed

Arguments

seed

Random seed (optional). Default: NULL.


session_names

Description

session_names

Arguments

session_names

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.


Organize data for Bayesian GLM

Description

Transforms the usual TxV BOLD data matrix Y into vector form, and the usual TxK design matrix X into big sparse matrix form for use in Bayesian GLM.

Usage

sparse_and_PW(
  BOLD,
  design,
  spatial,
  spde,
  field_names,
  design_type,
  valid_cols,
  nT,
  sqrtInv_all
)

Arguments

BOLD, design, spatial, spde

See fit_bayesglm.

field_names, design_type

See fit_bayesglm.

valid_cols, nT, sqrtInv_all

See fit_bayesglm.

Details

The Bayesian GLM requires y (a vector of length TV containing the BOLD data) and X_k (a sparse TVxV matrix corresponding to the kth field regressor) for each field k. The design matrices are combined as A=cbind(X_1,...,X_K).

The Bayesian GLM requires y (a vector of length TV containing the BOLD data) and X_k (a sparse TVxV matrix corresponding to the kth field regressor) for each field k. The design matrices are combined as A=cbind(X_1,...,X_K).

Value

A list containing fields y and A (see Details)


Calculate the SPDE covariance

Description

Calculate the SPDE covariance

Usage

spde_Q_phi(kappa2, phi, spde)

Arguments

kappa2

A scalar

phi

A scalar

spde

An object containing the information about the mesh structure for the SPDE prior

Value

The SPDE prior matrix


Summarize a "BGLM" object

Description

Summary method for class "BGLM"

Usage

## S3 method for class 'BGLM'
summary(object, ...)

## S3 method for class 'summary.BGLM'
print(x, ...)

## S3 method for class 'BGLM'
print(x, ...)

Arguments

object

Object of class "BGLM".

...

further arguments passed to or from other methods.

x

Object of class "summary.BGLM".

Value

A "summary.BGLM" object, a list summarizing the properties of object.

NULL, invisibly.

NULL, invisibly.


Summarize a "BGLM2" object

Description

Summary method for class "BGLM2"

Usage

## S3 method for class 'BGLM2'
summary(object, ...)

## S3 method for class 'summary.BGLM2'
print(x, ...)

## S3 method for class 'BGLM2'
print(x, ...)

Arguments

object

Object of class "BGLM2".

...

further arguments passed to or from other methods.

x

Object of class "summary.BGLM2".

Value

A "summary.BGLM2" object, a list summarizing the properties of object.

NULL, invisibly.

NULL, invisibly.


Summarize a "act_BGLM" object

Description

Summary method for class "act_BGLM"

Usage

## S3 method for class 'act_BGLM'
summary(object, ...)

## S3 method for class 'summary.act_BGLM'
print(x, ...)

## S3 method for class 'act_BGLM'
print(x, ...)

Arguments

object

Object of class "act_BGLM".

...

further arguments passed to or from other methods.

x

Object of class "summary.act_BGLM".

Value

A "summary.act_BGLM" object, a list summarizing the properties of object.

NULL, invisibly.

NULL, invisibly.


Summarize a "act_fit_bglm" object

Description

Summary method for class "act_fit_bglm"

Usage

## S3 method for class 'act_fit_bglm'
summary(object, ...)

## S3 method for class 'summary.act_fit_bglm'
print(x, ...)

## S3 method for class 'act_fit_bglm'
print(x, ...)

Arguments

object

Object of class "act_fit_bglm".

...

further arguments passed to or from other methods.

x

Object of class "summary.act_fit_bglm".

Value

A "summary.act_fit_bglm" object, a list summarizing the properties of object.

NULL, invisibly.

NULL, invisibly.


Summarize a "fit_bglm" object

Description

Summary method for class "fit_bglm"

Usage

## S3 method for class 'fit_bglm'
summary(object, ...)

## S3 method for class 'summary.fit_bglm'
print(x, ...)

## S3 method for class 'fit_bglm'
print(x, ...)

Arguments

object

Object of class "fit_bglm".

...

further arguments passed to or from other methods.

x

Object of class "summary.fit_bglm".

Value

A "summary.fit_bglm" object, a list summarizing the properties of object.

NULL, invisibly.

NULL, invisibly.


Summarize a "fit_bglm2" object

Description

Summary method for class "fit_bglm2"

Usage

## S3 method for class 'fit_bglm2'
summary(object, ...)

## S3 method for class 'summary.fit_bglm2'
print(x, ...)

## S3 method for class 'fit_bglm2'
print(x, ...)

Arguments

object

Object of class "fit_bglm2".

...

further arguments passed to or from other methods.

x

Object of class "summary.fit_bglm2".

Value

A "summary.fit_bglm2" object, a list summarizing the properties of object.

NULL, invisibly.

NULL, invisibly.


Summarize a "prev_BGLM" object

Description

Summary method for class "prev_BGLM"

Summary method for class "prev_BGLM"

Usage

## S3 method for class 'prev_BGLM'
summary(object, ...)

## S3 method for class 'summary.prev_BGLM'
print(x, ...)

## S3 method for class 'prev_BGLM'
print(x, ...)

## S3 method for class 'prev_BGLM'
summary(object, ...)

## S3 method for class 'summary.prev_BGLM'
print(x, ...)

## S3 method for class 'prev_BGLM'
print(x, ...)

Arguments

object

Object of class "prev_BGLM".

...

further arguments passed to or from other methods.

x

Object of class "summary.prev_BGLM".

Value

A "summary.prev_BGLM" object, a list summarizing the properties of object.

NULL, invisibly.

NULL, invisibly.

A "summary.prev_BGLM" object, a list summarizing the properties of object.

NULL, invisibly.

NULL, invisibly.


Summarize a "prev_fit_bglm" object

Description

Summary method for class "prev_fit_bglm"

Summary method for class "prev_fit_bglm"

Usage

## S3 method for class 'prev_fit_bglm'
summary(object, ...)

## S3 method for class 'summary.prev_fit_bglm'
print(x, ...)

## S3 method for class 'prev_fit_bglm'
print(x, ...)

## S3 method for class 'prev_fit_bglm'
summary(object, ...)

## S3 method for class 'summary.prev_fit_bglm'
print(x, ...)

## S3 method for class 'prev_fit_bglm'
print(x, ...)

Arguments

object

Object of class "prev_fit_bglm".

...

further arguments passed to or from other methods.

x

Object of class "summary.prev_fit_bglm".

Value

A "summary.prev_fit_bglm" object, a list summarizing the properties of object.

NULL, invisibly.

NULL, invisibly.

A "summary.prev_fit_bglm" object, a list summarizing the properties of object.

NULL, invisibly.

NULL, invisibly.


surfaces

Description

surfaces

Arguments

surfL, surfR

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.


trim_INLA

Description

trim_INLA

Arguments

trim_INLA

(logical) should the INLA_model_obj within the result be trimmed to only what is necessary to use activations? Default: TRUE.


Trim INLA object

Description

Trim an INLA object to only include what is necessary for activations or BayesGLM2.

Usage

trim_INLA_model_obj(INLA_model_obj, minimal = FALSE)

Arguments

INLA_model_obj

An object of class "inla".

minimal

Just keep the two parameters needed for BayesGLM2? Default: FALSE. !minimal is required for activations, but minimal is sufficient for BayesGLM2.

Value

A trimmed "inla" object.


Unmask data

Description

maskMdat is a subset of maskIn. x is masked by maskMdat; this function makes it masked by maskIn instead and uses fill_val for the previously masked-out locations.

Usage

unmask_Mdat2In(x, maskIn, maskMdat, fill_val = NA)

Arguments

x

The matrix to unmask

maskIn, maskMdat

The input mask, and the modeled locations mask, as logical vectors. The latter should be a subset of the former. They should have the same length.

fill_val

The fill value. Default: NA.


Validate spatial

Description

Validate spatial

Usage

validate_spatial(spatial)

Arguments

spatial

spatial

Value

NULL, invisibly


verbose

Description

verbose

Arguments

verbose

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.


Surface area of each vertex

Description

Compute surface areas of each vertex in a triangular mesh.

Usage

vertex_areas(mesh)

Arguments

mesh

An "inla.mesh" object (see make_mesh for surface data).

Value

Vector of areas

INLA Requirement

This function requires the INLA package, which is not a CRAN package. See https://www.r-inla.org/download-install for easy installation instructions.


vertices

Description

vertices

Arguments

vertices

A V \times 3 matrix, where each row contains the Euclidean coordinates at which a given vertex in the mesh is located. V is the number of vertices in the mesh


Construct a triangular mesh from a 3D volumetric mask

Description

Construct a triangular mesh from a 3D volumetric mask

Usage

vol2spde(mask, res, nbhd_order = 1, buffer = c(1, 1, 3, 4, 4))

Arguments

mask

An array of 0s and 1s representing a volumetric mask

res

The spatial resolution in each direction, in mm. For example, c(2,2,2) indicates 2mm isotropic voxels.

nbhd_order

For volumetric data, what order neighborhood around data locations to keep? (0 = no neighbors, 1 = 1st-order neighbors, 2 = 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.

buffer

For volumetric data, size of extra voxels layers around the bounding box, in terms of voxels. Set to NULL for no buffer.

Value

An inla.spde2 object.

INLA Requirement

This function requires the INLA package, which is not a CRAN package. See https://www.r-inla.org/download-install for easy installation instructions.