Linear Mixedeffects Model for fMRI data
Description
Group maps are directly estimated from the BOLD time series data of all subjects using lme
from R package nlme to fit a Linear Mixedeffects Model with temporally correlated and heteroscedastic withinsubject errors. Voxelwise regression analysis is accelerated by optional parallel processing using R package parallel.
Usage
1 2 3 
Arguments
bold 
a large 4DArray with the aggregated fMRI data of all subjects that were previously registered to a common brain atlas. Be careful with the assembly of this array, the order of the data sets has to be compatible with the design matrix: 
z 
a design matrix for a multisubject and/or multisession fMRIstudy of class 
fixed 
optionally, a onesided linear formula describing the fixedeffects part of the model. Default settings are:

random 
optionally, a onesided formula of the form Default is always the basic model without covariates, i.e. 
mask 
if available, a logical 3DArray of dimensionality of the data (without time component) describing a brain mask. The computation is restricted to the selected voxels. 
ac 
if available, a numeric 3DArray of dimensionality of the data (without time component) with spatially smoothed autocorrelation parameters should be used in the AR(1) models fitted in each voxel, e.g. locally estimated and smoothed AR(1)coefficients from 
vtype 
a character string choosing the residual variance model. If 
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. 
Details
fmri.lmePar()
fits the configured Linear Mixedeffects Model separately at each voxel and extracts estimated BOLD contrasts, corresponding squared standard errors and degrees of freedom as well as the residuals from resulting lme()
objects to produce a statistical parametric map (SPM) for the group(s). Voxelbyvoxel analysis is performed by either the function apply
or parApply
from parallel package, which walks through the bold
array.
If one group is analyzed, from each fitted model the first fixedeffects coefficient and corresponding parameters are stored in results object. This should be the first specified predictor in the fixedeffects part of the model (verify the attribute of "df"
in returned object). However, in twosample case this principle does not work. The order changes, estimated sessionspecific intercepts now comes first and the number of these coefficients is not fixed. Therefore in current version it has explicitly been looked for the coefficient names: "hrf:group1" and "hrf:group2". Available functions within the nlme package to extract estimated values from lme()
objects do not operate at contrast matrices.
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, step 1); especially for more advanced designs! 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.
Value
An object of class "fmrispm"
and "fmridata"
, basically a list
with components:
cbeta, cbeta2 
estimated BOLD contrast parameters separated for the groups 1 and 2 
var, var2 
estimated variance of the contrast parameters separated for the groups 1 and 2 
mask 
brain mask 
res, res2 
raw (integer size 2) vector containing residuals of the estimated Linear Mixedeffects Model up to scale factor 
resscale, resscale2 

arfactor 
autocorrelation parameters used in AR(1)model 
rxyz, rxyz2 
array of smoothness from estimated correlation for each voxel in resel space separated for the groups 1 and 2 (for analysis without smoothing) 
scorr, scorr2 
array of spatial correlations with maximal lags 5, 5, 3 in x, y and zdirection separated for the groups 1 and 2 
bw, bw2 
vector of bandwidths (in FWHM) corresponding to the spatial correlation within the data separated for the groups 1 and 2 
weights 
ratio of voxel dimensions 
dim, dim2 
dimension of the data cube and residuals separated for the groups 1 and 2 
df, df2 
degrees of freedom for tstatistics reported in 
subjects 
number of subjects in the study 
subj.runs 
number of repeated measures within subjects 
sessions 
number of total sessions that were analyzed 
groups 
number of groups in the study 
fixedModel 
fixedeffects model 
randomModel 
randomeffects model 
VarModel 
assumption about the subject error variances 
cluster 
number of threads run in parallel 
attr(*,"design") 
design matrix for the multisubject fMRIstudy 
attr(*,"approach") 
onestage estimation method 
Note
Maybe the computing power is insufficient to carry out a whole brain analysis. You have two opportunities: either select and analyze a certain brain area or switch to a twostage model.
Current Limitations
The function cannot handle experimental designs with:
more than two independent groups
more than one stimulus (task)
paired samples with varying tasks
user defined contrasts
Author(s)
Sibylle Dames
References
Pinheiro J. and Bates D. (2000). MixedEffects Models in S and SPlus. Springer.
Pinheiro J., Bates D., DebRoy S., Sarkar D. and the R Core team (2014). nlme: Linear and Nonlinear Mixed Effects Models R package version 3.1117.
See Also
lme
, fmri.designG
,
fmri.design
, fmri.stimulus
,
fmri.metaPar
Examples
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 66 67 68 69 70 71 72 73 74 75 76 77  ## 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"
## Construct a design matrix for a multisubject study
subj < 4
runs < 1
z <fmri.designG(hrf, subj = subj, runs = runs)
## Assembly of the aggregated BOLDArray
Bold < array(0, dim = c(dx,dy,dz,subj*runs*dt))
Bold[1:dx,1:dy,1:dz,1:(dt*1)] < extract.data(ds1)
Bold[1:dx,1:dy,1:dz,(dt*1+1):(dt*2)] < extract.data(ds2)
Bold[1:dx,1:dy,1:dz,(dt*2+1):(dt*3)] < extract.data(ds3)
Bold[1:dx,1:dy,1:dz,(dt*3+1):(dt*4)] < extract.data(ds4)
## Step 1: Check the model
y < Bold[16, 16, 16, ] # choose one voxel
M1.1 < lme(fixed = y ~ 0 + hrf + session + drift1:session + drift2:session,
random = ~ 0 + hrfsubj,
correlation = corAR1(value = 0.3, form = ~ 1subj/session, fixed=TRUE),
weights = varIdent(form =~ 1subj),
method ="REML",
control = lmeControl(rel.tol=1e6, returnObject = TRUE),
data = z)
summary(M1.1)
# Residual plots
plot(M1.1, resid(.,type = "response") ~ scansubj)
qqnorm(M1.1, ~resid(.,type = "normalized")subj, abline = c(0,1))
# Testing the assumption of homoscedasticity
M1.2 < update(M1.1, weights = NULL, data = z)
anova(M1.2, M1.1)
# Model fit: observed and fitted values
fitted.values < fitted(M1.1)
plot(y[1:dt], type="l", main = "Subject 1", xlab = "scan",
ylab = "BOLDsignal", ylim = c(5,5))
lines(fitted.values[names(fitted.values)==1],lty=1,lwd=2)
plot(y[(dt+1):(2*dt)], type="l", main = "Subject 2", xlab = "scan",
ylab = "BOLDsignal", ylim = c(5,5))
lines(fitted.values[names(fitted.values)==2],lty=1,lwd=2)
plot(y[(2*dt+1):(3*dt)], type="l", main = "Subject 3", xlab = "scan",
ylab = "BOLDsignal", ylim = c(5,5))
lines(fitted.values[names(fitted.values)==3],lty=1,lwd=2)
plot(y[(3*dt+1):(4*dt)], type="l", main = "Subject 4", xlab = "scan",
ylab = "BOLDsignal", ylim = c(5,5))
lines(fitted.values[names(fitted.values)==4],lty=1,lwd=2)
## Step 2: Estimate a group map
## without parallelizing
spm.group1a < fmri.lmePar(Bold, z, mask = mask, cluster = 1)
# same with 4 parallel threads
spm.group1b < fmri.lmePar(Bold, z, mask = mask, cluster = 4)
## Example for two independent groups
group < c(1,1,4,4)
z2 < fmri.designG(hrf, subj = subj, runs = runs, group = group)
spm.group2 < fmri.lmePar(Bold, z2, mask = mask, cluster = 4)
## End(Not run)
