Statistical Parametric Maps for Functional MRI Experiments in R : The Package fmri

The purpose of the package fmri is the analysis of single subject functional magnetic resonance imaging (fMRI) data. It provides fMRI analysis from time series modeling by a linear model to signal detection and publication quality images. Speciﬁcally, it implements structural adaptive smoothing methods with signal detection for adaptive noise reduction which avoids blurring of activation areas. Within this paper we describe the complete pipeline for fMRI analysis using fmri . We describe data reading from various medical imaging formats and the linear modeling used to create the statistical parametric maps. We review the rationale behind the structural adaptive smoothing algorithms and explain their usage from the package fmri . We demonstrate the results of such analysis using two experimental datasets. Finally, we report on the usage of a graphical user interface for some of the package functions.


Introduction
Neuroscience is an active field that combines challenging scientific questions with interesting methodological developments.It draws expertise from diverse fields as biology, medicine, physics, mathematics, statistics, and computer science.In fact, the emergence of medical imaging techniques has triggered a boost of developments for medical or scientific applications in particular for the examination of the human brain.One of them is functional magnetic resonance imaging (fMRI, Friston et al. 2007a;Lazar 2008) which, due to its non-invasive character, has become the most informative tool for in-vivo examination of human brain function on small spatial scales.It is utilized both in research as well as clinical applications such as diagnosis and treatment of brain lesions (Vlieger et al. 2004;Haberg et al. 2004;Henson et al. 2005).
As fMRI suffers from a low signal to noise ratio (SNR), noise reduction plays a very important role.Gaussian filtering is applied in standard analysis mostly to increase the SNR and the sensitivity of statistical inference (Friston et al. 2007a).At the same time smoothing reduces the number of independent decisions, relaxes the severe multiple test problem and leads to a situation where critical values for signal detection can be assigned using random field theory (RFT, Worsley 1994;Adler 2000).The inherent blurring effect can be ignored as long as the precise shape and extent of the activation area is not crucial.However, experiments are becoming more and more sophisticated and explore columnar functional structures in the brain (Cheng et al. 2001) or functions near brain lesions.Here, adaptive noise reduction methods are essential.For a comprehensive introduction into statistical issues in fMRI we refer to Lazar (2008).
In the R package fmri (Tabelow and Polzehl 2011) two algorithms from a special class of structural adaptive smoothing methods are implemented together with appropriate signal detection: (a) structural adaptive smoothing with RFT for signal detection and (b) structural adaptive segmentation based on multiscale tests.The package fmri (version 1.0-0) was first reported in Polzehl and Tabelow (2007b).Since then, the package has evolved significantly and now includes additional features, advances in methodology, and a graphical user interface (GUI).Here, we will report on these updates which refer to the package version 1.4-6.For an introduction to the R environment for statistical computing (R Development Core Team 2010) see Dalgaard (2008) or the materials at http://www.R-project.org/.For the use of R in medical imaging problems we refer the reader to the CRAN task view on Medical Imaging (Whitcher 2010a).

Data handling and linear model
The package fmri analyzes fMRI time series using the general linear modeling approach (Friston et al. 2007b).The focus of the package is on the use of structural adaptive smoothing methods (Section 3).These methods increase sensitivity of signal detection while limiting the loss in spatial resolution.For an alternative modeling of fMRI data using independent component analysis (ICA) we refer the reader to the package AnalyzeFMRI (Marchini and Lafaye de Micheaux 2010; Bordier et al. 2011).Other analysis methods for fMRI data are implemented for example in the package arf3ds4 (activated region fitting, Weeda 2010; Weeda et al. 2011) and cudaBayesreg (Bayesian multilevel model, Ferreira da Silva 2010, 2011).

