Convert Between fmridata and oro.nifti Objects

Description

NIfTI data can be converted between fmridata S3 objects (from the fmri package) and nifti S4 objects.

Usage

oro2fmri(from, value = NULL, level = 0.75, mask=NULL, setmask = TRUE)
fmri2oro(from, value = NULL, verbose = FALSE, reorient = FALSE,
         call = NULL)

Arguments

Details

These functions enhance the capabilities of fmri by allowing the exchange of data objects between nifti and fmridata classes.

Value

The function oro2fmri produces an S3 object of class fmridata. The function fmri2oro produces an S4 object of class nifti.

Author(s)

Brandon Whitcher bwhitcher@gmail.com

See Also

read.NIFTI

IC fingerprinting

Description

Implements ICA fingerprinting mainly following De Martino et.al., Neuroimage 2007.

Usage

ICAfingerprint(icaobj, nbin = 256, plot = FALSE)

Arguments

Details

For some characteristics normalization of values differs from De Martino et. al.. Frequency bands are obtained from periodogram estimated instead of using Welch's method.

Value

object of class ”fmriICA” list with components

Author(s)

Joerg Polzehl polzehl@wias-berlin.de

References

De Martino et. al., Classification of fMRI independent components using IC-fingerprints and support vector machine classifiers, Neuroimage 34 (2007) 177-194.

See Also

fmri.sICA, plot.fmriICA, fastICA

I/O function

Description

This functions cuts a region-of-interest (ROI) from input data.

Usage

cutroi(data, xind = 1:data$dim[1], yind = 1:data$dim[2], 
             zind = 1:data$dim[3], tind = 1:data$dim[4])

Arguments

Details

Cut a region of interest from frmidata.

Value

Corresponding cut fmridata object.

Author(s)

Karsten Tabelow tabelow@wias-berlin.de

See Also

read.AFNI, read.ANALYZE, read.NIFTI

Examples

##---- Should be DIRECTLY executable !! ----
##-- ==>  Define data, use random,
##--	or do  help(data=index)  for the standard data sets.

Extract data or residuals from a fmridata object

Description

The function extracts data stored as raw within an object of class 'fmridata'.

Usage

extractData(z, what = "data", maskOnly = FALSE)
expandfMRI(z)
condensefMRI(z, mask)

Arguments

Details

The function extractData extracts data stored as raw within an object of class 'fmridata'. Functions expandfMRI and condensefMRI change the way data and residuals are stored between full 3D data and data within a brain mask. condensefMRI can also be used to set a more restrictive brain mask.

Value

In case of function extractData an array of dimension data$dim containing either the fmri-data or residuals. The other two functions return an object of class 'fmridata'.

Author(s)

Joerg Polzehl polzehl@wias-berlin.de

See Also

fmri.lm

Cluster thresholding.

Description

Detection of activated regions using cluster thresholding.

Usage

fmri.cluster(spm, alpha = 0.05, ncmin = 2, ncmax=ncmin, 
             minimum.signal = 0, verbose = FALSE)

Arguments

Details

Approximate thresholds for the existence of a cluster with spm-values exceeding a 1-beta threshold k_{nc,na:ne} for cluster size nc are based on a simulation study under the hypothesis and adjusted for number of voxel in mask and spatial correlation. beta is chosen such that under the hypothesis the combined (over cluster sizes ncmin:ncmax) test has approximate significance level alpha.

Value

Object with class attributes "fmripvalue" and "fmridata"

Author(s)

Joerg Polzehl polzehl@wias-berlin.de

See Also

fmri.lm, fmri.pvalue, fmri.searchlight

Examples

  ## Not run: fmri.cluster(fmrispmobj)

Linear Model for FMRI Data

Description

Return a design matrix for a linear model with given stimuli and possible polynomial drift terms.

Usage

  fmri.design(stimulus, order = 2, cef = NULL, verbose = FALSE)

Arguments

Details

The stimuli given in stimulus are used as first columns in the design matrix.

The order of the polynomial drift terms is given by order, which defaults to 2.

Confounding effects can be included in a matrix cef.

The polynomials are defined orthogonal to the stimuli given in stimulus.

Value

design matrix of the linear model

Author(s)

Karsten Tabelow tabelow@wias-berlin.de, Joerg Polzehl polzehl@wias-berlin.de

References

Polzehl, J. and Tabelow, K.(2007). fmri: A Package for Analyzing fmri Data, R News, 7:13-17 .

See Also

fmri.stimulus, fmri.lm

Examples

  # Example 1
  hrf <- fmri.stimulus(107, c(18, 48, 78), 15, 2)
  z <- fmri.design(hrf, 2)
  par(mfrow=c(2, 2))
  for (i in 1:4) plot(z[, i], type="l")

