fmri.lmePar: Linear Mixed-effects Model for fMRI data
In neuroconductor/fmri: Analysis of fMRI Experiments

Description Usage Arguments Details Value Note Author(s) References See Also Examples

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.

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))

`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.

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.

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

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

Sibylle Dames

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.

lme, fmri.designG, fmri.design, fmri.stimulus, fmri.metaPar

## 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)] <- extractData(ds1)
Bold[1:dx,1:dy,1:dz,(dt*1+1):(dt*2)] <- extractData(ds2)
Bold[1:dx,1:dy,1:dz,(dt*2+1):(dt*3)] <- extractData(ds3)
Bold[1:dx,1:dy,1:dz,(dt*3+1):(dt*4)] <- extractData(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)