Linear Mixed-effects 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 Mixed-effects Model with temporally correlated and heteroscedastic within-subject errors. Voxel-wise regression analysis is accelerated by optional parallel processing using R package parallel.

Usage

1
2
3
fmri.lmePar(bold, z, fixed = NULL, random = NULL, mask = NULL, 
            ac = 0.3, vtype = "individual", cluster = 2, 
            wghts = c(1, 1, 1))

Arguments

bold

a large 4D-Array 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". If not the whole brain but a region is analyzed, vectors with region-indices can be preserved by adding as attributes (e.g. attr(bold, "xind") <- xind).

z

a design matrix for a multi-subject and/or multi-session fMRI-study of class "data.frame" specifying the expected BOLD response(s) and additional components for trend and other effects. Typically a fmri.designG object. This data frame contains all variables named in the model. There are some indispensable variables: "group", "subj", "session" and "run", which define the different strata. That information will be used for setting up the residual variance structure.

fixed

optionally, a one-sided linear formula describing the fixed-effects part of the model. Default settings are: fixed <- ~ 0 + hrf + session + drift1:session + drift2:session in case of one detected group, and the same but "hrf" replaced with "hrf:group" if two group levels in z are found. Since an intercept would be a linear combination of the session factor-variable modeling session-specific intercepts, it is excluded.

random

optionally, a one-sided formula of the form ~ x1 + ... + xn | g1/.../gm, with ~ x1 + ... + xn specifying the model for the random effects and g1/.../gm the grouping structure.

Default is always the basic model without covariates, i.e.
random <- ~ 0 + hrf|subj if no repeated measures in z are found (nlevels(z$run)==1),
random <- ~ 0 + hrf|subj/session if repeated measures and
random <- ~ 0 + hrf|session if repeated measures but one subject only.
In case of two independent groups:
random <- list(subj = pdDiag(~ 0 + hrf:group)) is used.

mask

if available, a logical 3D-Array 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 3D-Array 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 fmri.lm applied to the first subject. Alternatively, a global approach with uniform value can be used. In this case enter a number between 0 and 1. Default is 0.3 applied to all voxels.

vtype

a character string choosing the residual variance model. If "equal", homoscedastic variance across subjects is assumed setting weights argument in function lme() to zero, whereas "individual" allows different within-subject variances. Default method is "individual" that means subject-specific error variances using formula: weights <- varIdent(form =~ 1|subj).

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.

Details

fmri.lmePar() fits the configured Linear Mixed-effects 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). Voxel-by-voxel 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 fixed-effects coefficient and corresponding parameters are stored in results object. This should be the first specified predictor in the fixed-effects part of the model (verify the attribute of "df" in returned object). However, in two-sample case this principle does not work. The order changes, estimated session-specific 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 Mixed-effects Model up to scale factor resscale separated for the groups 1 and 2

resscale, resscale2

resscale*extract.data(object,"residuals") are the residuals of group 1 and group 2 respectively.

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 z-direction 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 t-statistics reported in lme() objects for the extracted regression coefficients separated for the groups 1 and 2. The name of the coefficient belonging to this df-value appears as attribute.

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

fixed-effects model

randomModel

random-effects model

VarModel

assumption about the subject error variances

cluster

number of threads run in parallel

attr(*,"design")

design matrix for the multi-subject fMRI-study

attr(*,"approach")

one-stage 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 two-stage 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). Mixed-Effects Models in S and S-Plus. 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.1-117.

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 multi-subject study
subj <- 4
runs <- 1
z <-fmri.designG(hrf, subj = subj, runs = runs)

## Assembly of the aggregated BOLD-Array
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 + hrf|subj,           
            correlation = corAR1(value = 0.3, form = ~ 1|subj/session, fixed=TRUE),
            weights = varIdent(form =~ 1|subj), 
            method ="REML",
            control = lmeControl(rel.tol=1e-6, returnObject = TRUE),
            data = z)
summary(M1.1)

# Residual plots
plot(M1.1, resid(.,type = "response") ~ scan|subj)
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 = "BOLD-signal", 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 = "BOLD-signal", 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 = "BOLD-signal", 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 = "BOLD-signal", 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)