Type:	Package
Title:	Bayesian Hierarchical Multi-Subject Multiscale Analysis of Functional MRI (fMRI) Data
Version:	2.2
Date:	2024-06-12
Description:	Package BHMSMAfMRI performs Bayesian hierarchical multi-subject multiscale analysis of fMRI data as described in Sanyal & Ferreira (2012) <doi:10.1016/j.neuroimage.2012.08.041>, or other multiscale data, using wavelet-based prior that borrows strength across subjects and provides posterior smoothed images of the effect sizes and samples from the posterior distribution.
License:	GPL-2 | GPL-3 [expanded from: GPL (≥ 2)]
Imports:	Rcpp (≥ 1.0.9), oro.nifti, wavethresh
LinkingTo:	Rcpp, RcppArmadillo
Suggests:	knitr, rmarkdown, bookdown
VignetteBuilder:	knitr
URL:	https://nilotpalsanyal.github.io/BHMSMAfMRI/
BugReports:	https://github.com/nilotpalsanyal/BHMSMAfMRI/issues
Repository:	CRAN
NeedsCompilation:	yes
Packaged:	2024-06-12 08:14:26 UTC; nsanyal
Author:	Nilotpal Sanyal [aut, cre], Marco A.R. Ferreira [aut]
Maintainer:	Nilotpal Sanyal <nilotpal.sanyal@gmail.com>
Date/Publication:	2024-06-12 08:30:02 UTC

Bayesian hierarchical multi-subject multiscale analysis (BHMSMA) of functional MRI data or other multiscale data

Description

BHMSMA is the main function that performs BHMSMA (Sanyal & Ferreira, 2012) of fMRI data, or other multiscale data, using wavelet-based prior that borrows strength across subjects and provides posterior smooth estimates of the effect sizes. Specifically, BHMSMA constructs a 2D regression coefficient map (e.g., corresponding to a single brain slice) of a given regressor and returns its posterior smoothed version based on multi-subject or single subject analyses.

Usage

BHMSMA(n, grid, data, designmat, k, analysis="multi", 
  truecoef=NULL, wave.family="DaubLeAsymm", filter.number=6, 
  bc="periodic")

Arguments

Details

The wavelet computations are performed by using the R package wavethresh.

Value

A list containing the following.

Author(s)

Nilotpal Sanyal, Marco Ferreira

Maintainer: Nilotpal Sanyal <nilotpal.sanyal@gmail.com>

References

Sanyal, Nilotpal, and Ferreira, Marco A.R. (2012). Bayesian hierarchical multi-subject multiscale analysis of functional MRI data. Neuroimage, 63, 3, 1519-1531.

See Also

readfmridata, glmcoef, waveletcoef, hyperparamest, postmixprob, postwaveletcoef, substituteWaveletCoef, postglmcoef, imwd, imwr

Examples

# BHMSMA multi-subject analysis for simulated (fMRI) 
# data at 4 timepoints over an 8x8 grid (of a brain 
# slice) for 3 subjects
set.seed(1)
n <- 3
grid <- 8
ntime <- 4
data <- array(rnorm(n*grid*grid*ntime),
  dim=c(n,grid,grid,ntime))
designmat <- cbind(c(1,1,1,1),c(1,0,1,0))
k <- 2
analysis <- "multi"
BHMSMAmulti <- BHMSMA(n, grid, data, designmat, 
  k, analysis)

zlim = c(0,max(abs(BHMSMAmulti$GLMCoefStandardized)))
par(mfrow=c(1,2))
image( abs(BHMSMAmulti$GLMCoefStandardized[1,,,k]),
  col=heat.colors(12),zlim=zlim,main="GLM coef map")
image( abs(BHMSMAmulti$GLMcoefposterior[1,,]),
  col=heat.colors(12),zlim=zlim,main="GLM coef posterior map")


## Not run: 
# BHMSMA multi-subject analysis for simulated (fMRI) 
# data at 100 timepoints over an 64x64 grid (of a 
# brain slice) for 15 subjects
# (takes ~12s in a 2.8 GHz Quad-Core Intel Core i7 processor)
set.seed(1)
n <- 15
grid <- 64
ntime <- 100
data <- array(rnorm(n*grid*grid*ntime),
          dim=c(n,grid,grid,ntime))
