Package 'BHMSMAfMRI' reference manual

Title:	Bayesian Hierarchical Multi-Subject Multiscale Analysis of Functional MRI (fMRI) Data
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.
Authors:	Nilotpal Sanyal [aut, cre] , Marco A.R. Ferreira [aut]
Maintainer:	Nilotpal Sanyal <[email protected]>
License:	GPL (>= 2)
Version:	2.1
Built:	2024-11-20 04:04:50 UTC
Source:	https://github.com/nilotpalsanyal/bhmsmafmri

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

Package:	BHMSMAfMRI
Type:	Package
Version:	2.1
Date:	2022-10-01
License:	GPL (>= 2)

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 <[email protected]>, Marco Ferreira <[email protected]>

Maintainer: Nilotpal Sanyal <[email protected]>

References

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

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")
BHMSMA(n, grid, data, designmat, k, analysis="multi", 
  truecoef=NULL, wave.family="DaubLeAsymm", filter.number=6, 
  bc="periodic")

Arguments

`n`	Number of subjects.
`grid`	The number of voxels in one row (or column) of the brain slice of interest. Must be a power of 2. The total number of voxels is `grid^2`. The maximum value of `grid` for this package is 512.
`data`	The data in the form of an array with dimension `(n,grid,grid,ntime)`, where ntime is the size of the time series for each voxel. The rows (first dimension) must correspond to `n`.
`designmat`	The design matrix used to generate the data. The rows must correspond to `n`. It should include an intercept column if desired. No intercept column will be added by the function. The GLM fit will include all the columns of `designmat`.
`k`	Index of the regressor chosen for analysis, consistently with `designmat`. After fitting GLM, the regression coefficients of this regressor will be considered for the rest of the analysis.
`analysis`	"MSA" or "SSA", depending on whether performing multi-subject analysis or single subject analysis.
`truecoef`	If available, the true GLM coefficients in the form of an array with dimension `(n,grid,grid)`. NULL by default.
`wave.family`	The family of wavelets to use - "DaubExPhase" or "DaubLeAsymm". Default is "DaubLeAsymm".
`filter.number`	The number of vanishing moments of the wavelet. Default is 6.
`bc`	The boundary condition to use - "periodic" or "symmetric". Default is "periodic".

Details

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

Value

A list containing the following.

`GLMCoefStandardized`	An array of dimension `(n,grid,grid)`, containing for each subject the standardized GLM coefficients obtained by fitting GLM to the time-series corresponding to the voxels.
`GLMCoefSE`	An array of dimension `(n,grid,grid)`, containing for each subject the estimated standard errors of the GLM coefficients.
`WaveletCoefficientMatrix`	A matrix of dimension `(n,grid^2-1)`, containing for each subject the wavelet coefficients of all levels stacked together (by the increasing order of resolution level).
`hyperparam`	A vector containing the estimates of the six hyperparameters.
`hyperparamVar`	Estimated covariance matrix of the hyperparameters.
`posteriorMixProb`	A matrix of dimension `(n,grid^2-1)`, containing the piklj.bar values (see Reference for details).
`PostMeanWaveletCoef`	A matrix of size `(n,grid^2-1)`, containing for each subject the posterior mean of the wavelet coefficients of all levels stacked together (by the increasing order of resolution level).
`GLMcoefposterior`	An array of dimension `(n,grid,grid)`, containing for each subject the posterior means of the standardized GLM coefficients.
`MSE`	MSE of the posterior estimates of the GLM coefficients, if the true values of the GLM coefficients are available.

Author(s)

Nilotpal Sanyal, Marco Ferreira

Maintainer: Nilotpal Sanyal <[email protected]>

References

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

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

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)
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)
glmcoef(n, grid, data, designmat)

Arguments

`n`	Number of subjects.
`grid`	The number of voxels in one row (or, one column) of the brain slice of interest. Must be a power of 2. The total number of voxels is `grid^2`. The maximum value of `grid` for this package is 512.
`data`	The data in the form of an array with dimension `(n,grid,grid,ntime)`, where ntime is the size of the time series for each voxel.
`designmat`	The design matrix used to generate the data. An intercept column should be included unless not desired.

