fmri (version 1.9.12.1)

fmri.metaPar: Linear Mixed-effects Meta-Analysis model for fMRI data

Description

Group maps are estimated from BOLD effect estimates and their variances previously determined for each subject. The function rma.uni from R package metafor is used to fit mixed-effects meta-analytic models at group level. Voxel-wise regression analysis is accelerated by optional parallel processing using R package parallel.

Usage

fmri.metaPar(Cbold, Vbold, XG = NULL, model = NULL, method = "REML",
             weighted = TRUE, knha = FALSE, mask = NULL, cluster = 2,
             wghts = c(1, 1, 1))

Value

An object of class "fmrispm" and "fmridata", basically a list with components:

beta

estimated regression coefficients

se

estimated standard errors of the coefficients

cbeta

estimated BOLD contrast parameters for the group. Always the first regression coefficient is taken.

var

estimated variance of the BOLD contrast parameters

mask

brain mask

residuals

raw (integer size 2) vector containing residuals of the estimated linear mixed-effects meta-analytic model up to scale factor resscale

resscale

resscale*extractData(object,"residuals") are the residuals.

tau2

estimated amount of (residual) heterogeneity. Always 0 when method = "FE".

rxyz

array of smoothness from estimated correlation for each voxel in resel space (for analysis without smoothing).

scorr

array of spatial correlations with maximal lags 5, 5, 3 in x, y and z-direction

bw

vector of bandwidths (in FWHM) corresponding to the spatial correlation within the data

weights

ratio of voxel dimensions

dim

dimension of the data cube and residuals

df

degrees of freedom for t-statistics, df = (n-p-1)

sessions

number of observations entering the meta-analytic model, n

coef

number of coefficients in the meta-analytic model (including the intercept, p+1)

method

estimator used to fit the meta-analytic model. In case of "FE", it is weighted or unweighted least squares.

weighted

estimation with inverse-variance weights

knha

Knapp and Hartung adjustment

model

meta-analytic regression model

cluster

number of threads running in parallel

attr(*,"design")

group-level design matrix

attr(*,"approach")

two-stage estimation method

Arguments

Cbold

a 4D-Array with the aggregated individual BOLD contrast estimates in standard space, e.g. all cbeta maps obtained from single-session analysis with fmri.lm may put together. Dimensions 1 to 3 define the voxel space, dimension 4 indicates a subject. If not the whole brain but a region is analyzed, vectors with region-indices can be preserved by adding as attributes (e.g. attr(Cbold, "xind") <- xind).

Vbold

a 4D-Array with the aggregated variance estimates for the contrast parameters in Cbold, e.g. all var maps obtained from single-session analysis with fmri.lm may put together. Dimensions 1 to 3 define the voxel space, dimension 4 indicates a subject.

XG

optionally, a group-level design matrix of class "data.frame" to include one or more moderators in the model. By default, an intercept is added to the model.

model

optionally, a one-sided formula of the form: model <- ~ mod1 + mod2 + mod3 describing a model with moderator variables. Adding "-1" removes the intercept term.

method

a character string specifying whether a fixed- (method = "FE") or a random/mixed-effects model (method = "REML", default) should be fitted. Further estimators for random/mixed-effects models are available, see documentation of rma.uni function for more details.

weighted

logical indicating whether weighted (weighted = TRUE, default) or unweighted estimation should be used to fit the model.

knha

logical specifying whether the method by Knapp and Hartung (2003) should be used for adjusting standard errors of the estimated coefficients (default is FALSE). The Knapp and Hartung adjustment is only meant to be used in the context of random- or mixed-effects models.

mask

if available, a logical 3D-Array of dimensionality of the data (without 4th subject component) describing a brain mask. The computation is restricted to the selected voxels.

cluster

number of threads for parallel processing, which is limited to available multi-core CPUs. If you do not know your CPUs, try: detectCores() from parallel package. Presets are 2 threads. cluster = 1 does not use parallel package.

wghts

a vector of length 3 specifying ratio of voxel dimensions. Isotropic voxels (e.g. MNI-space) are set as default.

Author

Sibylle Dames

Details

fmri.metaPar() fits the configured linear mixed-effects meta-analytic (MEMA) model separately at each voxel and extracts the first regression coefficient (usually the overall group mean), corresponding squared standard errors and degrees of freedom as well as the residuals from resulting rma.uni objects, to obtain a statistical parametric map (SPM) for the group. Voxel-by-voxel analysis is performed by either the function apply or parApply from parallel package, which walks through the Cbold array.

This two-stage approach reduces the computational burden of fitting a full linear mixed-effects (LME) model, fmri.lmePar would do. It assumes first level design is same across subjects and normally distributed not necessarily homogeneous within-subject errors. Warping to standard space has been done before first-stage analyses are carried out. Either no masking or a uniform brain mask should be applied at individual subject analysis level, to avoid loss of information at group level along the edges.