designmat <- cbind(rep(1,ntime),runif(ntime))
k <- 2
analysis <- "multi"
system.time({BHMSMAmulti <- BHMSMA(n,grid,data, 
  designmat,k,analysis)})

## End(Not run)

Bayesian Hierarchical Multi-Subject Multiscale Analysis (BHMSMA) of Functional MRI Data

Description

The BHMSMAfMRI package performs BHMSMA (Sanyal & Ferreira, 2012) of fMRI data, or other multiscale data, using wavelet based prior that borrows strength across subjects and provides posterior smoothed images of the effect sizes and samples from the posterior distribution. The package currently considers analysis of 2D slices/grids only.

Details

Import fMRI data using:
readfmridata

The main analysis function, which provides subject-specific posterior estimates, is:
BHMSMA

The main function sucessively calls the following functions:
glmcoef (get regression coefficients)
waveletcoef (get wavelet coefficients)
hyperparamest (estimate model hyperparameters)
postmixprob (estimate posterior mixture probabilities of wavelet coefficients)
postwaveletcoef (compute posterior estimates of wavelet coefficients)
postglmcoef (compute posterior estimates of regression coefficients)

For posterior group estimates of regression coefficients use:
postgroupglmcoef

For posterior uncertainty estimates use:
postsamples

Internal sample data:
fmridata

Miscellaneous:
substituteWaveletCoef

Author(s)

Nilotpal Sanyal <nilotpal.sanyal@gmail.com>, Marco Ferreira <marf@vt.edu>

Maintainer: Nilotpal Sanyal <nilotpal.sanyal@gmail.com>

References

Sanyal, Nilotpal, and Ferreira, Marco A.R. (2012). Bayesian hierarchical multi-subject multiscale analysis of functional MRI data. Neuroimage, 63, 3, 1519-1531.

A simulated fMRI data for 3 subjects

Description

A simulated fMRI data containing true regression coefficients images for three subjects and design matrix

Usage

data(fmridata)

Format

A list containing the following.

grid

=32. The image dimension is 32 by 32.

nsubject

=3.

TrueCoeff

An array of dimension (3,32,32), containing the true regression coefficients for the 3 subjects.

DesignMatrix

A matrix with 9 columns and 2 rows. The first column is a column of ones.

Details

This dataset contains only the true coefficients. The noisy fMRI data, which are generated by adding Gaussian random noise to these true coefficients, are included in the extdata directory within the package directory. The function readfmridata can be used to read those data files. The true coefficients and the noisy data both are generated using the R package neuRosim. The following specifications were used to generate the data: totaltime=18, onsets=seq(1,18,by=8), durations=1, TR=2, effectsize=1, hrf="double-gamma", regions=3, radius=c(1,1,1), form="sphere", fading=1, SNR=1.5, noise="white". The centers of the activation regions were chosen manually. For information regarding the specifications, see neuRosim help.

Fit GLM (general linear model) to the fMRI time-series of all voxels within a single 2D brain slice

Description

glmcoef fits a GLM to the fMRI time-series of all voxels within a single 2D brain slice for each subject, and returns standardized GLM coefficients along with their standard error for the included regressors (it does not add any intercept by itself).

Usage

glmcoef(n, grid, data, designmat)

Arguments

Value

A list containing the following.

Author(s)

Nilotpal Sanyal, Marco Ferreira

Maintainer: Nilotpal Sanyal <nilotpal.sanyal@gmail.com>

References

Friston, K.J., Holmes, A.P., Worsley, K.J., Poline, J., Frith, C.D., Frackowiak, R.S.J., 1994. Statistical parametric maps in functional imaging: a general linear approach. Hum. Brain Mapp. 2 (4), 189-210.

See Also

lm, waveletcoef

Examples

set.seed(1)
n <- 3
grid <- 8
ntime <- 10
designmat <- cbind(rep(1,10),c(rep(c(1,0),5)))
data <- array(dim=c(n,grid,grid,ntime),
  rnorm(n*grid*grid*ntime))
glm.fit <- glmcoef(n,grid,data,designmat)
dim(glm.fit$GLMCoefStandardized)
#[1] 3 8 8

Obtain estimates of the hyperparameters of the BHMSMA model

Description

hyperparamest computes the MLEs (maximum likelihood estimates) of the hyperparameters of the BHMSMA model using an empirical Bayes approach for multi-subject or single subject analyses, and returns the hyperparameters estimates along with their covariance matrix estimate (see References).