Design matrix for fMRI group analysis

Description

This function returns a design matrix for multi-subject fMRI data to fit a Linear Mixed-effects Model (one-stage procedure) with given stimuli, polynomial drift terms and a set of known population parameters.

Usage

fmri.designG(hrf, subj = 1, runs = 1, group = NULL, XG = NULL)

Arguments

Details

Based on the dimensionality of the hrf object, which provides the total number of scans (time-points) within each session, the entered number of subjects and repeated measures the auxiliary variables: "subj", "run", "scan" and "session" are generated as first part of the returned design matrix.

If no group argument is specified, only one population will be assumed; otherwise group labels are replicated within sessions of the same subject.

First a design matrix for a single run is created by calling: x <- fmri.design(hrf, order = 2). Hence the polynomial drift terms are defined orthogonal to the stimuli (see fmri.design). This matrix is replicated blockwise to all sessions assuming the same experimental design for all runs. The first drift term, a column of ones, is called "drift0" and models an intercept.

If given, further subject characteristics are filled in the design matrix.

Value

A design matrix as a data frame, which contains the following variables:

Author(s)

Sibylle Dames

References

Polzehl, J. and Tabelow, K.(2007). fmri: A Package for Analyzing fmri Data, R News, 7:13-17.

See Also

fmri.stimulus, fmri.design, fmri.lmePar

Examples

subj <- 6
runs <- 1
scans <- 121
times <- c(12, 48, 84, 120, 156, 192, 228, 264) 
duration <- 24          
tr <- 2.5

hrf <- fmri.stimulus(scans, times, duration, tr, times = TRUE)
x.group <- fmri.designG(hrf, subj = subj, runs = runs)
# View(x.group)

Detrend fMRI time series

Description

Detrend fMRI dataset with a polynomial of given degree

Usage

fmri.detrend(data, degree = 1, nuisance=NULL, accoef = 0)

Arguments

Details