At the second stage, observed individual BOLD effects from each study are combined in a meta-analytic model. There is the opportunity of weighting the fMRI studies by the precision of their respective effect estimate to take account of first level residual heterogeneity (weighted = TRUE). This is how to deal with intra-subject variability. The REML estimate of cross-subject variability (tau-squared) assumes that each of these observations is drawn independently from the same Gaussian distribution. Since correlation structures cannot be modeled, multi-subject fMRI studies with repeated measures cannot be analyzed in this way.

Spatial correlation among voxels, e.g. through the activation of nearby voxels, is ignored at this stage, but corrects for it, when random field theory define a threshold for significant activation at inference stage.

It is recommended to check your model syntax and residuals choosing some distinct voxels before running the model in loop (see Example). Error handling default is to stop if one of the threads produces an error. When this occurs, the output will be lost from any voxel, where the model has fitted successfully.

References

Chen G., Saad Z.S., Nath A.R., Beauchamp M.S., Cox R.W. (2012). FMRI group analysis combining effect estimates and their variances. NeuroImage, 60: 747-765.

Knapp G. and Hartung J. (2003). Improved tests for a random effects meta-regression with a single covariate. Statistics in Medicine, 22: 2693-2710.

Viechtbauer W. (2005). Bias and efficiency of meta-analytic variance estimators in the random-effects model. Journal of Educational and Behavioral Statistics, 30: 261-293.

Viechtbauer W. (2010). Conducting meta-analyses in R with the metafor package. Journal of Statistical Software, 36(3): 1-48

Viechtbauer W. (2015). metafor: Meta-Analysis Package for R R package version 1.9-7.

See Also

rma.uni, fmri.lm, fmri.lmePar

Examples

if (FALSE) ## Generate some fMRI data sets: noise + stimulus
dx <- dy <- dz <- 32
dt <- 107
hrf <- fmri.stimulus(dt, c(18, 48, 78), 15, 2)
stim <- matrix(hrf, nrow= dx*dy*dz, ncol=dt, byrow=TRUE)
mask <- array(FALSE, c(dx, dy, dz))
mask[12:22,12:22,12:22] <- TRUE

ds1 <- list(ttt=writeBin(1.0*rnorm(dx*dy*dz*dt) + as.vector(5*stim),
           raw(), 4), mask = mask, dim = c(dx, dy, dz, dt))
ds2 <- list(ttt=writeBin(1.7*rnorm(dx*dy*dz*dt) + as.vector(3*stim),
           raw(), 4), mask = mask, dim = c(dx, dy, dz, dt))
ds3 <- list(ttt=writeBin(0.8*rnorm(dx*dy*dz*dt) + as.vector(1*stim),
           raw(), 4), mask = mask, dim = c(dx, dy, dz, dt))
ds4 <- list(ttt=writeBin(1.2*rnorm(dx*dy*dz*dt) + as.vector(2*stim),
           raw(), 4), mask = mask, dim = c(dx, dy, dz, dt))
class(ds1) <- class(ds2) <- class(ds3) <- class(ds4) <- "fmridata"

## Stage 1: single-session regression analysis
x <- fmri.design(hrf, order=2)
spm.sub01 <- fmri.lm(ds1, x, mask, actype = "smooth", verbose = TRUE)
spm.sub02 <- fmri.lm(ds2, x, mask, actype = "smooth", verbose = TRUE)
spm.sub03 <- fmri.lm(ds3, x, mask, actype = "smooth", verbose = TRUE)
spm.sub04 <- fmri.lm(ds4, x, mask, actype = "smooth", verbose = TRUE)

## Store observed individual BOLD effects and their variance estimates
subj <- 4
Cbold <- array(0, dim = c(dx, dy, dz, subj))
Cbold[,,,1] <- spm.sub01$cbeta
Cbold[,,,2] <- spm.sub02$cbeta
Cbold[,,,3] <- spm.sub03$cbeta
Cbold[,,,4] <- spm.sub04$cbeta

Vbold <- array(0, dim = c(dx, dy, dz, subj))
Vbold[,,,1] <- spm.sub01$var
Vbold[,,,2] <- spm.sub02$var
Vbold[,,,3] <- spm.sub03$var
Vbold[,,,4] <- spm.sub04$var

## Stage 2: Random-effects meta-regression analysis
## a) Check your model
library(metafor)
M1.1 <- rma.uni(Cbold[16,16,16, ],
                Vbold[16,16,16, ],
                method = "REML",
                weighted = TRUE,
                knha = TRUE,
                verbose = TRUE,
                control = list(stepadj=0.5, maxiter=2000, threshold=0.001))

# Control list contains convergence parameters later used
# at whole data cube. Values were adjusted to fMRI data.

summary(M1.1)
forest(M1.1)
qqnorm(M1.1)

## b) Estimate a group map
## without parallelizing
spm.group1a <- fmri.metaPar(Cbold, Vbold, knha = TRUE,
                            mask = mask, cluster = 1)
## same with 4 parallel threads
spm.group1b <- fmri.metaPar(Cbold, Vbold, knha = TRUE,
                            mask = mask, cluster = 4)