Usage

hyperparamest(n, grid, waveletcoefmat, analysis)

Arguments

Value

A list containing the following.

Author(s)

Nilotpal Sanyal, Marco Ferreira

Maintainer: Nilotpal Sanyal <nilotpal.sanyal@gmail.com>

References

Sanyal, Nilotpal, and Ferreira, Marco A.R. (2012). Bayesian hierarchical multi-subject multiscale analysis of functional MRI data. Neuroimage, 63, 3, 1519-1531.

See Also

waveletcoef, nlminb, postmixprob

Examples

set.seed(1)
n <- 3
grid <- 8
waveletcoefmat <- array(dim=c(n,grid^2-1),
  rnorm(n*(grid^2-1)))
analysis <- "multi"
hyperest <- hyperparamest(n,grid,waveletcoefmat,analysis)
hyperest$hyperparam
# [1]  1.00000  1.00000  1.00000  1.00000  0.00000 28.37678

Obtain posterior estimate of a 2D GLM coefficients map of a regressor

Description

postglmcoef computes posterior mean (or median) of a 2D GLM coefficients map (e.g., corresponding to a single brain slice) of a regressor using the posterior mean (or median) of the corresponding wavelet coefficients in the inverse discrete wavelet transform for each subject based on multi-subject or single subject analyses (see References).

Usage

postglmcoef(n, grid, glmcoefstd, postmeanwaveletcoef, 
wave.family="DaubLeAsymm", filter.number=6, bc="periodic")

Arguments

Details

The wavelet transformation and reconstruction are performed by using the functions imwd and imwr, respectively.

Value

A list containing the following.

Author(s)

Nilotpal Sanyal, Marco Ferreira

Maintainer: Nilotpal Sanyal <nilotpal.sanyal@gmail.com>

See Also

glmcoef, postwaveletcoef, substituteWaveletCoef, imwr, postgroupglmcoef, postsamples

Examples

set.seed(1)
n <- 3
grid <- 8
glmcoefstd <- array(rnorm(n*grid*grid),
  dim=c(n,grid,grid))
postmeanwaveletcoef <- array(rnorm(n*(grid^2-1)),
  dim=c(n,(grid^2-1)))
postmeanglmcoef <- postglmcoef(n,grid,glmcoefstd,
  postmeanwaveletcoef)
dim(postmeanglmcoef$GLMcoefposterior)
#[1] 3 8 8

Obtain posterior group estimate of a 2D GLM coefficients map of a regressor

Description

postgroupglmcoef computes posterior group mean (or group median) of a 2D GLM coefficients map (e.g., corresponding to a single brain slice) of a regressor using the posterior means (or medians) of the corresponding wavelet coefficients from all subjects in the inverse discrete wavelet transform based on multi-subject or single subject analyses (see References).

Usage

postgroupglmcoef( n, grid, glmcoefstd, postmeanwaveletcoef, 
wave.family="DaubLeAsymm", filter.number=6, bc="periodic" )

Arguments

Details

The wavelet transformation and reconstruction are performed by using the functions imwd and imwr, respectively.

Value

A list containing the following.

Author(s)

Nilotpal Sanyal, Marco Ferreira

Maintainer: Nilotpal Sanyal <nilotpal.sanyal@gmail.com>

References

Sanyal, Nilotpal, and Ferreira, Marco A.R. (2012). Bayesian hierarchical multi-subject multiscale analysis of functional MRI data. Neuroimage, 63, 3, 1519-1531.

See Also

readfmridata, glmcoef, postglmcoef, postwaveletcoef, substituteWaveletCoef, imwd, imwr

Examples

set.seed(1)
n <- 3
grid <- 8
glmcoefstd <- array(rnorm(n*grid*grid),
  dim=c(n,grid,grid))
postmeanwaveletcoef <- array(rnorm(n*(grid^2-1)),
  dim=c(n,grid^2-1))
post.groupcoef <- postgroupglmcoef(n,grid,glmcoefstd,
  postmeanwaveletcoef)
dim(post.groupcoef$groupcoef)
#[1] 8 8

Obtain estimates of the mixture probabilities defining the BHMSMA posterior wavelet coefficients distributions

Description