The function can be used to detrend the time series of an fMRI dataset data (of class ”fmridata” using polynomials. If the argument degree is larger than 0 (default: 1) the polynomial trends up to the given degree are removed from the data. If the argument accoef is larger than 0 (default: 0) prewhitening using an AR(1) model is performed.

Value

Detrended data object of class ”fmridata”.

Author(s)

Joerg Polzehl polzehl@wias-berlin.de

References

Polzehl, J. and Tabelow, K. (2007). fmri: A Package for Analyzing fmri Data, R News, 7:13-17.

See Also

fmri.lm

Examples

  # Example 1
  data <- list(ttt=writeBin(rnorm(32*32*32*107),raw(),4),
               mask=array(1,c(32,32,32)),dim=c(32,32,32,107))
  class(data) <- "fmridata"
  data <- fmri.detrend(data,2)

Linear Model for fMRI data

Description

Estimate the parameters and variances in a linear model.

Usage

  fmri.lm(ds, z, mask = NULL,
          actype = c("smooth", "noac", "ac", "accalc"),
          contrast = c(1), verbose = FALSE)

Arguments

Details

This function performs parameter estimation in the linear model. It implements a two step procedure. After primary estimation of the parameters in the first step residuals are obtained. If actype %in% c("ac", "accalc", "smooth") an AR(1) model is fitted, in each voxel, to the time series of residuals. The estimated AR-coefficients are corrected for bias. If actype=="smooth" the estimated AR-coefficients are spatially smoothed. If actype %in% c("ac", "smooth") the linear model is pre-whitened using the estimated (and possibly smoothed) AR-coefficients. Parameter and variance estimates are then obtained from the pre-whitened data. The argument keep describes the amount of data which is returned. The estimated effects

\tilde{\gamma}_i = C^T\tilde{\beta}_i

and their estimated variances are returned as well as the residuals and temporal autocorrelation. cbeta then contains the corresponding parameter estimates and thus is a vector of corresponding length in each voxel.

If z is an 3-dimensional array the third component is assumed to code the design matrix information for the corresponding slice, i.e. design matrices to differ with respect to slice timing effects. Note that if motion correction needs to be performed in preprocessing slice time correction may be better carried out on the data before image registration using, e.g., function slicetiming.

If warning "Local smoothness characterized by large bandwidth" occurs, check scorr elements. If correlation drops with lag towards zero, data has been pre-smoothed. Adaptive smoothing the SPM can then only be of limited use. If correlation does not go to zero, check the residuals of the linear model for unexplained structure (spin saturation in first scans? discard them!).

Value

object with class attributes "fmrispm" and "fmridata"

Author(s)

Karsten Tabelow tabelow@wias-berlin.de, Joerg Polzehl polzehl@wias-berlin.de

References

Worsley, K.J. (2005). Spatial smoothing of autocorrelations to control the degrees of freedom in fMRI analysis. NeuroImage, 26:635-641.

Worsley, K.J., Liao, C., Aston, J., Petre, V., Duncan, G.H., Morales, F., Evans, A.C. (2002). A general statistical analysis for fMRI data. NeuroImage, 15:1-15.

Tabelow, K., Polzehl, J., Voss, H.U., and Spokoiny, V. (2006). Analysing fMRI experiments with structure adaptive smoothing procedures, NeuroImage, 33:55-62.

See Also

fmri.design, fmri.stimulus

Examples

  ## Not run: 
  # Example 1
  data <- list(ttt=writeBin(rnorm(32*32*32*107), raw(), 4),
               mask=array(TRUE, c(32, 32, 32)), dim=c(32, 32, 32, 107))
  class(data) <- "fmridata"
  hrf <- fmri.stimulus(107, c(18, 48, 78), 15, 2)
  z <- fmri.design(hrf,2)
  model <- fmri.lm(data, z, verbose=TRUE)
  plot(extractData(data)[16, 16, 16,])
  lines(extractData(data)[16, 16, 16, ] - extractData(model, "residuals")[16, 16, 16, ], col=2)
  
## End(Not run)

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

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

Arguments

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:

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

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

Linear Mixed-effects Meta-Analysis model for fMRI data

Description

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 mixed-effects meta-analytic models at group level. Voxel-wise regression analysis is accelerated by optional parallel processing using R package parallel.

Usage

fmri.metaPar(Cbold, Vbold, XG = NULL, model = NULL, method = "REML",
             weighted = TRUE, knha = FALSE, mask = NULL, cluster = 2,
             wghts = c(1, 1, 1))

Arguments

Details

fmri.metaPar() fits the configured linear mixed-effects meta-analytic (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. Voxel-by-voxel analysis is performed by either the function apply or parApply from parallel package, which walks through the Cbold array.

This two-stage approach reduces the computational burden of fitting a full linear mixed-effects (LME) model, fmri.lmePar would do. It assumes first level design is same across subjects and normally distributed not necessarily homogeneous within-subject errors. Warping to standard space has been done before first-stage 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 meta-analytic 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 intra-subject variability. The REML estimate of cross-subject variability (tau-squared) assumes that each of these observations is drawn independently from the same Gaussian distribution. Since correlation structures cannot be modeled, multi-subject 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.

Value

An object of class "fmrispm" and "fmridata", basically a list with components:

Note

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 between-subject variance (tau-squared) 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 one-stage 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 within-subject (repeated measures) factor

• paired samples with varying tasks, unless the contrast of the two conditions is used as input

Author(s)

Sibylle Dames

References

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: 747-765.

Knapp G. and Hartung J. (2003). Improved tests for a random effects meta-regression with a single covariate. Statistics in Medicine, 22: 2693-2710.

Viechtbauer W. (2005). Bias and efficiency of meta-analytic variance estimators in the random-effects model. Journal of Educational and Behavioral Statistics, 30: 261-293.

Viechtbauer W. (2010). Conducting meta-analyses in R with the metafor package. Journal of Statistical Software, 36(3): 1-48

Viechtbauer W. (2015). metafor: Meta-Analysis Package for R R package version 1.9-7.

See Also

rma.uni, fmri.lm, fmri.lmePar

Examples

## 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: single-session 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: Random-effects meta-regression 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)

P-values

Description

Determine p-values.

Usage

fmri.pvalue(spm, mode="basic", na.rm=FALSE, minimum.signal = 0, alpha= 0.05)

Arguments

Details

If only a contrast is given in spm, we simply use a t-statistic and define p-values according to random field theory for the resulting gaussian field (sufficiently large number of df - see ref.). If spm is a vector of length larger than one for each voxel, a chisq field is calculated and evaluated (see Worsley and Taylor (2006)). If delta is given, a cone statistics is used.

The parameter mode allows for different kinds of p-value calculation. mode="voxelwise" refers to voxelwise tests while mode="Bonferroni" adjusts the significance level for multiple testing. An alternative is mode="FDR" specifying signal detection by False Discovery Rate (FDR) with proportion of false positives level specified by alpha. The other choices apply results on excursion sets of random fields (Worsley 1994, Adler 2003) for smoothed SPM's. "basic" corresponds to a global definition of the resel counts based on the amount of smoothness achieved by an equivalent Gaussian filter. The propagation condition ensures, that under the hypothesis

\hat{\Theta} = 0

adaptive smoothing performs like a non adaptive filter with the same kernel function which justifies this approach. "local" corresponds to a more conservative setting, where the p-value is derived from the estimated local resel counts that has been achieved by adaptive smoothing. In contrast to "basic", "global" takes a global median to adjust for the randomness of the weighting scheme generated by adaptive smoothing. "global" and "local" are more conservative than "basic", that is, they generate slightly larger p-values.

Value

Object with class attributes "fmripvalue" and "fmridata"

Note

Unexpected side effects may occur if spm does not meet the requirements, especially if a parameter estimate vector of length greater than 2 through argument vvector in fmri.lm has been produced for every voxel.

Author(s)

Karsten Tabelow tabelow@wias-berlin.de

References

Polzehl, J. and Tabelow, K. (2007) fmri: A Package for Analyzing fmri Data, R News, 7:13-17 .

Tabelow, K., Polzehl, J., Voss, H.U., and Spokoiny, V. (2006). Analysing fMRI experiments with structure adaptive smoothing procedures, NeuroImage, 33:55-62.

Worsley, K.J., and Taylor, J.E., Detecting fMRI activation allowing for unknown latency of the hemodynamic response, NeuroImage 29:649-654 (2006).

See Also

fmri.lm, fmri.smooth, plot.fmridata, fmri.cluster, fmri.searchlight

Examples

## Not run: fmri.pvalue(smoothresult)

Spacial ICA for fmri data

Description

Uses fastICA to perform spatial ICA on fMRI data.

Usage

fmri.sICA(data, mask=NULL, ncomp=20,
  alg.typ=c("parallel","deflation"), fun=c("logcosh","exp"),
  alpha=1, detrend=TRUE, degree=2, nuisance= NULL, ssmooth=TRUE,
  tsmooth=TRUE, bwt=4, bws=8, unit=c("FWHM","SD"))

Arguments

Details

If specified polynomial trends and effects due to nuisance parameters, e.g., motion parameters, are removed. If smooth==TRUE the resulting residual series is spatially smoothed using a Gaussian kernel with specified bandwidth. ICA components are the estimated using fastICA based on data within brain mask. The components of the result are related as XKW=scomp[mask,] and X=scomp[mask,]*A.

Value

object of class ”fmriICA” list with components

Author(s)

Joerg Polzehl polzehl@wias-berlin.de

See Also

plot.fmriICA,ICAfingerprint, fastICA

Searchlight signal detection

Description

Detection of activated regions using searchlights.

Usage

fmri.searchlight(spm, alpha = 0.05, radius, minimum.signal = 0,
   kind = c("abs", "squared"))

Arguments

Details

The function computes mean statistics (depending on kind) over a searchlight region of radius radius. Approximate voxelwise p-values are determined with respect an empirical (simulated) distribution of the searchlight statistics under the null hypothesis a central t-distributed spm. Thresholding used FDR with rate alpha.

Value

Object with class attributes "fmripvalue" and "fmridata"

Author(s)

Joerg Polzehl polzehl@wias-berlin.de

References

Kriegeskorte, N.; Goebel, R. & Bandettini, P. (2006) Information-based functional brain mapping, PNAS 103:3863-3868.

See Also

fmri.lm, fmri.pvalue, fmri.cluster

Examples

  ## Not run: fmri.searchlight(fmrispmobj)

Spatial group ICA for fmri

Description

Combine ICA results from multiple runs or multiple subjects in group fMRI studies

Usage

fmri.sgroupICA(icaobjlist, thresh = 0.75, minsize=2)

Arguments

Details

The fMRI time series need to be preprocessed and registered before thr ICA decomposition is performed.

The function employs a hierarchical clustering algorithm (complete linkage) on the combined set of spatial independent components obtained from the individual time series. A distance matrix is obtained from correlations of the independent component images. Aggregation of two components from the same fmri series is prevented in the algorithm.

Value

An object of class ”fmrigroupICA” with components

Author(s)

Joerg Polzehl polzehl@wias-berlin.de

References

F. Esposito et al (2005) Independent component analysis of fMRI group studies by self-organizing clustering, Neuroimage, pp. 193-205.

See Also

fmri.sICA, plot.fmrigroupICA, hclust

Smoothing Statistical Parametric Maps

Description

Perform the adaptive weights smoothing procedure

Usage

fmri.smooth(spm, hmax = 4, adaptation="aws",
            lkern="Gaussian", skern="Plateau", weighted=TRUE,...)

Arguments

Details

This function performs the smoothing on the Statistical Parametric Map spm.

hmax is the (maximal) bandwidth used in the last iteration. Choose adaptation as "none" for non adaptive smoothing. lkern can be used for specifying the localization kernel. For comparison with non adaptive methods use "Gaussian" (hmax times the voxelsize in x-direction will give the FWHM bandwidth in mm), for better adaptation use "Plateau" or "Triangle" (default, hmax given in voxel). For lkern="Plateau" and lkern="Triangle" thresholds may be inaccurate, due to a violation of the Gaussian random field assumption under homogeneity. lkern="Plateau" is expected to provide best results with adaptive smoothing.

skern can be used for specifying the kernel for the statistical penalty. "Plateau" is expected to provide the best results, due to a less random weighting scheme.

The function handles zero variances by assigning a large value (1e20) to these variances. Smoothing is restricted to voxel with spm$mask.

Value

object with class attributes "fmrispm" and "fmridata", or "fmrisegment" and "fmridata" for segmentation choice

Author(s)

Joerg Polzehl polzehl@wias-berlin.de, Karsten Tabelow tabelow@wias-berlin.de

References

Polzehl, J., Voss, H.U., and Tabelow, K. (2010). Structural Adaptive Segmentation for Statistical Parametric Mapping, NeuroImage, 52:515-523.

Tabelow, K., Polzehl, J., Voss, H.U., and Spokoiny, V. (2006). Analysing fMRI experiments with structure adaptive smoothing procedures, NeuroImage, 33:55-62.

Polzehl, J. and Spokoiny, V. (2006). Propagation-Separation Approach for Local Likelihood Estimation, Probab. Theory Relat. Fields 135:335-362.

Polzehl, J. and Tabelow, K. (2007) fmri: A Package for Analyzing fmri Data, R News, 7:13-17 .

Examples

## Not run: fmri.smooth(spm, hmax = 4, lkern = "Gaussian")

Linear Model for FMRI Data

Description

Create the expected BOLD response for a given task indicator function.

Usage

  fmri.stimulus(scans = 1, onsets = c(1), durations = c(1), TR = 2,
                times = FALSE, sliceorder = NULL,
                type = c("canonical", "gamma", "boxcar", "user"),
                par = NULL, scale = 10, hrf = NULL, verbose = FALSE)

Arguments

Details

The functions calculates the expected BOLD response for the task indicator function given by the argument as a convolution with the hemodynamic response function.

If sliceorder provides an ordering of slice acquisitions a matrix of expected Bold responses with columns corresponding to the slice number is computed.

For type is "canonical" the latter is modelled by the difference between two gamma functions as given in the reference (with the defaults for a1, a2, b1, b2, cc given therein):

\left(\frac{t}{d_1}\right)^{a_1} \exp \left(-\frac{t-d_1}{b_1}\right) - c \left(\frac{t}{d_2}\right)^{a_2} \exp \left(-\frac{t-d_2}{b_2}\right)

The parameters a1, a2, b1, b2, cc of this function can be changed through the argument par in this order.

Other choices are a simple gamma function

\frac{1}{k\tau_h (k-1)!} \left( \frac{t}{\tau_h} \right)^k \exp \left( - \frac{t}{\tau_h} \right)

or the "boxcar" stimulus, or a user defined function hrf.

The dimension of the function value is set to c(scans, 1).

If !is.null(times) durations are specified in seconds.

Value

Vector with dimension c(scans, 1) or a matrix with dimension c(scans, number of slices).

Author(s)

Karsten Tabelow tabelow@wias-berlin.de, Joerg Polzehl polzehl@wias-berlin.de

References

Worsley, K.J., Liao, C., Aston, J., Petre, V., Duncan, G.H., Morales, F., Evans, A.C. (2002). A general statistical analysis for fMRI data. NeuroImage, 15:1-15.

Polzehl, J. and Tabelow, K. (2007) fmri: A Package for Analyzing fmri Data, R News, 7:13-17 .

See Also

fmri.design, fmri.lm

Examples

  # Example 1
  hrf <- fmri.stimulus(107, c(18, 48, 78), 15, 2)
  z <- fmri.design(hrf, 2)
  par(mfrow=c(2, 2))
  for (i in 1:4) plot(z[, i], type="l")

Generate fmridata example

Description

Generate fmridata example

Usage

gen_fmridata(signal = 1.5, noise = 20, arfactor = 0.3)

Arguments

Value

Object of class fmridata

Examples

gen_fmridata()

Extract searchlight pattern from a SPM

Description

For a provided spm object and a mask of voxel the function extracts the values of the parameter estimates within the searchlight region and for all voxel in the mask.

Usage

getSearchlightPattern(spm, voxel, radius)

Arguments

Value

an array of dimension c(nb, nsl, nvox) with nb the number of estimated parameters in spm$beta, nsl the number of voxel in the searchlight and nvox the number of voxel in the mask provided as second argument

Author(s)

Joerg Polzehl polzehl@wias-berlin.de

See Also

fmri.searchlight, fmri.lm~~~

Translation between smoothness and bandwidth for Gaussian kernel

Description

Translation table between smoothness and bandwidth for Gaussian kernel

Usage

data(hvred)

Format

The format is: num [1:500, 1:2] 0.101 0.102 0.103 0.104 0.105 ...

Examples

data(hvred)
## maybe str(hvred) ; plot(hvred) ...

Create fmridata object from niftiImage

Description

Transforms a niftiImage (created by readNifti from package RNiftyReg) into an object with class fmridata

Usage

niftiImage2fmri(niftiobj, level = 0.75, mask=NULL, setmask = TRUE, indx = NULL,
   indy = NULL, indz = NULL, avoidnegs = FALSE)

Arguments

Details

This function can be used in connection with readNifti from package RNiftyReg to read large fMRI series from nifti files. The resulting fmridata-object stores the image data as 2 byte integer in raw format, in contrast for the 4 byte real used with other functions.

Value

an object of class fmridata

Author(s)

Joerg Polzehl polzehl@wias-berlin.de

See Also

read.AFNI, read.DICOM, read.ANALYZE, read.NIFTI

Diagnostics plots for objects of class ”fmriICA”

Description

The function generates plots for inspecting independent components obtained by spatial independent component analysis.

Usage

## S3 method for class 'fmriICA'
plot(x, comp = 1, center = NULL, thresh = 1.5, ...)
## S3 method for class 'fmrigroupICA'
plot(x, comp = 1, center = NULL, thresh = 1.5, ...)

Arguments

Details

The function generates diagnostic plots for the independent component specified in comp. It provides axial, coronal and sagittal images as determined by center. Values exceeding the threshold are displayed using a color scale. An IC fingerprint is given as a star plot. Additionally the time series corresponding to the spatial IC and its spectral density are plotted.

Value

nothing returned.

Author(s)

Joerg Polzehl polzehl@wias-berlin.de

References

De Martino et. al., Classification of fMRI independent components using IC-fingerprints and support vector machine classifiers, Neuroimage 34 (2007) 177-194.

See Also

fmri.sICA, ICAfingerprint, fastICA

I/O functions

Description

Visualize fMRI data and (intermediate) results.

Usage

## S3 method for class 'fmridata'
plot(x, anatomic = NULL, maxpvalue = 0.05,
              spm = TRUE, pos = c(-1, -1, -1), type = "slice",
              slice =  1, view = "axial" ,zlim.u =
              NULL, zlim.o = NULL,col.o = heat.colors(256), col.u =
              grey(0:255/255), cutOff = c(0, 1), ...)
## S3 method for class 'fmrisegment'
plot(x, anatomic = NULL,
              slice =  1, view = c( "axial", "coronal", "sagittal") ,zlim.u =
              NULL, zlim.o = NULL,col.o = c( rainbow( 64, start = 2/6, end = 4/6),
              rainbow( 64, start = 0, end = 1/6)),
              col.u = grey(0:127/127), verbose = FALSE, ...)

Arguments

Details

Provides a slicewise view of "fmridata" objects with anatomic overlay (if appropriate, that is for class "fmripvalue"). For objects of class "fmrispm" it plots the t-statistics for the estimated effects if spm is TRUE, or the estimated effect otherwise. For objects of class "fmridata" only a plot of the data slices itself is produced. If device is specified as "png", "jpeg", "ppm" output is done to a file. A grey/color scale is provided in the remaining space.

For objects of class "fmrisegment" the smoothed signal size is shown in the activation segments (two-sided test!).

If type is "3d" a 3 dimensional interactive view opens. Sliders to move in the data cube are given ("x", "y", "z", and "t" if class is "fmridata" only). Time series are shown if available. For objects of class "fmrispm" a slider is created to remove information for voxels with smaller signals than a cut-off value from the plot. Use pvalues for statistical evaluation. If spm is FALSE the estimated BOLD response together with a confidence interval corresponding to maxpvalue is drawn. For objects of class "fmripvalue" the pvalues with overlay are shown.

Value

If 'type' is "3d" the Tk-object is returned. (Remove the display with tkdestroy(object))

Note

3 dimensional plotting requires the tkrplot package.

Author(s)

Karsten Tabelow tabelow@wias-berlin.de

References

Polzehl, J. and Tabelow, K. (2007) fmri: A Package for Analyzing fmri Data, R News, 7:13-17 .

See Also

fmri.pvalue

Examples

## Not run: plot(pvalue)

Visualize fMRI p-value maps

Description

Visualize objects created by function fmri.pvalue

Usage

## S3 method for class 'fmripvalue'
plot(x, template = NULL, mask = NULL,
      view = c("axial", "coronal", "sagittal", "orthographic"),
      slices = NULL, ncol = 1, nrow = 1, center = NULL, ...)

Arguments

Value

list with components

Author(s)

Joerg Polzehl polzehl@wias-berlin.de

See Also

fmri.pvalue, ~~~

I/O functions

Description

'print' method for class '"fmridata"'.

Usage

## S3 method for class 'fmridata'
print(x, ...)

Arguments

Details

The method tries to print information on data, like data dimension, voxel size, value range.

Value

none

Author(s)

Karsten Tabelow tabelow@wias-berlin.de

References

Polzehl, J. and Tabelow, K. (2007) fmri: A Package for Analyzing fmri Data, R News, 7:13-17 .

See Also

summary.fmridata

Examples

## Not run: print(data)

I/O function

Description

Read HEAD/BRIK file.

Usage

read.AFNI(filename,vol=NULL,level=0.75,mask=NULL,setmask=TRUE)

Arguments

Details

The function reads a HEAD/BRIK file. If vol is given (defaults to NULL), only volumes in this vector are read, in order to save memory.

Value

Object of class "fmridata" with the following list entries:

Author(s)

Karsten Tabelow tabelow@wias-berlin.de

References

R. W. Cox (1996). AFNI: Software for analysis and visualization of functional magnetic resonance neuroimages. Computers and Biomed. Res. 29:162-173.

Polzehl, J. and Tabelow, K. (2007) fmri: A Package for Analyzing fmri Data, R News, 7:13-17 .

See Also

write.AFNI, read.ANALYZE

Examples

## Not run: afni <- read.AFNI("afnifile")

I/O Functions

Description

Read fMRI data from ANALYZE file(s).

Usage

read.ANALYZE(prefix = "", numbered = FALSE, postfix = "",
             picstart = 1, numbpic = 1, level = 0.75, mask=NULL, setmask=TRUE)

Arguments

Details

This function reads fMRI data files in ANALYZE format. If numbered is FALSE, only the vector of strings in prefix is used for file name (default).

If numbered is TRUE, it takes the first string in prefix and postfix and a number of the form "007" in between to create the file name.

The number is assumed to be 3 digits (including leading zeros). First number is given in picstart, while numbpic defines the total number of images to be read. Data in multiple files will be combined into a four dimensional datacube.

Value

Object of class "fmridata" with the following list entries:

Note

Since numbering and naming of ANALYZE files widely vary, this function may not meet your personal needs. See Details section above for a description.

Author(s)

Karsten Tabelow tabelow@wias-berlin.de

References

Biomedical Imaging Resource (2001). Analyze Program. Mayo Foundation.

Polzehl, J. and Tabelow, K. (2007) fmri: A Package for Analyzing fmri Data, R News, 7:13-17 .

See Also

write.ANALYZE, read.AFNI

Examples

 ## Not run: analyze <- read.ANALYZE("analyze",TRUE,"file",31,107)

I/O function

Description

Read DICOM file.

Usage

read.DICOM(filename,includedata = TRUE)

Arguments

Details

The function reads a DICOM file.

Value

Object with the following list entries:

Note

Since the DICOM standard is rather complicated, there may be cases where this function cannot read a DICOM file. Known issue: it cannot read header with implicit VR. Return value may change in future version!

Author(s)

Karsten Tabelow tabelow@wias-berlin.de

References

http://medical.nema.org

Polzehl, J. and Tabelow, K. (2007) fmri: A Package for Analyzing fmri Data, R News, 7:13-17 .

See Also

read.AFNI, read.ANALYZE

Examples

## Not run: dicom <- read.DICOM("dicomfile")

I/O Functions

Description

Read fMRI data from NIFTI file(s).

Usage

read.NIFTI(filename, level = 0.75, mask=NULL, setmask=TRUE)

Arguments

Details

This function reads fMRI data files in NIfTI format.

The filename can be given with or without extension. If extension is not included, the function searches for the ".nii" file and then for the "hdr/img" pair.

Value

Object of class "fmridata" with the following list entries:

Author(s)

Karsten Tabelow tabelow@wias-berlin.de

References

Polzehl, J. and Tabelow, K. (2007) fmri: A Package for Analyzing fmri Data, R News, 7:13-17 .

See Also

read.ANALYZE, read.AFNI

Examples

 ## Not run: analyze <- read.NIFTI("niftifile.nii")

Add or replace mask in an fmridata object

Description

The function replaces the information in the mask component of an fmridata object.

Usage

setmask(fmriobj, mask)

Arguments

Details

The dimensions of both objects supplied as arguments need to be compatible.

Value

on object of class 'fmridata'.

Author(s)

Joerg Polzehl polzehl@wias-berlin.de

See Also

oro2fmri, niftiImage2fmri, read.NIFTI, read.AFNI, read.ANALYZE

A function for sinc-interpolation

Description

Performs sinc interpolation for a equidistant time series x to times t.

Usage

sincfilter(t, x, wr=8)

Arguments

Value

a vector of interpolated values of the time series at time points given in t.

Author(s)

Joerg Polzehl polzehl@wias-berlin.de

See Also

slicetiming

Examples

  x <- 1:107
  y <- rnorm(x)
  z <- sincfilter(seq(1,107,.01),y)
  plot(x, y, ylim=range(y,z))
  lines(seq(1,107,.01),z,col=2)

slicetiming for fmridata-objects

Description

Perform slicetiming for fMRI data, ideally before preprocessing (registration). Recording times for slices are assumed to be equispaced between scans with argument sliceorder providing the order of slice acquisitions. Interpolation between slices is performed using a sinc filter.

Usage

slicetiming(fmridataobj, sliceorder = NULL)

Arguments

Value

an object of class fmridata

Author(s)

Joerg Polzehl polzehl@wias-berlin.de

See Also

fmri.stimulus, fmri.design,fmri.lm,~~~

Examples

  ## Not run: 
# Example 1
  data <- list(ttt=writeBin(rnorm(32*32*32*107), raw(), 4),
               mask=array(TRUE, c(32, 32, 32)), dim=c(32, 32, 32, 107))
  class(data) <- "fmridata"
  data <- slicetiming(data,sliceorder=1:32)
  ## provides data corrected for sequential slice acquisition in linear order

## End(Not run)

I/O functions

Description

'summary' method for class '"fmridata"'.

Usage

## S3 method for class 'fmridata'
summary(object, ...)

Arguments

Details

The method tries to print information on data, like data dimension, voxel size, value range.

Value

A list with the following elements:

Author(s)

Karsten Tabelow tabelow@wias-berlin.de

See Also

print.fmridata

Examples

## Not run: summary(data)

I/O functions

Description

Write BRIK/HEAD files.

Usage

write.AFNI(filename, ttt, label = NULL, note = NULL, origin = NULL, 
               delta = NULL, idcode = NULL, header = NULL, taxis = FALSE)

Arguments

Details

Write out BRIK/HEAD files as required by AFNI. Most arguments correspond to entries in the HEAD file, but use is depreciated. Use header and taxis instead!

Value

Nothing is returned.

Author(s)

Karsten Tabelow tabelow@wias-berlin.de

References

Polzehl, J. and Tabelow, K. (2007) fmri: A Package for Analyzing fmri Data, R News, 7:13-17 .

See Also

read.AFNI,write.ANALYZE

Examples

 ## Not run: write.AFNI(tempfile(), array(as.integer(65526*runif(10*10*10*20)),
     c(10,10,10,20)), c("signal"), note="random data",
     origin=c(0,0,0), delta=c(4,4,5), idcode="unique ID")