Value

A list containing the following.

`GLMCoefStandardized`	An array of dimension `(n,grid,grid)`, containing for each subject the standardized GLM coefficients obtained by fitting GLM to the time-series corresponding to the voxels.
`GLMCoefSE`	An array of dimension `(n,grid,grid)`, containing for each subject the estimated standard errors of the GLM coefficients.

Author(s)

Nilotpal Sanyal, Marco Ferreira

Maintainer: Nilotpal Sanyal <[email protected]>

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.

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
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)
hyperparamest(n, grid, waveletcoefmat, analysis)

Arguments

`n`	Number of subjects.
`grid`	The number of voxels in one row (or, one column) of the brain slice of interest. Must be a power of 2. The total number of voxels is `grid^2`. The maximum value of `grid` for this package is 512.
`waveletcoefmat`	A matrix of dimension `(n,grid^2-1)`, containing for each subject the wavelet coefficients of all levels stacked together (by the increasing order of resolution level).
`analysis`	"multi" or "single", depending on whether performing multi-subject analysis or single subject analysis.

Value

A list containing the following.

`hyperparam`	A vector containing the estimates of the six hyperparameters of the BHMSMA model.
`hyperparamVar`	Estimated covariance matrix of the hyperparameters.

Author(s)

Nilotpal Sanyal, Marco Ferreira

Maintainer: Nilotpal Sanyal <[email protected]>

References

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

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
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")
postglmcoef(n, grid, glmcoefstd, postmeanwaveletcoef, 
wave.family="DaubLeAsymm", filter.number=6, bc="periodic")

Arguments

`n`	Number of subjects.
`grid`	The number of voxels in one row (or, one column) of the brain slice of interest. Must be a power of 2. The total number of voxels is `grid^2`. The maximum value of `grid` for this package is 512.
`glmcoefstd`	An array of dimension `(n,grid,grid)`, containing for each subject the standardized GLM coefficients obtained by fitting GLM to the time-series corresponding to the voxels.
`postmeanwaveletcoef`	A matrix of size `(n,grid^2-1)`, containing for each subject the posterior mean of the wavelet coefficients of all levels stacked together (by the increasing order of resolution level).
`wave.family`	The family of wavelets to use - "DaubExPhase" or "DaubLeAsymm". Default is "DaubLeAsymm".
`filter.number`	The number of vanishing moments of the wavelet. Default is 6.
`bc`	The boundary condition to use - "periodic" or "symmetric". Default is "periodic".

Details

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

Value

A list containing the following.

GLMcoefposterior

An array of dimension (n,grid,grid), containing for each subject the posterior means of the standardized GLM coefficients.

Author(s)

Nilotpal Sanyal, Marco Ferreira

Maintainer: Nilotpal Sanyal <[email protected]>

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
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" )
postgroupglmcoef( n, grid, glmcoefstd, postmeanwaveletcoef, 
wave.family="DaubLeAsymm", filter.number=6, bc="periodic" )

Arguments

`n`	Number of subjects.
`grid`	The number of voxels in one row (or, one column) of the brain slice of interest. Must be a power of 2. The total number of voxels is `grid^2`. The maximum value of `grid` for this package is 512.
`glmcoefstd`	An array of dimension `(n,grid,grid)`, containing for each subject the standardized GLM coefficients obtained by fitting GLM to the time-series corresponding to the voxels.
`postmeanwaveletcoef`	A matrix of size `(n,grid^2-1)`, containing for each subject the posterior mean of the wavelet coefficients of all levels stacked together (by the increasing order of resolution level).
`wave.family`	The family of wavelets to use - "DaubExPhase" or "DaubLeAsymm". Default is "DaubLeAsymm".
`filter.number`	The number of vanishing moments of the wavelet. Default is 6.
`bc`	The boundary condition to use - "periodic" or "symmetric". Default is "periodic".

Details

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

Value

A list containing the following.

groupcoef

A matrix of dimension (grid, grid), containing the posterior group coefficients obtained by BHMSMA methodology.

Author(s)

Nilotpal Sanyal, Marco Ferreira