postmixprob computes the mixture probabilities (piklj.bar), which define the marginal posterior distribution of the wavelet coefficients of the BHMSMA model, using Newton Cotes algorithm for each subject based on multi-subject or single subject analyses, and returns the same (see References).

Usage

postmixprob(n, grid, waveletcoefmat, hyperparam, analysis)

Arguments

Value

A list containing the following.

Author(s)

Nilotpal Sanyal, Marco Ferreira

Maintainer: Nilotpal Sanyal <nilotpal.sanyal@gmail.com>

References

Sanyal, Nilotpal, and Ferreira, Marco A.R. (2012). Bayesian hierarchical multi-subject multiscale analysis of functional MRI data. Neuroimage, 63, 3, 1519-1531.

See Also

waveletcoef, hyperparamest, postwaveletcoef

Examples

set.seed(1)
n <- 3
grid <- 8
waveletcoefmat <- matrix(nrow=n,ncol=grid^2-1)
for(i in 1:n) waveletcoefmat[i,] <- rnorm(grid^2-1)
hyperparam <- rep(.1,6)
analysis <- "multi"
pkljbar <- postmixprob(n,grid,waveletcoefmat,hyperparam,
  analysis)
dim(pkljbar$pkljbar)
#[1]  3 63

Obtain samples from the posterior distribution of a 2D GLM coefficient map.

Description

postsamples generates samples from the posterior distribution of a 2D GLM coefficient map (e.g., corresponding to a single brain slice) of a regressor in the BHMSMA model for each subject based on multi-subject or single subject analyses (see References).

Usage

postsamples(nsample, n, grid, glmcoefstd, waveletcoefmat, 
hyperparam, pkljbar, analysis, wave.family="DaubLeAsymm", 
filter.number=6, bc="periodic", seed)

Arguments

Details

The wavelet computations are performed by using the R package wavethresh.

Value

A list containing the following.

Author(s)

Nilotpal Sanyal, Marco Ferreira

Maintainer: Nilotpal Sanyal <nilotpal.sanyal@gmail.com>

References

Sanyal, Nilotpal, and Ferreira, Marco A.R. (2012). Bayesian hierarchical multi-subject multiscale analysis of functional MRI data. Neuroimage, 63, 3, 1519-1531.

Morris, J.S. et al. (2011). Automated analysis of quantitative image data using isomorphic functional mixed models, with application to proteomic data. Ann. Appl. Stat. 5, 894-923.

See Also

readfmridata, glmcoef, waveletcoef, hyperparamest, postmixprob, postwaveletcoef, substituteWaveletCoef, postglmcoef, imwd, imwr

Examples

set.seed(1)
n <- 3
grid <- 8
nsample <- 5
glmcoefstd <- array(rnorm(n*grid*grid),
  dim=c(n,grid,grid))
waveletcoefmat <- array(rnorm(n*(grid^2-1)),
  dim=c(n,(grid^2-1)))
hyperparam <- rep(.2,6)
pkljbar <- array(runif(n*(grid^2-1)),
  dim=c(n,(grid^2-1)))
analysis <- "multi"
postsample <- postsamples(nsample,n,grid,glmcoefstd, 
waveletcoefmat, hyperparam,pkljbar,analysis,seed=1)
dim(postsample$samples)
#[1] 3 8 8 5

Obtain posterior estimates of the BHMSMA wavelet coefficients

Description

postwaveletcoef computes posterior mean and posterior median of the wavelet coefficients of the BHMSMA model for each subject based on multi-subject or single subject analyses (see References).

Usage

postwaveletcoef(n, grid, waveletcoefmat, hyperparam, 
pkljbar, analysis)

Arguments

Value

A list containing the following.

Author(s)

Nilotpal Sanyal, Marco Ferreira

Maintainer: Nilotpal Sanyal <nilotpal.sanyal@gmail.com>

References

Sanyal, Nilotpal, and Ferreira, Marco A.R. (2012). Bayesian hierarchical multi-subject multiscale analysis of functional MRI data. Neuroimage, 63, 3, 1519-1531.

See Also

waveletcoef, hyperparamest, postmixprob, postglmcoef

Examples

set.seed(1)
n <- 3
grid <- 8
nsample <- 5
waveletcoefmat <- array(rnorm(n*(grid^2-1)),
  dim=c(n,grid^2-1))
