Affiliations F. M. Kirby Research Center for Functional Brain Imaging, Kennedy Krieger Institute, Baltimore, Maryland, United States of America, Department of Radiology, The Johns Hopkins School of Medicine, Baltimore, Maryland, United States of America ⨯
Affiliations F. M. Kirby Research Center for Functional Brain Imaging, Kennedy Krieger Institute, Baltimore, Maryland, United States of America, Department of Neurology and Psychiatry, The Johns Hopkins School of Medicine, Baltimore, Maryland, United States of America ⨯
Affiliation Department of Biostatistics, Bloomberg School of Public Health, Johns Hopkins University, Baltimore, Maryland, United States of America ⨯
Functional magnetic resonance imaging (fMRI) is a thriving field that plays an important role in medical imaging analysis, biological and neuroscience research and practice. This manuscript gives a didactic introduction to the statistical analysis of fMRI data using the R project, along with the relevant R code. The goal is to give statisticians who would like to pursue research in this area a quick tutorial for programming with fMRI data. References of relevant packages and papers are provided for those interested in more advanced analysis.
Citation: Eloyan A, Li S, Muschelli J, Pekar JJ, Mostofsky SH, Caffo BS (2014) Analytic Programming with fMRI Data: A Quick-Start Guide for Statisticians Using R. PLoS ONE 9(2): e89470. https://doi.org/10.1371/journal.pone.0089470
Editor: Essa Yacoub, University of Minnesota, United States of America
Received: July 12, 2012; Accepted: January 22, 2014; Published: February 28, 2014
Copyright: © 2014 Eloyan et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Funding: The project was supported by grants P41EB015909 and R01EB012547 from the National Institute of Biomedical Imaging and Bioengineering, grant R01NS060, from the National Institute of Neurological Disorders and Stroke. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
This primer was created to give statisticians who are new to the field of imaging data analysis a quick overview of functional magnetic resonance imaging (fMRI) data and to describe tools for programming their own analyses of fMRI. It is targeted at those who would like to pursue programming using the R project [1]. The key benefit of programming on one's own is the ability to extend analyses and build and create new tools. Of course, the key benefits of using the R project are that it is open source, cross-platform, free (as in cost), and is an award-winning lingua franca of statisticians. Matlab (www.mathworks.com) is the standard scripting language in neuroimaging and signal processing, and has a widely used integrated development environment and GUI creation software. Both R and Matlab are cross-platform and have excellent graphics capabilities. Both are high level scripting languages that are easy to learn and have a large collection of add-ons and subroutines.
For those wanting more automated methods for the analysis of fMRI data, several programs exist to perform most standard analyses. In addition, power users should learn these programs for comparison purposes to avoid reinventing existing and well established tools. We give a nonexhaustive list of some of our favorite freely available software below.
SPM SPM (http://www.fil.ion.ucl.ac.uk/spm/) is a collection of open source Matlab scripts, often calling compiled code. SPM is arguably the most popular software for the analysis of fMRI data, in no small part due to a well developed GUI. It has methods for single-subject and multi-subject analyses and tools for displaying the results and data visualization. It can perform all of the basic preprocessing procedures of the imaging data. Moreover, it has an active user community with several user contributed add-ons. The SPM scripts are well documented, and are easily understood.
FSL FSL (http://www.fmrib.ox.ac.uk/fsl/) is a UNIX-based software that performs single- and multi-subject analyses, preprocessing and visualization of results. It has a graphical user interface, though is most effectively used at the command line, tying routines together using shell scripts. FSL is open source and can be compiled from scratch. In addition, pre-compiled binaries are available, notably for OSX. For Windows, FSL can be run in a virtual machine.
MIPAV MIPAV (http://mipav.cit.nih.gov/) is a JAVA software created by the US National Institutes of Health. It has very broad functionality for preprocessing, analysis and visualization of imaging data. However, it is most useful through an active collection of modules.
AFNI AFNI (http://afni.nimh.nih.gov/afni/) is a UNIX-based software. It offers a full analysis and processing suite. We often use it in conjunction with FSL and BASH scripting.
FreeSurfer FreeSurfer (http://surfer.nmr.mgh.harvard.edu/) is a UNIX-based software for registration (often atlas based) and analysis of medical imaging data. It is often discussed as having excellent registration and visualization tools.
ANTS ANTS [3] is a UNIX-based software for preprocessing of imaging data. It provides methods for linear and nonlinear registration of images.
In addition to the above tools for analysis and processing of fMRI data, there are R packages related to fMRI. We list a few packages here primarily focused on analysis postregistration.
AnalyzeFMRI [3] Analysis package for reading and writing of fMRI including a graphical user interface. It can also be used for performing temporal and spatial ICA.
oro.nifti [4] Package for reading in NIfTI files, a common fMRI imaging format. It can be used for reading compressed files.
oro.dicom [5] Package for reading in DICOM files, another common fMRI imaging format.
fmri [6] Package for reading in and analyzing fMRI data. Includes graphical user interface for fMRI analysis and adaptive smoothing for fMRI along with inference.
neuRosim [7] Package for simulating fMRI data. As a consequence, has tools for creating and investigating designs with HRF convolution models.
arf3DS4 [8] Package for activated region fitting for fMRI data in NIfTI format.
Rniftilib [9] Package for loading and writing the 3D or 4D images.
RNiftyReg [10] Package for 2D or 3D image registration.
brainwaver [11] Package for functional connectivity analysis.
tractor.base [12] Package for reading, writing, and graphical representation of images.
waveslim [13] Package for wavelet analysis of 1D, 2D and 3D images.
cudaBayesreg [14] Package for Compute Unified Device Architecture (CUDA) based Bayesian multilevel analysis of fMRI data.
We will use some of the packages above in the examples in this manuscript. In addition, a frequently updated more exhaustive list on CRAN of packages for medical image analysis can be found here:
The remainder of the manuscript is organized as follows. Sections 1, 2 and 3 describe the structure of the fMRI data, discuss ways of obtaining the data and give a brief overview of the preprocessing steps. Section 4 describes the different formats that can be used for viewing the fMRI data in R along with examples of obtaining a Gaussian smoother of the data matrix and the mask of the three dimensional image. A discussion of the hemodynamic response (HR) function and the role of the design matrix in task-based fMRI studies is presented in Section 5. The between-subject random effect models are discussed in Section 6. Section 7 gives an overview of the connectivity based analysis of fMRI data including methods based on seed-voxel analysis, singular value decomposition and independent component analysis. As the main goal of the manuscript is centered on the statistical analysis of fMRI data via R, we do not present visualization tools in detail, however, we briefly discuss a few of the visualization tools in R or other software in Section 8. Finally, the last section completes the article by pointing out useful further reading material.
FMRI is a modern brain imaging measurement technology. As its name suggests, the technology is used to explore brain function (activity) by obtaining several images of the brain over time using an MRI machine. Standard fMRI images the so-called blood oxygen level dependent (BOLD) signal, described further below. An extensive overview of the statistical methods developed for the analysis of fMRI data is presented by [15], [16], [17] and [18]. Herein, we will discuss the implementation of some of these commonly used techniques in R, while providing brief discussions as to the methods and the interpretation of the results. BOLD fMRI is not the only functional brain imaging technology available using an MRI scanner (see [19] for an overview). For example, arterial spin labeling [20] is a closely related functional MRI technology. However, the term fMRI, when used without qualifying statements, refers to BOLD fMRI. We note that positron emissions tomography (PET), single photon emissions computed tomography, electroencephalograms (EEG), magnetoencephalograms (MEG) and other non-magnetic resonance based measurement devices can be used for functional brain imaging, each with different goals, limitations, strengths and scientific interpretations.
We will not thoroughly discuss the details of the methods by which the MRI scanners use the principles of magnetic resonance to achieve different images. However, newcomers to the area may be surprised to find the amazing variety of biological signals that one can visualize using MRI. These include separate modalities to study white matter, gray matter, overall structure, white matter tracts, tracers injected into the body and lesions in the brain. There are multiple types of structural and functional images, all acquired using the principles of MR. See [19] and [21] for a general introduction to MR imaging.
The BOLD signal ([22], [23]) measures the cerebral hemodynamic changes concomitant to neuronal activation of the brain. A localized increase in neural activity results in an increase in blood flow (“reactive hyperemia” [24],[25]) to the activated region when excess of oxygenated hemoglobin is delivered to the region. Consequently, a reduction in deoxyhemoglobin concentration follows, resulting in an increase in magnetic field homogeneity. The gradient echo EPI sequence ([26], [27]) allows to study the activation of specific regions of the brain during a task. Deoxyhemoglobin serves as an endogenous susceptibility contrast agent allowing MRI to report on local hemodynamic changes. Hence, the BOLD signal in response to a stimulus is given a name - the hemodynamic response. Figure 1 shows an example hemodynamic response function (HRF). In Figure 1, an initial delay of the response can be observed since it takes time for the vasculature of the brain to respond to the need for oxygen after the stimulus. This is followed by a subsequent brief undershoot [28]. The origin of the undershoot is controversial. One explanation of the undershoot is that there is ongoing oxygen consumption after reaching the point of origin. The other is that the excess volume of deoxygenated blood results in delayed vascular compliance.
PowerPoint slide larger image original image Figure 1. Anatomy of the HRF.In BOLD fMRI the scanner records images approximately every second (the so-called TR or time of repetition). At each time point, a 3D image of the brain is obtained. Often, in BOLD fMRI, the image is collected one slice at a time for each time point instead of scanning the 3D image immediately. Hence, a slice time correction (e.g., by Fourier transformation [28]) is applied to the resulting images as part of the preprocessing. Empirical face validity of the methodology can be demonstrated by simple motor or visual tasks. For example, consider a “finger tapping” task where a subject is asked to tap their finger while in the scanner and then rest, repeating this sequence for a certain time interval. When the images are averaged over the times at rest and compared to the times when tapping the finger, and appropriate statistical tests are performed, motor areas of the brain are clearly activated (see [29] for overview of finger tapping task-related analysis in the literature). By virtue of the well-defined motor area of activation and ease of performing the task, finger tapping tasks are commonly used for scanner and methodological validation. They are also of intrinsic interest in the study of motor areas. Various other tasks can be designed depending on the goals of the study. More recently, resting state fMRI has been used to explore the brain function when the subject is at rest. Section 5.1 provides details on how the experimental design can inform the statistical analysis of the data.
The development of the BOLD fMRI technology facilitated research in the analysis of cognitive function of the brain [15]. For instance, one may be interested in identifying which regions are involved in performing a mental operation or a cognitive task, how active is a brain region during a task, what is the shape of the time course of activation or quantifying the extent of connectivity between brain regions during tasks. Since there are several sources of variability introduced during data collection, statistical methods are required in order to accurately obtain the regions of interest. For instance, in Section 5.1 we discuss two types of experimental designs used frequently in fMRI analysis: block design and event-related design. It can be shown that the block designs have high detection power (i.e., the ability of the method to detect activation), whereas the event related designs have better estimation efficiency in that they can be used to estimate the shape of the HRF [30]. In Section 5.1, we show several statistical methods that can be used to answer the questions posed above.
For the purpose of statistical analysis the observed fMRI dataset contains a 3D array of intensities for each time point. The number of timepoints can vary depending on the length of the scan and is usually in hundreds. Axial slices of the 3D images at each time point for a subject at rest are shown in Figure 2 top panel. A voxel (or volumetric pixel) is defined as a 3D unit element in the image. The background of the image is generally removed and the voxels in the brain are used for analysis. Depending on the research question one may choose to organize the data differently. For instance, the time courses of each voxel (as shown in Figure 2 bottom panel) may be analyzed (voxel-wise) to identify the voxels that are activated during a task. The 3D images at each time point may be vectorized and concatenated to obtain a matrix that can be further analyzed to discover brain networks. The correlations of time courses for each pair of voxels in the brain are often used in connectivity analysis, that is in identifying areas of the brain that activate and deactivate together.
PowerPoint slide larger image original image Figure 2. Resting state fMRI image for 1 subject.The top panel shows one axial slice from the 3D image acquared at each of the 152 time points. denotes a voxel in the brain with coordinates [46, 64, 37]. The plot below shows the intensity of the image at voxel (on the y-axis) at each time point (the x-axis corresponds to time).
As many neuroscience researchers are joining open science initiatives promoting public access for data, numerous websites offer datasets that can be downloaded freely. However, we have a few favorite canonical datasets that are good starting points for analysis. The first are SPM's example datasets that can be used for single-subject and inter-subject analyses. The second dataset of interest is the 1000 Functional Connectomes Project (FCP)/INDI [31] resting-state data. Most of the SPM datasets are preprocessed and ready for analysis; however, the FCP data need extensive preprocessing. UNIX shell scripts using FSL and AFNI can be found on their website for those who are interested in preprocessing.
The SPM example datasets can be downloaded at
The 1,000 FCP dataset can be downloaded at
As mentioned above, the 1,000 FCP data involves the most preprocessing. We therefore put a processed image in .nii format along with the corresponding 2D matrix in an R data file format with the extension .rda (the .rda files can be loaded by using the load() function in R) for download on our servers at
For illustrating the analysis in this manuscript, we choose images from the Kennedy Kreiger Institute (the home institution of some of the coauthors of this manuscript).
We also note that FSL has an evaluation suite, the FEEDS dataset, that we have found very useful, both for learning FSL as well as debugging new methods. The FSL example dataset can be downloaded at:
Preprocessing is an essential component in the analysis of fMRI data. Standard preprocessing requires, for example, skull stripping, evaluation and corrections for motion, coregistering, spatial smoothing and registration to a template. We discuss two preprocessing steps in particular: spatial smoothing and registration to a template. However, for newcomers to the area, we suggest leaving all the steps, except for possibly the spatial smoothing to standard software, or acquiring data where preprocessing has already been completed in an acceptable manner. For those readers interested in pursuing preprocessing as a research endeavor, good theoretical foundations are given by [32], [33] and [34].
Warping the images of a subject to a template, often called spatial registration or spatial normalization, is the key step of attempting to put subject images into a common space. That is, each registered image is transformed to the same space and hence the spatial locations of the voxels in the brain are ideally interpretable across subjects, in other words, voxel 1 for subject 1 is the same as voxel 1 for subject 2. In functional neuroimaging, this process relies on numerous assumptions. For example, it presumes anatomical similarity across subjects at the functional resolution of interest. This assumption becomes especially problematic when studying diseased brains, where anatomy can be drastically different than the healthy control template brain. Also, it presumes localization of the brain function of interest to common anatomical locales across subjects.
Setting these issues aside, registration is accomplished by mathematically warping each subject's collection of images over time to a so-called “template”. A template is typically a highly detailed structural image obtained by imaging a subject repeatedly or warping several subjects into one common image. The template is useful for overlaying results on top of it which then helps in contextualizing the results. It is important to know what template was used to warp the images one is analyzing. Because it is often a highly accurate structural image (or there is an associated image in the same space), it is useful for interpretation, as neuroscientists can visually relate findings to important structural landmarks. Moreover, some templates are in wide use and key regions have been tabulated (see [35] for more details). We normally use templates created by the Montreal Neurological Institute (MNI) at:
The International Consortium for Brain Mapping (http://www.loni.ucla.edu/ICBM/) releases labeled brain atlases. We should also mention the famous Talairach Atlas, which is a coordinate system created before MR imaging was in place. Relating results to Talairach coordinates is often done to put results in a historical perspective familiar to neuroscientists. Some templates have conversion utilities to convert their coordinates to Talairach coordinates.
Finally, spatial smoothing using a filter is a step that can be easily done in R. The choice of whether spatial smoothing is necessary to perform is generally a designed component of a statistician's analysis plan. Hence, of all the preprocessing steps, this is the one we suggest to be performed in R if deemed necessary.
Two key data representations are the array format and the matrix format. The first is the more natural format for the data. Suppose I is an array representing a 3D image of the brain, where I[i,j,k] is the intensity of voxel i,j,k. If the image is collected repeatedly over time (as in fMRI), then 4D arrays ( I[i,j,k,t]) are useful, though often are very large.
Medical images in neuroimaging are usually saved in one of a few available formats. However, all formats have common properties, for instance, they are binary and the data are stored as a long floating point vector. Header information includes descriptive features of the image such as the physical dimensions of a voxel and the size of each of the three dimensions, as well as other information.
DICOM Digital Imaging and Communications in Medicine (DICOM) is a standard image format from National Electrical Manufacturers Association often used as scanner output. DICOM images have the header information built into the binary file. 3D DICOM files are sometimes split into several files of 2D slices, but can be combined into a single file.
Analyze The Analyze format includes two files for each 3D image: a header file with the extension .hdr and an image file with the extension .img. The header and image files have to have the same filename and have to be saved in the same folder. The 4D fMRI images are often stored as separate files, where the 3D image for each time point is saved as a separate file. In Analyze 7.5 format, the 4D images can be saved using one .img and one .hdr file.
NIfTI NIfTI files typically have the extension .nii and are a current standard for the analysis of fMRI data. NIfTI files can have more than 3 dimensions. Often NIfTI files are compressed using gzip and so the file extension is .nii.gz. If your analysis program cannot gunzip the file, simply uncompress it using any of the standard compression software tools and the result will be a standard (uncompressed) NIfTI file. The R package oro.nifti can be used to read in the gzipped files directly.
In this example, we read in a sample 3D so-called contrast image and smooth it using a 3D Gaussian filter [36] included in the package AnalyzeFMRI. For illustration, we are using one of the SPM datasets (see [37], for more information about the dataset) described in Section 2. First we read in the data.
Here we assume that imageFile.img is an Analyze image file with a header and image pair: imageFile.hdr and imageFile.img. Notice this data is 3D, yet the program reads it in as 4D, where the first three dimensions correspond to the size of the 3D image and the fourth dimension is set to 1. Hence the array subscripting drops the dimension by 1. It is often useful to read in the header information.
We can visualize this image using the image function.
# visualize the image
for(k in 1: dim(img)[3])
image(tempImg[,,k], axes = FALSE)
Note that the creation of the temporary image tempImg in this example is only useful as the background of the raw image is stored as NA. Note also by default the image command does not scale the image relative to the actual voxel dimensions. The parameter asp can be used to specify the aspect ratio y/x which would result in the correct voxel dimensions. The following code smooths the image using the function GaussSmoothArray.
for(k in 1: dim(img)[3])
image(simg[,,k], axes = FALSE)
Here the mask statement contains the smoother to the non-background areas, i.e., the smoother is applied to the non-background voxels which have the value of 1 in the mask. The voxdim variable contains the voxel dimensions, in this case three millimeters cubed. The voxel dimensions can be found using visualization programs or inspecting the hdr variable. The ksize and sigma parameters are smoothing parameters. The sigma parameter is the standard deviation of the Gaussian kernel (specified in mm) while ksize is the truncation width. See the documentation for GaussSmoothArray for more information. Figure 3 shows an axial slice of the original image along with the corresponding smoothed image (see Section 8 for more details on displaying images).
PowerPoint slide larger image original imageFigure 3. A slice of the fMRI image plotted by using the image() function in R (left) along with the corresponding smoothed slice (right).
Here the red color corresponds to high intensity values, followed by yellow and white as the values decrease.
The array format, while useful for smoothing, is not a convenient data structure for most other analyses. Most importantly, it is not a parsimonious structure, since it includes all of the voxels outside of the brain, i.e., the background of the image. In some fMRI analysis, where the spatial structure of the voxels is disregarded, a common mask is applied to the 3D data and the data are then saved in a long 1D vector. In fMRI, if there is one vector per time or subject, these can be stacked into a matrix for analysis. Under this structure, any analysis of the resulting matrix must be invariant to spatial location. The spatial location information can be kept elsewhere for back reconstruction of the image. For example, if there are 30,000 nonbackground voxels, one could retain the 3 30,000 matrix of location indices for those voxels. Masks can be created in a variety of ways. They can be obtained from the template or from the data itself. The mask is an important structure to retain, as it can map a vector back to an array format for display.
Here, we will take a collection of 3D images for one subject (one 3D image per time point) and create a mask. Using the mask, we will then create the data matrix. First, we read in the list of the files (with absolute paths).
# the directory where the images are stored
# the collection of images
It is useful to have the image dimensions assigned to an object. We are assuming that all images are registered and have the same dimensions.
# obtain the image dimensions by loading the first image
Now we loop over the time courses and collect them into a matrix by grabbing the relevant indices across the time courses. We then find the corresponding non-background indices and place them in a vector called mask.
for (file in files)
# the mask is a list of indices
Next, we loop over scans and collect them into a matrix by grabbing the relevant indices across scans.
# now cycle through the data and get the data into a matrix format
for (file in files)
# load file and make 3D
When designing the experiment in task-based fMRI the aim is to maximize the statistical power in detecting activation during the task while assuring the validity of the psychological results. In task-based fMRI, there are two common designs: event-related and block design. Event-related refers to multiple stimuli that are assumed to occur instantaneously and have randomized time between events. A simple task would be to have the participant push a button approximately every 20 seconds to look at motor cortex activation though the events do not have to be equidistant. If a “continuous” task is performed, such as sequential finger tapping, for a certain time interval followed by a short block of rest, the design is referred to as block design. Examples of stimulus functions in each design are presented in Figure 4. The choice of event-related design versus a block design depends on the final goal of the experiment. The block designs lead to higher detection power, however, the subjects may learn the patterns of the experiment which may affect the interpretation of the results. Event-related designs reduce the effects of learning, boredom and other events unrelated to the task while exhibiting loss in detection power (see [17] for more discussion on designs and [38] on optimization of design parameters).
PowerPoint slide larger image original image Figure 4. Stimulus vectors for an event-related design (top) and a block design (bottom).The spikes correspond to the onsets of stimuli. There is a sustained period of activity/task in the block design. In both cases, the spacings between events can be unequal.
Suppose we have fMRI data for one subject during a task. The goal is to identify the area of the brain that activates during the task, which is performed by the subject. Typically, the fMRI data can be modeled as a sum of response, drift and noise:
The response (i.e., the BOLD response) is often modeled using a linear time invariant system (see [17] for other ways of modeling the BOLD response). Assume the response is a linear combination of responses from different stimuli. That is, the response for voxel , at time , , can be written as , where is the response from stimulus at time in voxel , is the number of stimuli, is the number of scans and is the number of non-background voxels in each scan. The response from each stimulus depends on the stimulus function and the HRF, often modeled as a convolution of the two functions , where is the HRF at voxel and is the stimulus function corresponding to stimulus . The modeling of the HRF function is an important part of the analysis and is discussed in more detail in Section 5.2.
In addition to the random variability, periodic noise may be present in fMRI data unrelated to the true BOLD response. This can result, for example, from the patient's heartbeat or other systematic effects. As a result, some voxels in the brain may show considerable drifts over time. The drift can be linear or nonlinear, hence a flexible polynomial model is often employed to allow for nonlinear effects in the drift. For example, a polynomial drift of order for scan at voxel can be written as , where . Finally, a random error is added to the observed data, , at voxel scan yielding
In Section 4.2, we discussed how to construct a matrix from the 4D fMRI array. The resulting matrix is dimensional: the columns of which correspond to time courses in each voxel and the rows are the spacial images for each scan. We define the vectorized 3D image for each time point t as and and
The general linear model (GLM) can be written in matrix notation as
Note that, in the above equation, is derived based on specific hemodynamic assumptions, and is interpreted as the size of the hemodynamic response (see [39] for more details). The random errors are generally not independent in time. The correlation structure can be specified based on assumptions one is willing to make about the dependence structure over time. Misspecification of the correlation structure in GLM can lead to biased estimation of coefficients. For example, the temporal correlations can be modeled using a first order autoregressive model [40], that is, , where and are independent and identically distributed. After estimating the parameters in the GLM model, the p-value maps of the coefficients are thresholded to find the voxels where the corresponding is significantly different from 0.
In Section 1 we discussed the HRF function, an example of which is shown in Figure 1. The estimation of the HRF function is often of interest in fMRI analysis. The HRF model in SPM implies that the HRF is a discretized difference of two gamma functions (see [15], [40])
for . Typically, it is also normalized, either divided by its integral or maximum, to have an average or peak value of 1 respectively. The logic behind this specification of the HRF is that the initial post-stimulus increase in blood oxygenation is represented by , with the subsequent depletion represented by .
Notice that is a gamma density with shape and rate , and is a gamma density with shape and rate . The SPM documentation defines as the TR, as the delay of the response to the onset, as the delay of the undershoot relative to the onset, as the dispersion of the response, as the dispersion of the undershoot and as the ratio of the response to the undershoot. Conceptually, these are illustrated in Figure 1. We briefly give motivation for these parameters. Note, the mean of density is while the mean of density is . Hence, if and are given in seconds and is seconds per scan, can be thought of as the delay to the mean (not modal) response in TR units while is the delay to the mean undershoot in TR units. Given that the gamma variance is , and the mean is , is the factor by which the response gamma distribution mean is scaled to obtain its variance. Similarly, is the factor by which the undershoot gamma distribution mean is scaled to obtain its variance. Since the integral of and are 1, is the ratio of the two functions comprising the HRF. Note, it is not, as its name suggests, the ratio of the peak value of to the peak value of .
The following example illustrates the use of the fmri.stimulus function in R to create the expected HRF as described by [40] and [6]. More appropriate choices of an HRF, such as spatially varying HRF and dynamic sets of basis functions may also be used [41]. For instance, finite impulse response models and Fourier basis sets are more flexible alternatives to the canonical HRF.
The R package fmri can be used to analyze fMRI data from one subject and identify activation during a task. In the following example, suppose imageFile.nii is a NIfTI data file. We first extract the 4D numeric array using the function read.NIFTI. The onset times of the stimulus (5.25, 21.45, 100.12, 223.5) and the durations of ON stimuli in scans (5, 5, 10, 10) are known from the experimental design. Next, we create the expected response for each stimulus (defined by above, assuming that the HRF is the same in each voxel) using the function fmri.stimulus. The function fmri.design is invoked to create the design matrix. Finally, the coefficients of the hemodynamic response to the stimulus are evaluated in the GLM [42].