Maintainer: Nilotpal Sanyal <[email protected]>

References

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

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
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)
postmixprob(n, grid, waveletcoefmat, hyperparam, analysis)

Arguments

`n`	Number of subjects.
`grid`	The number of voxels in one row (or, one column) of the brain slice of interest. Must be a power of 2. The total number of voxels is `grid^2`. The maximum value of `grid` for this package is 512.
`waveletcoefmat`	A matrix of dimension `(n,grid^2-1)`, containing for each subject the wavelet coefficients of all levels stacked together (by the increasing order of resolution level).
`hyperparam`	A vector containing the estimates of the six hyperparameters.
`analysis`	"MSA" or "SSA", depending on whether performing multi-subject analysis or single subject analysis.

Value

A list containing the following.

pkljbar

A matrix of dimension (n,grid^2-1), containing the piklj bar values.

Author(s)

Nilotpal Sanyal, Marco Ferreira

Maintainer: Nilotpal Sanyal <[email protected]>

References

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

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
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)
postsamples(nsample, n, grid, glmcoefstd, waveletcoefmat, 
hyperparam, pkljbar, analysis, wave.family="DaubLeAsymm", 
filter.number=6, bc="periodic", seed)

Arguments

`nsample`	Number of samples to be generated.
`n`	Number of subjects.
`grid`	The number of voxels in one row (or, one column) of the brain slice of interest. Must be a power of 2. The total number of voxels is `grid^2`. The maximum value of `grid` for this package is 512.
`glmcoefstd`	An array of dimension `(n,grid,grid)`, containing for each subject the standardized GLM coefficients obtained by fitting GLM to the time-series corresponding to the voxels.
`waveletcoefmat`	A matrix of dimension `(n,grid^2-1)`, containing for each subject the wavelet coefficients of all levels stacked together (by the increasing order of resolution level).
`hyperparam`	A vector containing the estimates of the six hyperparameters.
`pkljbar`	A matrix of dimension `(n,grid^2-1)`, containing the piklj bar values (see References for details).
`analysis`	"MSA" or "SSA", depending on whether performing multi-subject analysis or single subject analysis.
`wave.family`	The family of wavelets to use - "DaubExPhase" or "DaubLeAsymm". Default is "DaubLeAsymm".
`filter.number`	The number of vanishing moments of the wavelet. Default is 6.
`bc`	The boundary condition to use - "periodic" or "symmetric". Default is "periodic".
`seed`	Must be a positive integer. Provide to set random number generation seed for reproducibility.

Details

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

Value

A list containing the following.

`samples`	An array of dimension `(n,grid,grid,nsample)`, containing for each subject the posterior samples of the GLM coefficients.
`postdiscovery`	An array of dimension `(n,grid,grid)`, containing for each subject the posterior discovery maps of the GLM coefficients (for details see Morris et al. (2011)).

Author(s)

Nilotpal Sanyal, Marco Ferreira

Maintainer: Nilotpal Sanyal <[email protected]>

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.

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
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)
postwaveletcoef(n, grid, waveletcoefmat, hyperparam, 
pkljbar, analysis)

Arguments

`n`	Number of subjects.
`grid`	The number of voxels in one row (or, one column) of the brain slice of interest. Must be a power of 2. The total number of voxels is `grid^2`. The maximum value of `grid` for this package is 512.
`waveletcoefmat`	A matrix of dimension `(n,grid^2-1)`, containing for each subject the wavelet coefficients of all levels stacked together (by the increasing order of resolution level).
`hyperparam`	A vector containing the estimates of the six hyperparameters.
`pkljbar`	A matrix of dimension `(n,grid^2-1)`, containing the piklj bar values.
`analysis`	"MSA" or "SSA", depending on whether performing multi-subject analysis or single subject analysis.

Value

A list containing the following.

`PostMeanWaveletCoef`	A matrix of size `(n,grid^2-1)`, containing for each subject the posterior mean of the wavelet coefficients of all levels stacked together (by the increasing order of resolution level).
`PostMedianWaveletCoef`	A matrix of size `(n,grid^2-1)`, containing for each subject the posterior median of the wavelet coefficients of all levels stacked together.

