Linear Mixed-effects Model for fMRI data

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.

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*extractData(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

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.

Author

Sibylle Dames

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.

more than two independent groups
more than one stimulus (task)
paired samples with varying tasks
user defined contrasts

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)
#> Error in matrix(hrf, nrow = dx * dy * dz, ncol = dt, byrow = TRUE) object 'dx' not found
mask <- array(FALSE, c(dx, dy, dz))
#> Error in array(FALSE, c(dx, dy, dz)) object 'dx' not found
mask[12:22,12:22,12:22] <- TRUE
#> Error in mask[12:22, 12:22, 12:22] <- TRUE object 'mask' not found

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))
#> Error in rnorm(dx * dy * dz * dt) object 'dx' not found
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))
#> Error in rnorm(dx * dy * dz * dt) object 'dx' not found
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))
#> Error in rnorm(dx * dy * dz * dt) object 'dx' not found
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))
#> Error in rnorm(dx * dy * dz * dt) object 'dx' not found
class(ds1) <- class(ds2) <- class(ds3) <- class(ds4) <- "fmridata"
#> Error in class(ds4) <- "fmridata" object 'ds4' not found

## 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))
#> Error in array(0, dim = c(dx, dy, dz, subj * runs * dt)) object 'dx' not found
Bold[1:dx,1:dy,1:dz,1:(dt*1)] <- extractData(ds1)
#> Error in extractData(ds1) object 'ds1' not found
Bold[1:dx,1:dy,1:dz,(dt*1+1):(dt*2)] <- extractData(ds2)
#> Error in extractData(ds2) object 'ds2' not found
Bold[1:dx,1:dy,1:dz,(dt*2+1):(dt*3)] <- extractData(ds3)
#> Error in extractData(ds3) object 'ds3' not found
Bold[1:dx,1:dy,1:dz,(dt*3+1):(dt*4)] <- extractData(ds4)
#> Error in extractData(ds4) object 'ds4' not found

## Step 1: Check the model
y <- Bold[16, 16, 16, ] # choose one voxel
#> Error in eval(expr, envir, enclos) object 'Bold' not found
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)
#> Error in 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-06,         returnObject = TRUE), data = z) could not find function "lme"
summary(M1.1)
#> Error in summary(M1.1) object 'M1.1' not found

# Residual plots
plot(M1.1, resid(.,type = "response") ~ scan|subj)
#> Error in plot(M1.1, resid(., type = "response") ~ scan | subj) object 'M1.1' not found
qqnorm(M1.1, ~resid(.,type = "normalized")|subj, abline = c(0,1))
#> Error in qqnorm(M1.1, ~resid(., type = "normalized") | subj, abline = c(0,     1)) object 'M1.1' not found

# Testing the assumption of homoscedasticity
M1.2 <- update(M1.1, weights = NULL, data = z)
#> Error in update(M1.1, weights = NULL, data = z) object 'M1.1' not found
anova(M1.2, M1.1)
#> Error in anova(M1.2, M1.1) object 'M1.2' not found

# Model fit: observed and fitted values
fitted.values <- fitted(M1.1)
#> Error in fitted(M1.1) object 'M1.1' not found
plot(y[1:dt], type="l", main = "Subject 1", xlab = "scan",
     ylab = "BOLD-signal", ylim = c(-5,5))
#> Error in plot(y[1:dt], type = "l", main = "Subject 1", xlab = "scan",     ylab = "BOLD-signal", ylim = c(-5, 5)) object 'y' not found
lines(fitted.values[names(fitted.values)==1],lty=1,lwd=2)
#> Error in fitted.values[names(fitted.values) == 1] object of type 'closure' is not subsettable

plot(y[(dt+1):(2*dt)], type="l", main = "Subject 2", xlab = "scan",
     ylab = "BOLD-signal", ylim = c(-5,5))
#> Error in plot(y[(dt + 1):(2 * dt)], type = "l", main = "Subject 2", xlab = "scan",     ylab = "BOLD-signal", ylim = c(-5, 5)) object 'y' not found
lines(fitted.values[names(fitted.values)==2],lty=1,lwd=2)
#> Error in fitted.values[names(fitted.values) == 2] object of type 'closure' is not subsettable

plot(y[(2*dt+1):(3*dt)], type="l", main = "Subject 3", xlab = "scan",
     ylab = "BOLD-signal", ylim = c(-5,5))
#> Error in plot(y[(2 * dt + 1):(3 * dt)], type = "l", main = "Subject 3",     xlab = "scan", ylab = "BOLD-signal", ylim = c(-5, 5)) object 'y' not found
lines(fitted.values[names(fitted.values)==3],lty=1,lwd=2)
#> Error in fitted.values[names(fitted.values) == 3] object of type 'closure' is not subsettable

plot(y[(3*dt+1):(4*dt)], type="l", main = "Subject 4", xlab = "scan",
     ylab = "BOLD-signal", ylim = c(-5,5))
#> Error in plot(y[(3 * dt + 1):(4 * dt)], type = "l", main = "Subject 4",     xlab = "scan", ylab = "BOLD-signal", ylim = c(-5, 5)) object 'y' not found
lines(fitted.values[names(fitted.values)==4],lty=1,lwd=2)
#> Error in fitted.values[names(fitted.values) == 4] object of type 'closure' is not subsettable

## Step 2: Estimate a group map
## without parallelizing
spm.group1a <- fmri.lmePar(Bold, z, mask = mask, cluster = 1)
#> fmri.lme: entering function
#> Error in fmri.lmePar(Bold, z, mask = mask, cluster = 1) object 'Bold' not found
# same with 4 parallel threads
spm.group1b <- fmri.lmePar(Bold, z, mask = mask, cluster = 4)
#> fmri.lme: entering function
#> Error in fmri.lmePar(Bold, z, mask = mask, cluster = 4) object 'Bold' not found
## 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)
#> fmri.lme: entering function
#> Error in fmri.lmePar(Bold, z2, mask = mask, cluster = 4) object 'Bold' not found