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## $blocklens
## [1] 200 200 200 200 200
##
## $TR
## [1] 2 2 2 2 2
##
## $start_time
## [1] 1 1 1 1 1
##
## $blockids
## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [38] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [75] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [112] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [149] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [186] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [223] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [260] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [297] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [334] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [371] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3
## [408] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
## [445] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
## [482] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
## [519] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
## [556] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
## [593] 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
## [630] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
## [667] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
## [704] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
## [741] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
## [778] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5
## [815] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
## [852] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
## [889] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
## [926] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
## [963] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
## [1000] 5
##
## $time
## [1] 1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35
## [19] 37 39 41 43 45 47 49 51 53 55 57 59 61 63 65 67 69 71
## [37] 73 75 77 79 81 83 85 87 89 91 93 95 97 99 101 103 105 107
## [55] 109 111 113 115 117 119 121 123 125 127 129 131 133 135 137 139 141 143
## [73] 145 147 149 151 153 155 157 159 161 163 165 167 169 171 173 175 177 179
## [91] 181 183 185 187 189 191 193 195 197 199 201 203 205 207 209 211 213 215
## [109] 217 219 221 223 225 227 229 231 233 235 237 239 241 243 245 247 249 251
## [127] 253 255 257 259 261 263 265 267 269 271 273 275 277 279 281 283 285 287
## [145] 289 291 293 295 297 299 301 303 305 307 309 311 313 315 317 319 321 323
## [163] 325 327 329 331 333 335 337 339 341 343 345 347 349 351 353 355 357 359
## [181] 361 363 365 367 369 371 373 375 377 379 381 383 385 387 389 391 393 395
## [199] 397 399 1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31
## [217] 33 35 37 39 41 43 45 47 49 51 53 55 57 59 61 63 65 67
## [235] 69 71 73 75 77 79 81 83 85 87 89 91 93 95 97 99 101 103
## [253] 105 107 109 111 113 115 117 119 121 123 125 127 129 131 133 135 137 139
## [271] 141 143 145 147 149 151 153 155 157 159 161 163 165 167 169 171 173 175
## [289] 177 179 181 183 185 187 189 191 193 195 197 199 201 203 205 207 209 211
## [307] 213 215 217 219 221 223 225 227 229 231 233 235 237 239 241 243 245 247
## [325] 249 251 253 255 257 259 261 263 265 267 269 271 273 275 277 279 281 283
## [343] 285 287 289 291 293 295 297 299 301 303 305 307 309 311 313 315 317 319
## [361] 321 323 325 327 329 331 333 335 337 339 341 343 345 347 349 351 353 355
## [379] 357 359 361 363 365 367 369 371 373 375 377 379 381 383 385 387 389 391
## [397] 393 395 397 399 1 3 5 7 9 11 13 15 17 19 21 23 25 27
## [415] 29 31 33 35 37 39 41 43 45 47 49 51 53 55 57 59 61 63
## [433] 65 67 69 71 73 75 77 79 81 83 85 87 89 91 93 95 97 99
## [451] 101 103 105 107 109 111 113 115 117 119 121 123 125 127 129 131 133 135
## [469] 137 139 141 143 145 147 149 151 153 155 157 159 161 163 165 167 169 171
## [487] 173 175 177 179 181 183 185 187 189 191 193 195 197 199 201 203 205 207
## [505] 209 211 213 215 217 219 221 223 225 227 229 231 233 235 237 239 241 243
## [523] 245 247 249 251 253 255 257 259 261 263 265 267 269 271 273 275 277 279
## [541] 281 283 285 287 289 291 293 295 297 299 301 303 305 307 309 311 313 315
## [559] 317 319 321 323 325 327 329 331 333 335 337 339 341 343 345 347 349 351
## [577] 353 355 357 359 361 363 365 367 369 371 373 375 377 379 381 383 385 387
## [595] 389 391 393 395 397 399 1 3 5 7 9 11 13 15 17 19 21 23
## [613] 25 27 29 31 33 35 37 39 41 43 45 47 49 51 53 55 57 59
## [631] 61 63 65 67 69 71 73 75 77 79 81 83 85 87 89 91 93 95
## [649] 97 99 101 103 105 107 109 111 113 115 117 119 121 123 125 127 129 131
## [667] 133 135 137 139 141 143 145 147 149 151 153 155 157 159 161 163 165 167
## [685] 169 171 173 175 177 179 181 183 185 187 189 191 193 195 197 199 201 203
## [703] 205 207 209 211 213 215 217 219 221 223 225 227 229 231 233 235 237 239
## [721] 241 243 245 247 249 251 253 255 257 259 261 263 265 267 269 271 273 275
## [739] 277 279 281 283 285 287 289 291 293 295 297 299 301 303 305 307 309 311
## [757] 313 315 317 319 321 323 325 327 329 331 333 335 337 339 341 343 345 347
## [775] 349 351 353 355 357 359 361 363 365 367 369 371 373 375 377 379 381 383
## [793] 385 387 389 391 393 395 397 399 1 3 5 7 9 11 13 15 17 19
## [811] 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49 51 53 55
## [829] 57 59 61 63 65 67 69 71 73 75 77 79 81 83 85 87 89 91
## [847] 93 95 97 99 101 103 105 107 109 111 113 115 117 119 121 123 125 127
## [865] 129 131 133 135 137 139 141 143 145 147 149 151 153 155 157 159 161 163
## [883] 165 167 169 171 173 175 177 179 181 183 185 187 189 191 193 195 197 199
## [901] 201 203 205 207 209 211 213 215 217 219 221 223 225 227 229 231 233 235
## [919] 237 239 241 243 245 247 249 251 253 255 257 259 261 263 265 267 269 271
## [937] 273 275 277 279 281 283 285 287 289 291 293 295 297 299 301 303 305 307
## [955] 309 311 313 315 317 319 321 323 325 327 329 331 333 335 337 339 341 343
## [973] 345 347 349 351 353 355 357 359 361 363 365 367 369 371 373 375 377 379
## [991] 381 383 385 387 389 391 393 395 397 399
##
## $precision
## [1] 0.1
##
## attr(,"class")
## [1] "sampling_frame" "list"