Overview of the fmrireg package
The analysis of Functional Magnetic Resonance Imaging (fMRI) has a lot of complexities due to the large amount of data, the fact that the data is in the form of multi-dimensional images (usually, 4-dimensional – x,y,z, and time), the fact that the hemodynamic response is sluggish and delayed (4-6s) with respect to stimulus presentation or task events, the existence of physiological noise and other “nuisance” signals, issues related to multiple comparisons, and the existence of spatially and temporally auto-correlated errors in the signal.
However, if we push all that complexity to the side for a moment,
standard “univariate” fMRI analysis involves simple multiple regression
analysis applied to every 3d location (“voxel”) in the brain. The goal
of the fmrireg package is to provide a convenient
interface, in the spirit of other R regression modeling tools, for
analyzing fMRI data. We therefore use a R’s formula interface to specify
regression models in a way that makes for convenient and readable
code.
The building blocks of an fMRI regression model
The basic, or lowest level, unit of an fMRI regression model is the hemodynamic reesponse function (HRF). The HRF is a model of the mapping, or transfer function, between an external stimulus (or neural event) and a BOLD signal. A commonly used HRF is the so-called SPM canonical HRF:
Currently, fmrireg provides several builtin HRFs that
can be used to model BOLD activity. See the vignette “Hemodynamic
Response Functions” for a detailed overview. Choosing the correct HRF
depends on the use case and aims of the analysis, but the SPM canaonical
function is usually a good default choice.
Regressors
A “regressor” is simply an independent variable in a multiple regression model. In fMRI analysis, however, regressors are constructed by “convolving” a set of onset times–the vector of times at which certain task events occur–with a chosen HRF. Imagine we have a stimulus that is presented every 16 seconds with 8 repetitions. We can construct a “regressor” by convolving the chosen HRF with a vector onset times (in seconds), whereby each element of the vector corresponds to time at which a stimulus occurs.
Below we use the function gen_hrf to create a special
function that has the class “HRF”. An alternative is to use the pre-made
HRF objects, for example, HRF_SPMG1 if we want to use th
SPM canonical function. Once we create a regressor with the
regressor constructor function, we can evaluate it at any
arbitrary time points. In other words, just as the HRFs are functions of
time, so too are regressors. The difference is that
regressor, when evaluated at time t sum over all
the constituent HRFs, the domain of which are determined by the onset
vector.
In the code below, we evaluate the regressor at every time point between 1 and 140.
onsets <- seq(0, 16*8, by=16)
hrf <- gen_hrf(hrf_spmg1)
reg <- regressor(onsets, hrf=hrf)
grid <- seq(0,140, by=1)
plot(grid, evaluate(reg, grid), type='l', ylab="BOLD activity", xlab="time", main="ISI = 16")
To better illustrate how a regressor sums over its constituent HRFs, we can place the onsets closer together, so that the HRFs overlap in time:
onsets <- seq(0, 16*8, by=6)
hrf <- gen_hrf(hrf_spmg1)
reg <- regressor(onsets, hrf=hrf)
grid <- seq(0,140, by=1)
plot(grid, evaluate(reg, grid), type='l', ylab="BOLD activity", xlab="time", main="ISI = 6")
In fmrireg it is generally not necessary to manually
create regressors as we have done above. Rather, regressors are created
“behind the scenes” by the higher-level modeling functions, such as
event_model. However, the facility to create individual
regressors, for the purposes of learning fMRI analysis or in conducting
special analyses, may still be useful.
Event Models
An HRF is a function from time to BOLD activity: f(t) -> y. A regressor is a function of time that sums up BOLD activity over a sequence of spaced events. In conventional fMRI analysis, a regressor refers to a single column of a design matrix. For example, if we have two task conditions, face and scene, we would generate one face regressor and one scene regressor. However, in statistical modeling of experiments, we generally think not in terms of single columns of a design matrix, but rather in terms of the “factors”, which in this case might call “stim_type”, with levels “scene” and “face”. A complex regression model often involves a factorial combination of different categorical variables, which can be then be examined for main effects, interactions, and other contrasts.
Before we discuss event models, we need to briefly introduce the
fmrireg’s class for describing how scans are sampled over
time and organized in blocks (sometimes also referred
to as runs or scans). A “block” simply refers to a consecutive set of
images (usually a few hundred), collected back to back over the course
of one fMRI scan. An fMRI “session” consists of a number of these
scanning blocks. Each individual image in a block is separated by a
fixed interval, called the “TR” (repetition time). We encode this
information (the number of blocks, the number of images in each block,
and the TR) in a sampling_frame object. Below, we create a
sampling_frame for an experiment with 5 blocks, each
containing 200 images, and with a TR of 2s.
nblocks <- 5
sframe <- sampling_frame(rep(200, 5), TR=2)
sframe## sampling_frame:
## number of blocks: 5
## blocklens: 200, 200, 200, 200, 200
## TR: 2s 2s 2s 2s 2s
## start_time: 1s 1s 1s 1s 1s