hyperparam <- rep(.2,6)
pkljbar <- array(runif(n*(grid^2-1)),
  dim=c(n,grid^2-1))
analysis <- "multi"
postwavecoef <- postwaveletcoef(n,grid,waveletcoefmat, 
hyperparam,pkljbar,analysis)
dim(postwavecoef$PostMeanWaveletCoef)
#[1]  3 63

Import fMRI data from various fMRI image files

Description

readfmridata reads and imports fMRI data from various fMRI image files (Analyze, NIFTI and AFNI) into a 4D array. It is just a convenient wrapper around the data importing functions provided by the oro.nifti package.

Usage

readfmridata( directory, format, prefix, nimages, dim.image, 
  nii=TRUE )

Arguments

Details

The function uses package oro.nifti for reading from fMRI data files.

Value

A list containing the following:

Author(s)

Nilotpal Sanyal, Marco Ferreira Maintainer: Nilotpal Sanyal <nilotpal.sanyal@gmail.com>

See Also

readANALYZE, readNIfTI, readAFNI, BHMSMA

Examples

# import simmulated fMRI data from image files provided within this package
fpath <- system.file("extdata", package="BHMSMAfMRI")
untar(paste0(fpath,"/fmridata.tar"), exdir=tempdir())
data <- array(dim=c(3,32,32,9))
for(subject in 1:3)
{
  directory <- paste0(tempdir(),"/fmridata","/s0",subject,"/")
  a <- readfmridata(directory, format="Analyze", prefix=paste0("s0",subject,"_t"),
  nimages=9, dim.image=c(32,32,1))
  data[subject,,,] <- a[,,1,]
}
dim(a)

Substitute 2D wavelet transform coefficients with user-given values

Description

substituteWaveletCoef substitutes the wavelet coefficients stored wavelet object generated through 2D wavelet transform with user-given values and returns the modified wavelet object.

Usage

substituteWaveletCoef(grid, waveletobj, values)

Arguments

Details

The maximum value of grid for this package is 512.

Value

A wavelet object of class imwd.object with updated wavelet coefficients.

Author(s)

Nilotpal Sanyal, Marco Ferreira

Maintainer: Nilotpal Sanyal <nilotpal.sanyal@gmail.com>

See Also

imwd

Examples

set.seed(1)
n <- 3
grid <- 8
ntime <- 10
designmat <- cbind( rep(1,10), c(rep(c(1,0),5)) )
data <- array(dim=c(n,grid,grid,ntime),
  rnorm(n*grid*grid*ntime))
glm.fit <- glmcoef(n, grid, data, designmat)
glmcoefstd <- glm.fit$GLMCoefStandardized[,,,1]
dwt = wavethresh::imwd(glmcoefstd[1,,],type="wavelet",
  family="DaubLeAsymm",filter.number=6,bc="periodic")
dwt

values = rnorm(grid^2-1)
dwtnew = substituteWaveletCoef(grid,dwt,values)
dwtnew

Apply discrete wavelet transform (DWT) to a 2D GLM coefficient map of a regressor

Description

waveletcoef applies DWT to a 2D GLM coefficient map (e.g., corresponding to a single brain slice) of a regressor for each subject, and returns the wavelet coefficients at all resolution levels. This function wraps around the wavelet transformation function imwd of the wavethresh package.

Usage

waveletcoef(n, grid, glmcoefstd, wave.family="DaubLeAsymm", 
  filter.number=6, bc="periodic")

Arguments

Details

The wavelet decomposition is performed by using the function imwd.

Value

A list containing the following.

Author(s)

Nilotpal Sanyal, Marco Ferreira

Maintainer: Nilotpal Sanyal <nilotpal.sanyal@gmail.com>

See Also

imwd, hyperparamest

Examples

set.seed(1)
n <- 3
grid <- 8
ntime <- 10
designmat <- cbind( rep(1,10), c(rep(c(1,0),5)) )
data <- array(dim=c(n,grid,grid,ntime),
  rnorm(n*grid*grid*ntime))
glm.fit <- glmcoef(n,grid,data,designmat)
glmcoefstd <- glm.fit$GLMCoefStandardized[,,,1]
wavecoef <- waveletcoef(n,grid,glmcoefstd)
dim(wavecoef$WaveletCoefficientMatrix)
#[1]  3 63