## End(Not run)
 write.AFNI(tempfile(), array(as.integer(65526*runif(10*10*10*20)),
     c(10,10,10,20)), header=list(HISTORY_NOTE="random data",
     ORIGIN=c(0,0,0), DELTA=c(4,4,5), IDCODE_STRING="unique ID"),taxis=FALSE)

I/O Functions

Description

Write a 4 dimensional datacube in ANALYZE file format.

Usage

write.ANALYZE(ttt, header=NULL, filename)

Arguments

Details

Writes the datacube ttt to a file named file in ANALYZE file format. header is a list that contains the header information as documented by the Mayo Foundation. We give here a short summary. If a value is not provided, it will be tried to fill it with reasonable defaults, but do not expect fine results, if the entry has a special important meaning (h.i. pixdim).

See ANALYZE documentation for details.

Value

Nothing is returned.

Author(s)

Karsten Tabelow tabelow@wias-berlin.de

References

Polzehl, J. and Tabelow, K. (2007) fmri: A Package for Analyzing fmri Data, R News, 7:13-17 .

See Also

read.ANALYZE, write.AFNI

Examples

## Example 1
write.ANALYZE(array(as.integer(65526*runif(10*10*10*20)),c(10,10,10,20)),
              file=file.path(tempdir(),"analyzefile"))

I/O Functions

Description

Write a 4 dimensional datacube in NIfTI file format.

Usage

write.NIFTI(ttt, header=NULL, filename)

Arguments

Details

Writes the datacube ttt to a file named file in NIfTI file format. header is a list that contains the header information.

See NIfTI documentation for details.

Value

Nothing is returned.

Author(s)

Karsten Tabelow tabelow@wias-berlin.de

References

Polzehl, J. and Tabelow, K. (2007) fmri: A Package for Analyzing fmri Data, R News, 7:13-17 .

See Also

read.ANALYZE, write.AFNI

Examples

## Example 1
write.NIFTI(array(as.integer(65526*runif(10*10*10*20)),c(10,10,10,20)),
              file=file.path(tempdir(),"niftifile"))