Description Usage Arguments Details Value Note Author(s) References See Also Examples
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 mixedeffects metaanalytic models at group level. Voxelwise regression analysis is accelerated by optional parallel processing using R package parallel.
1 2 3 
Cbold 
a 4DArray with the aggregated individual BOLD contrast estimates in standard space, e.g. all 
Vbold 
a 4DArray with the aggregated variance estimates for the contrast parameters in 
XG 
optionally, a grouplevel design matrix of class 
model 
optionally, a onesided formula of the form: 
method 
a character string specifying whether a fixed (method = "FE") or a random/mixedeffects model (method = "REML", default) should be fitted. Further estimators for random/mixedeffects models are available, see documentation of 
weighted 
logical indicating whether weighted ( 
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 mixedeffects models. 
mask 
if available, a logical 3DArray 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 multicore CPUs. If you do not know your CPUs, try: 
wghts 
a vector of length 3 specifying ratio of voxel dimensions. Isotropic voxels (e.g. MNIspace) are set as default. 
fmri.metaPar()
fits the configured linear mixedeffects metaanalytic (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. Voxelbyvoxel analysis is performed by either the function apply
or parApply
from parallel package, which walks through the Cbold
array.
This twostage approach reduces the computational burden of fitting a full linear mixedeffects (LME) model, fmri.lmePar
would do. It assumes first level design is same across subjects and normally distributed not necessarily homogeneous withinsubject errors. Warping to standard space has been done before firststage 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 metaanalytic 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 intrasubject variability. The REML estimate of crosssubject variability (tausquared) assumes that each of these observations is drawn independently from the same Gaussian distribution. Since correlation structures cannot be modeled, multisubject 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.
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 
res 
raw (integer size 2) vector containing residuals of the estimated linear mixedeffects metaanalytic model up to scale factor 
resscale 

tau2 
estimated amount of (residual) heterogeneity. Always 0 when 
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 zdirection 
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 tstatistics, df = (np1) 
sessions 
number of observations entering the metaanalytic model, n 
coef 
number of coefficients in the metaanalytic model (including the intercept, p+1) 
method 
estimator used to fit the metaanalytic model. In case of "FE", it is weighted or unweighted least squares. 
weighted 
estimation with inversevariance weights 
knha 
Knapp and Hartung adjustment 
model 
metaanalytic regression model 
cluster 
number of threads running in parallel 
attr(*,"design") 
grouplevel design matrix 
attr(*,"approach") 
twostage estimation method 
Meta analyses tend to be less powerful for neuroimaging studies, because they only have as many degrees of freedom as number of subjects. If the number of subjects is very small, then it may be impossible to estimate the betweensubject variance (tausquared) with any precision. In this case the fixed effect model may be the only viable option. However, there is also the possibility of using a onestage model, that includes the full time series data from all subjects and simultaneously estimates subject and group levels parameters (see fmri.lmePar
). Although this approach is much more computer intensive, it has the advantage of higher degrees of freedom (> 100) at the end.
Current Limitations
The function cannot handle:
experimental designs with a withinsubject (repeated measures) factor
paired samples with varying tasks, unless the contrast of the two conditions is used as input
Sibylle Dames
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: 747765.
Knapp G. and Hartung J. (2003). Improved tests for a random effects metaregression with a single covariate. Statistics in Medicine, 22: 26932710.
Viechtbauer W. (2005). Bias and efficiency of metaanalytic variance estimators in the randomeffects model. Journal of Educational and Behavioral Statistics, 30: 261293.
Viechtbauer W. (2010). Conducting metaanalyses in R with the metafor package. Journal of Statistical Software, 36(3): 148
Viechtbauer W. (2015). metafor: MetaAnalysis Package for R R package version 1.97.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65  ## Not run: ## 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: singlesession 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: Randomeffects metaregression 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)
## End(Not run)

Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.