Data formats
fMRI experiments typically acquire time series of full three dimensional brain volumes when a specific stimulus or task is presented or performed, respectively.The package fmri provides functionality to read medical imaging data from several formats.It may occasionally be preferable to convert the data using external tools, e.g., AFNI (Cox 1996), MRIcro (Rorden and Brett 2000), to name only two, or use more specialized R packages, oro.dicom (Whitcher 2010b) or oro.nifti (Whitcher et al. 2010(Whitcher et al. , 2011) ) for data input.Objects containing fMRI data in the format used within fmri may then be easily generated, see documentation and the description below.
The standard format for data recorded by a MR scanner is DICOM (Digital Imaging and Communications in Medicine; http://medical.nema.org/)1 .The package fmri provides a function for reading DICOM files: where <filename> is the name of the file.A data description, without actually reading the data, can be obtained by setting the logical argument includedata to FALSE.The value returned by the function read.DICOM is a list ds containing the full header information as a list with the four-byte sequences (group,element) of the DICOM format as elements, the data (if read), and specific header information, e.g., voxel size, image dimension, which are of special interest.
Other formats are more common in the context of fMRI.The ANALYZE format developed at the Mayo Clinic is comprised of two files, the .imgand .hdrfiles, containing the image data and information about the image acquisition, respectively.Using the fmri package the file <filename> is read by where <filename> contains the data as a 4D volume.Data prepared for the analysis software SPM (Ashburner et al. 2008) generally comes as a series of 3D volume with file names of the form <prefix+3digits+postfix>. It is possible to read this series using where picstart is the first number in the digits series and <number> is the number of volumes.

Data preprocessing
Data from fMRI experiments acquisition need to be prepared ahead of the quantitative analysis.For example, head motion during the experiment leads to mis-registration of voxels within the data cube at different time points.For group analysis the individual brain has to be mapped (normalized) onto a "standard brain" to identify corresponding brain sections.
There are currently (version 1.4-6) no tools for preprocessing steps like motion-correction, registration and normalization within the package fmri.There exist many tools to perform these steps in advance, like SPM, AFNI, FSL (Smith et al. 2004), or others.Note, that there is an R package RNiftyReg to perform image registration (Clayden 2010).fMRI analysis with package fmri is currently restricted to single subject analysis with functions to perform group comparisons planned for future versions.Like all imaging modalities, fMRI suffers from a variety of noise sources, rendering subsequent signal detection difficult.In order to increase sensitivity smoothing of fMRI data is usually part of fMRI preprocessing.It is very important to note, that within the package fmri smoothing should not be applied as preprocessing step.Instead structural adaptive smoothing, which is the main feature of the package, is applied to the statistical parametric map (SPM) derived from the linear model for the time series.If the user intends to use the package fmri it is advisable not to smooth the data with third-party tools.Doing so will lead to a loss of essential information and impose a spatial correlation structure, which reduces the effect of structural adaptive smoothing.

Linear modeling
The observation that voxels with increased neuronal activity are characterized by a higher blood oxygenation level (Ogawa et al. 1990(Ogawa et al. , 1992) ) is known as the BOLD (blood oxygenation level dependent) effect.In fMRI it can be used as a endogenous contrast.Together with the fact that nuclear magnetic resonance is free of high energy radiation, this forms the basis of its appeal as a non-invasive biomarker.
The expected BOLD response corresponding to the stimulus or task of the fMRI experiment may be modeled by a convolution of the experimental design with the hemodynamic response function (HRF).This function characterizes the time delay and modified form of the response in blood oxygenation.Within the package fmri the hemodynamic response function h(t) is modeled as the difference of two gamma functions with default parameters a 1 = 6, a 2 = 12, b 1 = 0.9, b 2 = 0.9, and d i = a i b i , i = 1, 2, c = 0.35 and t the time in seconds (Glover 1999).Given the stimulus s(t) as a task indicator function, the expected BOLD response is calculated as the convolution of s(t) and h(t) The resulting function x(t) is evaluated at T scan acquisition times.Such expected BOLD signals may be created by the function R> hrf <-fmri.stimulus(scans= 1, onsets = 1, durations = 1, rt = 3, + times = NULL, mean = TRUE, a1 = 6, a2 = 12, b1 = 0.9, b2 = 0.9, where hrf is a vector of length T , onsets is the vector of scan numbers of the stimulus onsets, durations is the vector of stimulus durations in number of scans, and rt is the time between subsequent volumes.The design can be given in seconds by providing the argument times, in this case durations should also be given in seconds.One can change the parameters of h(t) to use a slightly different form for the HRF.The vector hrf contains the values of the expected BOLD response at the T scan acquisition times.Note, that slice time correction and local modeling of h(t) are currently not available within fmri.
Figure 1 shows the vector hrf created by fmri.stimulus() using R> hrf <-fmri.stimulus(105, c(16, 46, 76), 15, 2) and used for the analysis of the sports imagination dataset in Section 4.2.Alternatively, a vector of length T (or number of scans) describing an expected BOLD response may be supplied by the user.
For fmri we adopt the common view of a linear model for the time series Y i = (Y it ) in each voxel i after reconstruction of the raw data and motion correction where X denotes the design matrix, containing the expected bold responses hrf1, hrf2, . . ., hrfq corresponding to q experimental stimuli as columns (Friston et al. 1995;Worsley et al. 2002).Additional columns represent a mean effect and a polynomial trend of specified order.Using X we may also model s ≥ 0 confounding effects ce like respiration or heart beat.The columns of X that model the trend are chosen to be orthogonal to the stimulus effects.The components β ij (j = 1, . . ., q) of β i describe the effect of stimulus j for voxel i.The additive errors i = ( i1 , . . .iT ) are assumed to have zero mean and a variance depending on the underlying tissue in voxel i.
The design matrix X may be created using R> x <-fmri.design(cbind(hrf1,..., hrfq, ce1, ...), order = 2) where order is the order of the polynomial trend.Similar to the function hrf, x is simply a p × T matrix (p = q + s + order) and may also be provided directly by the user.The linear model is evaluated for the fMRI data in ds by R> spm <-fmri.lm(ds,x, contrast = 1) In order to access the variability of the parameter estimates βi correctly we need to model the structure of Σ i .In the package fmri we assume a temporal auto-correlation following an auto-regressive (AR) model of order 1.The auto-correlation coefficients ρ i are estimated from the residual vector r i = (r i1 , . . ., r iT ) of the fitted model 1 assuming independence, i.e., Σ i = σ 2 i I T .We apply the bias correction given by Worsley et al. (2002).Temporal correlation is usually assumed to vary smoothly in space and, by default, we smooth the auto-correlation coefficients using a Gaussian filter in order to stabilize the estimates and increase the number of degrees of freedom (Worsley 2005).Using the resulting estimates ρi we define a transformation matrix A i = R −1/2 (ρ i ), where R(ρ) denotes the correlation matrix under the assumption of an AR(1) process with parameter ρ for voxel i (Tabelow et al. 2006).
The model 1 is then transformed into a linear model where

The error variance σ 2
i is estimated from the residuals ri of the linear model 2 by σ2 i = T 1 r2 it /(T − p) leading to the voxel-wise estimated covariance matrix Usually special interest is in one of the parameters β i or a contrast of parameters γ i = c β i represented by the vector c.Hence the result of the parameter estimation in the linear model consists of two three dimensional arrays Γ and S containing the estimated effects γi = c βi and their estimated standard deviations si = c COV( βi )c.
The function fmri.lmgenerates a list object with class attribute "fmrispm" containing the arrays Γ and S2 as components beta and var, respectively, and optionally residual information in "raw" format, similar to the imaging data.One can use the extract.datafunction with argument what = "residuals" to extract the residuals from the object.Note, that changing the contrast of interest c requires re-estimation of the linear model using a new value for the argument contrast of fmri.lm.
The voxel-wise quotient θ i = γi /s i of both arrays forms a statistical parametric map Θ which is approximately a random t-field (Worsley 1994).Note that the number of degrees of freedom within the t-field depends on the temporal auto-correlation and the smoothing of the AR(1) coefficients (Worsley 2005).All arrays carry a correlation structure induced by the spatial correlation in the fMRI data.
The parameter estimates Θ may now be used for a voxel-wise signal detection.A voxel i is classified as activated if the corresponding value θi exceeds a critical value.In this case the estimated BOLD signal (or contrast) significantly deviates from zero.There are typically two intrinsic problems with signal detection.The first is related to the large number of multiple tests in the data cube leading to a large number of false positives if the significance level is not adjusted for multiple testing.Possible solutions are Bonferroni corrections, assuming independence, thresholds obtained from random field theory (Worsley 1994;Adler 2000), depending on sufficient spatial smoothness, or the false discovery rate (FDR, Benjamini and Hochberg 1995;Benjamini and Heller 2007).The second problem is related to small values in Θ caused by a large error variance.For both problems smoothing of Γ helps, as it reduces the number of independent tests and ideally reduces the variance of the estimated parameters or contrast.In the next section we will review two structural adaptive smoothing methods which are at the core of the package fmri.

Structural adaptive smoothing
The most common smoothing method for functional MRI data is the Gaussian filter in 3D.It may be easily applied using the fast Fourier transform and guarantees the assumption for RFT which requires a certain smoothness in the data.Unfortunately, non-adaptive smoothing leads to significant blurring and reduces specificity with respect to the spatial extent and form of the activation areas in fMRI.There are several applications where the blurring renders subsequent medical decisions or neuroscientific analysis difficult to interpret; e.g., pre-surgical planning for brain tumor resection (Vlieger et al. 2004;Haberg et al. 2004;Henson et al. 2005) or columnar functional structures (Cheng et al. 2001).We have proposed structural adaptive smoothing methods that overcome these drawbacks and allow for increased sensitivity of signal detection while preserving the borders of the activation areas (Tabelow et al. 2006;Polzehl et al. 2010).Smoothing is generally considered a preprocessing step for the data as is motion-correction etc.However, except for effects from pre-whitening, the order in which non-adaptive spatial smoothing and evaluation of the linear model are performed is arbitrary.Moreover, if the temporal correlations are spatially homogeneous the temporal modeling and spatial smoothing may be interchanged.
Structural adaptive smoothing uses local homogeneity tests for the adaptation (Polzehl and Spokoiny 2006) that are inefficient in the four-dimensional data space.Structural adaptive smoothing is therefore based on the estimates of Γ and S obtained by temporal modeling using fmri.lm.Parameter estimation in the linear model serves as a variance and dimension reduction step prior to spatial smoothing and therefore allows for a much better adaptation (Tabelow et al. 2006).
Structural adaptive smoothing requires an assumption on the spatial homogeneity structure of the true parameter γ i .We assume a local constant function for γ i .The null hypothesis for the statistical test of activation is H 0 : γ i = 0, i.e., we therefore assume that non-activated areas are characterized by a parameter value of zero.In areas which are activated during the acquisition the parameter values differ from zero and are similar, provided that the BOLD %-changes are similar.
Based on this assumption and the arrays Γ and S, we use an iterative smoothing algorithm for the SPM that is based on pairwise tests of homogeneity.We specify kernel functions For K l this choice is motivated by the kernels near efficiency in non-adaptive smoothing.Its compact support and simplicity reduces the computational effort.The second kernel should exhibit a plateau near zero, be compactly supported and monotone non-increasing on the positive axis.We denote by d(i, j) the distance between voxel i and j.Let H = {h k } k * k=1 , with h 0 = 1, be a series of bandwidths such that forms a geometric sequence with factor c h = 1.25.This provides, from step to step and in case of non-adaptive smoothing, a variance reduction by factor c h .At iteration step k we define for each voxel i and all voxels j within a distance d(i, j) < h k weights We call s The parameter λ in (3) can be determined by simulation depending on the error model (but not the actual data) using a propagation condition (Polzehl and Spokoiny 2006).
The term C(k, ρ s ) in ( 3) provides an adjustment for the effect of spatial correlation in the original data, characterized by a vector of first-order correlation coefficients ρ s , at iteration k.Optionally, specifying adaptation = "fullaws" for the function fmri.smooth,causes the term to be replaced by an estimated variance σ2;(k−1) . This leads to additional computational cost but usually improves the results.Finally, the variances of the smoothed estimates (contrasts) are estimated from smoothed residuals using the weighting scheme w Structural adaptive smoothing provides an intrinsic stopping criterion, since the quality of estimates is preserved at later steps of the algorithm.It leads to estimates Γ(k * ) and S2;(k * ) and a smoothed SPM Θ(k * ) .In these statistics shape and borders of the activation structure are preserved.As a consequence, compared to non-adaptive smoothing methods, the procedure reduces noise while preserving the resolution of the scan as required by many modern applications (Tabelow et al. 2009).
The number of iteration steps k * determines a maximum achievable variance reduction or equivalently a maximal achievable smoothness.It can be specified by selecting the maximum bandwidth h max in the series h k as the expected diameter of the largest area of activation.Over-smoothing structural borders is avoided in the algorithm as long as differences between the parameter values are statistically significant at any scale, i.e., for any bandwidth h k , visited within the iterations.The largest homogeneous region is expected to be the non-activation area, where parameter values do not significantly differ from zero.Within this area, under the hypothesis of no activation, structural adaptive smoothing due to the propagation condition behaves like its non-adaptive counterpart with λ = ∞.
Structural adaptive smoothing within the package fmri is performed by > spm.smooth <-fmri.smooth(spm,hmax = 4, adaptation = "aws") The choice of the adaptation method is specified using the argument adaptation.Specifically adaptation = "none" uses non-adaptive smoothing, adaptation = "aws" or "fullaws" corresponds to adaptive smoothing (with improved variance estimates for the latter at the cost of computation time), and adaptation = "segment" uses the structural adaptive segmentation described in the next section.Note, hmax is the only other parameter for the algorithm that should be specified by the user.Additional parameters of the function have only minor influence on the results and should be considered for experts usage and testing purpose only.
From the final estimates for k = k the random t-field Θ(k ) = Γ (k ) ( S2;(k * ) ) −1/2 can be constructed.Under the null hypotheses of no activation the propagation condition ensures a smoothness corresponding to the application of a non-adaptive filter with the bandwidth h max = h k .In order to define thresholds for signal detection on such a smooth t-field RFT is applied in the package fmri (Adler 2000;Tabelow et al. 2006).
Signal detection is performed by function fmri.pvalueapplied to an object of class "fmrispm".
R> pvalue <-fmri.pvalue(spm.smooth)creates a list object of class "fmripvalue" containing the p values at each voxel in the list element pvalues.For the determination of the p values we do not include effects from the nonstationarity of the random field (Hayasaka et al. 2004).This could be done using the objects spm and spm.smooth but is currently not implemented in fmri.pvalue.If no smoothing has been performed then signal detection by FDR may be specified using the argument mode = "FDR" in function fmri.pvalue.The result pvalue is displayed via R> plot(pvalue, anatomic = anatomic, maxpvalue = 0.05) with a possible anatomic image anatomic used as an underlay.The anatomic image must be supplied in NIfTI format using read.NIFTI or as an array of the same size as the functional data.The significance level for the multiple test corrected p values can be chosen by an appropriate value for maxpvalue.

Structural adaptive segmentation
Structural adaptive smoothing and subsequent signal detection by RFT forms a sequential procedure that depends on the similar behavior of the adaptive method and non-adaptive smoothing under the null hypothesis of no activation.It is desirable to combine adaptive smoothing and signal detection in a way that solves the multiple-comparison problem.Such integration is possible since the information used to generate weighting schemes in the structural adaptive procedure (Tabelow et al. 2006) can also be used for signal detection.The resulting method is called structural adaptive segmentation (Polzehl et al. 2010).This method leads to similar signal detection results, but is conceptually more coherent and requires fewer assumptions concerning the signal detection.For example, it automatically accounts for the spatial smoothness in the data and does not rely on the assumption of stationary errors.It also provides a more computationally efficient algorithm.
Let V be a region-of-interest (ROI), and construct a hypothesis test where δ corresponds to a suitable minimal signal size for γ.The probability to reject the hypothesis in any voxel i ∈ V should be less or equal a prescribed significance level α.
Using results from extreme value theory (Resnick 1987;Polzehl et al. 2010) and ideas from multiscale testing (Dümbgen and Spokoiny 2001) a suitable test statistics for this hypotheses is where the degrees of freedom ν of the t statistics θi is, in the non-adaptive case, the shape parameter of the limiting Fréchet distribution (Resnick 1987).The normalizing sequences a n(h) (ν) and b n(h) (ν) are selected to achieve, under the hypothesis H, δ = 0 and λ = ∞, a good approximation of the distribution of T ( ΓH ) by the Fréchet distribution Φ ν .
The adaptive segmentation algorithm employs in each iteration step a test statistic where n i (h) ≤ n(h) reflects the effect of adaptive weights, to decide whether or not to reject the null hypothesis at iteration k.Adaptive smoothing is performed as long as the hypothesis is not rejected.Otherwise non-adaptive smoothing is restricted to the set of detected deviance's from the hypothesis in either a positive or negative direction.Critical values are determined by simulation under the null hypothesis as quantiles of the distribution of for a wide range of n, ν and suitable values of λ.A minimum value of λ has again been selected to fulfill the propagation condition under the null hypothesis (Polzehl et al. 2010).
Structural adaptive segmentation combines the structural adaptive smoothing algorithm described in the preceding section with a test using T ( Γ) (k) at iteration k.The result provides three segments, one, where the hypotheses of no activation could not be rejected, and two segments for the two-sided alternatives.Note, no p values are determined but the method uses a specified significance level α for determining the critical values for the test.There are several parameters for the structural adaptive segmentation (Polzehl et al. 2010).Most, like the choice of kernel function or the rate of increase for the series of bandwidths for the iteration, have only minor influence.These are pre-defined in the package and cannot be changed by the user.The parameter λ, which controls the degree of adaptation, is chosen for any specified error distribution by simulation using the propagation-separation condition (Polzehl et al. 2010), and cannot be changed by the user.
In the fmri package structural adaptive segmentation may be performed on an object fmrispm from the linear model fmri.lmby choosing adaptation = "segment" in R> spm.segment <-fmri.smooth(spm,hmax = 4, adaptation = "segment", + alpha = 0.05) The (multiple test) significance level alpha may be chosen within the range 0.01 to 0.2.Additionally, an argument delta may be provided for the minimum signal size used in the test statistics (5).The result of the structural adaptive segmentation contains information on three segments that can be visualized using fmri: Statistical Parametric Maps for fMRI Experiments in R R> plot(spm.segment) Within activated areas the size of the estimated signal is provided as additional color-coded information.

Examples
4.1.Auditory experiment (dataset 1) In our first example we use an auditory dataset recorded by Geraint Rees under the direction of Karl Friston and the Functional Imaging Laboratory (FIL) methods group which is available from the SPM website (Ashburner et al. 2008).96 acquisitions were made with scan to scan repeat time RT=7s, in blocks of 6, resulting in 16 blocks of 42s duration.The condition for successive blocks alternated between rest and auditory stimulation, starting with a resting block.Auditory stimulation was with bi-syllabic words presented binaurally at a rate of 60 per minute.The functional data starts at acquisition 4, image fM00223_004.Echo planar images (EPI) were acquired on a modified 2T Siemens MAGNETOM Vision system.Each acquisition consisted of 64 contiguous slices with matrix size 64 × 64 and 27mm 3 isotropic voxel (http://www.fil.ion.ucl.ac.uk/spm/data/auditory/).
The following script (spm.R) can be used to process the first dataset.As an alternative the package provides a graphical user interface (GUI), see Section 5, to guide the user through the analysis.

Sports imagination experiment (dataset 2)
For the second example a sports imagination task fMRI scan was performed by one healthy adult female subject within a research protocol approved by the institutional review board of Weill Cornell Medical College.For functional MRI, a GE-EPI sequence with TE/TR = 40/2000 ms was used and 30 axial slices of 4 mm thickness were acquired on a 3T GE system.A field-of-view (FOV) of 24 cm with a matrix size 64 × 64, yielding in-plane voxel dimensions of 3.75 mm, respectively, was used.The excitation flip angle was 80 degrees.Task and rest blocks had a duration of 30 s and were played out in the following order: rest, task, rest, task, rest, task, rest, totaling 105 repetitions.Before the first block, 6 dummy scans where acquired to allow for saturation equilibrium.The task consisted of imagination of playing tennis.The following script (imagination.R) can be used to process the second dataset.

Results
The results of the script spm.R run on the first dataset can be seen in Figure 2. Some slices (33,34,36,37) with detected auditory activation have been chosen in the GUI invoked by the plot function in the script.The data has been smoothed using the structural adaptive segmentation algorithm with a maximum bandwidth of four voxels and a significance level of 0.05 (default).Figure 2 shows that there does not occur an apparent blurring effect as it would after applying a Gaussian filter with corresponding bandwidth.The color in the activation area codes the size of the estimated BOLD signal.Images obtained using the script imagination.R and the second dataset are shown in Figure 3.
The figure shows all slices of the dataset with detected signals using structural adaptive segmentation with a maximum bandwidth of four voxels and a significance level of 0.05.Activations are found primarily in the supplementary motor area and the premotor areas.
The structural adaptive smoothing leads to increased sensitivity without blurring the borders of the activation areas.The color again codes the size of the estimated BOLD signal after adaptive smoothing in arbitrary units using the full range of the color from red to white for the range of estimated signals.Note, that the MP-RAGE scan is not mapped to the EPI data.Therefore good matching between functional and anatomical data is limited to the upper brain slices with the motor activation, which is the main result of this experiment.Some spurious activation outside the brain is due to a pronounced EPI ghost signal as a result of the smoothing algorithm.

Organizing the computational work flow (GUI)
In order to provide a user friendly environment we created a graphical user interface (GUI) which guides the user through the work flow described above.The GUI does not provide the full functionality of the package but gives quick results for standard analysis scenarios.
After invoking the package fmri the GUI can be started R> library("fmri") R> fmrigui() and provides step by step instructions.At each step of the analysis help information is available.The following design information needs to be provided to perform the analysis for the auditory experiment.The information may be saved in a text file and re-used.Data may be accessed in AFNI, AN-ALYZE, or NIfTI files.In this example, use the file select box to navigate to the directory containing the ANALYZE files and select the first file within the directory, file fM00223_004.hdr in folder fM00223.The data are read and basic consistency checks against the experimental design are performed.
In the next step a threshold for mean image intensity is proposed.This threshold is used to define a mask of voxel with mean intensity larger than the threshold.The mask should contain voxel within the brain so that computations are restricted to voxel within this mask.The button View Mask provides images of the mask defined by the specified threshold together with density plots of image intensities for centered data cubes of varying sizes.This information is provided to assist the selection of an appropriate threshold.Next a contrast must be specified as a vector separated by commas or blanks.Trailing zeros may be omitted.Given this information the SPM and the corresponding variance estimates are computed following the steps described in Section 2.3 using the function fmri.lm.
Finally a significance level (default: 0.05) needs to be specified.This is used within the adaptive segmentation algorithm (Section 3.2) or to specify the threshold for multiple test correction using RFT after structural adaptive smoothing.Finally adaptive smoothing (Section 3.1) or segmentation (Section 3.2) using a specified bandwidth may be performed.and coronal slices with activations.The 2D view has more features (Figure 5).Sliders are available to navigate between the slices, while the view may be changed from axial to coronal, or sagittal.For convenience only four slices are shown by default.This may be changed using the number of slices and number of slices per page fields (Change View/Slices to apply the changes).
The Extract images button allows one to export images in JPEG or PNG format.Image processing is done using the package adimpro (Polzehl and Tabelow 2007a).

Conclusion
The R package fmri is a package for the analysis of single-subject fMRI data.Group analysis is currently not supported but will be included in a later version.
The package uses a linear model for the BOLD response.Based on the statistical parametric map noise reduction is performed by structural adaptive smoothing algorithms that reduce noise while not blurring the shape and extent of activated areas (Tabelow et al. 2006;Polzehl et al. 2010).
We demonstrated the usage of the main functions for fMRI analysis with two datasets which are publicly available.Finally we reported on the GUI which organizes the work-flow with slightly reduced functionality compared to the command-line.
Further development will include group analysis and a re-implementation using the advantages of S4 classes (Chambers 2008).One of the great strengths of the package is that it implements the complete fMRI analysis pipeline based on the general linear model (GLM) approach using new smoothing methods including signal detection.Other R packages for fMRI analysis are for example AnalyzeFMRI and arf3DS4.

Figure 1 :
Figure 1: Expected BOLD response for the sports imagination data set.

1Figure 2 :
Figure 2: Result of signal detection using structural adaptive segmentation and the dataset 1 (auditory experiment, slices 33, 34, 36, and 37).The color scheme codes the size of the estimated signal in a voxel where a activated segment has been detected by the algorithm.The underlay is the first volume of the time series.

Figure 3 :
Figure 3: Result of signal detection using structural adaptive segmentation and the second dataset described above.The color scheme codes the size of the estimated signal in a voxel where an activated segment has been detected by the algorithm.The underlay is a sagittal 3D MP-RAGE scan.

Figure 4 :
Figure 4: Snapshot of the fmri GUI after performing all steps.
Figure  4provides a snapshot of the GUI.Note that, if needed the GUI may be closed with or without saving the current results into the global environment of the R session.This enables one to continue the analysis using the full functionality of the package from the console.In a last step results are visualized.The results view in 3D mode shows axial, sagittal

Figure 5 :
Figure 5: Snapshot of the results window of the GUI.This is the same as invoked by the plot() function.
The input functions for ANALYZE, NIfTI, and AFNI files possess additional arguments setmask and level.These can be used to classify voxel within and outside the brain based on the mean voxel intensity.If setmask is TRUE (default) a threshold is set as the levelquantile of the mean voxel intensities.Voxel with mean intensity exceeding this threshold are considered to be within the brain.Further computations are restricted to such voxel.
The ANALYZE format has been further adapted to the NIfTI format by the DFWG (Data Format Working Group, http://nifti.nimh.nih.gov/Members/dfwg/) with improvements in extending the header information and the possibility to merge both files into one (.nii).Reading these files is performed by R> ds <-read.NIFTI(<filename>) Note, that currently (version 1.4-6) compressed NIfTI files cannot be handled by fmri.Last but not least R> ds <-read.AFNI(<filename>) reads file pairs (HEAD/BRIK) used by AFNI.