Author(s)

Nilotpal Sanyal, Marco Ferreira

Maintainer: Nilotpal Sanyal <[email protected]>

References

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

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
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 )
readfmridata( directory, format, prefix, nimages, dim.image, 
  nii=TRUE )

Arguments

`directory`	Location of the directory where the fMRI image files are stored. Insert within quotations ("").
`format`	The format of the data file. One of "Analyze" (.img/.hdr files), "Nifti" (.img/.hdr files or .nii files) or "Afni" (.HEAD/.BRIK files).
`prefix`	If format is "Analyze" or "Nifti", then the part of the fMRI image file name appearing before the image number. The image number is assumed to have four digit representation, that is, the lowest number is 0001 and the highest possible number is 9999. If format is "Afni", then the file name. Insert within quotations("").
`nimages`	If format is "Analyze", number of images to be read beginning from the first image. If format is "Afni", not necessary.
`dim.image`	Size of the 3D fMRI image. A vector with three elements.
`nii`	Necessary only for "Nifti" format. nii=TRUE (default) indicates the image files are in .nii files . nii=FALSE indicates the image files are .img/.hdr files.

Details

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

Value

A list containing the following:

fmridata

An array of dimension (dim.image, nimages), containing the image data for all images/time-points.

Author(s)

Nilotpal Sanyal, Marco Ferreira Maintainer: Nilotpal Sanyal <[email protected]>

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)
# 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)
substituteWaveletCoef(grid, waveletobj, values)

Arguments

`grid`	The number of voxels in one row (or, one column) of the brain slice of interest. Must be a power of 2. The total number of voxels is `grid^2`. The maximum value of `grid` for this package is 512.
`waveletobj`	A wavelet object of class `imwd.object`, usually obtained by performing 2D wavelet transformion through the function `imwd` of package wavethresh.
`values`	The values with which the wavelet coefficients are to be replaced. The order should be consistent with `imwd.object` class.

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 <[email protected]>

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
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")
waveletcoef(n, grid, glmcoefstd, wave.family="DaubLeAsymm", 
  filter.number=6, bc="periodic")

Arguments

`n`	Number of subjects.
`grid`	The number of voxels in one row (or, one column) of the brain slice of interest. Must be a power of 2. The total number of voxels is `grid^2`. The maximum value of `grid` for this package is 512.
`glmcoefstd`	An array of dimension `(n,grid,grid)`, containing for each subject the standardized GLM coefficients of a regressor obtained by fitting GLM to the time-series corresponding to the voxels.
`wave.family`	The family of wavelets to use - "DaubExPhase" or "DaubLeAsymm". Default is "DaubLeAsymm".
`filter.number`	The number of vanishing moments of the wavelet. Default is 6.
`bc`	The boundary condition to use - "periodic" or "symmetric". Default is "periodic".

Details

The wavelet decomposition is performed by using the function imwd.

Value

A list containing the following.

WaveletCoefficientMatrix

A matrix of dimension (n,grid^2-1), containing for each subject the wavelet coefficients of all levels stacked together (by the increasing order of resolution level).

Author(s)

Nilotpal Sanyal, Marco Ferreira

Maintainer: Nilotpal Sanyal <[email protected]>

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

Package 'BHMSMAfMRI'

Help Index

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

Description

Details

Author(s)

References

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

Description

Usage

Arguments

Details

Value

Author(s)

References

See Also

Examples

A simulated fMRI data for 3 subjects

Description

Usage

Format

Details

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

Description

Usage

Arguments

Value

Author(s)

References

See Also

Examples

Obtain estimates of the hyperparameters of the BHMSMA model

Description

Usage

Arguments

Value

Author(s)

References

See Also

Examples

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

Description

Usage

Arguments

Details

Value

Author(s)

See Also

Examples

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

Description

Usage

Arguments

Details

Value

Author(s)

References

See Also

Examples

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

Description

Usage

Arguments

Value

Author(s)

References

See Also

Examples

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

Description

Usage

Arguments

Details

Value

Author(s)

References

See Also

Examples

Obtain posterior estimates of the BHMSMA wavelet coefficients

Description