Permalink
Browse files

Basic code implementing simple methods for analyzing fMRI data.

  • Loading branch information...
johnmyleswhite committed Jan 13, 2010
0 parents commit c17824545286de30cac55cd094ab56cb6af6275d
Showing with 243 additions and 0 deletions.
  1. +1 −0 .gitignore
  2. +29 −0 README.markdown
  3. +5 −0 example_data.csv
  4. +60 −0 example_output.1D
  5. +63 −0 fmri_utilities.R
  6. +15 −0 practical_example.R
  7. +34 −0 tests.R
  8. +36 −0 toy_example.R
@@ -0,0 +1 @@
+.DS_Store
@@ -0,0 +1,29 @@
+First, we define two functions for converting between real time in seconds and
+the index of an image defined by a TR expressed in seconds:
+<pre>
+real.time.to.image.number
+image.number.to.real.time
+</pre>
+
+Next, we define a canonical HRF uses AFNI's version of Mark Cohen's gamma
+variate model.
+<pre>
+gamma.hrf
+hrf
+</pre>
+
+Then, we build up a representation of our experiment's events in the form
+that we will convolve with our canonical HRF:
+<pre>
+build.tr.representation
+</pre>
+
+Finally, we convolve this with the HRF to produce a finished parametric
+regressor.
+<pre>
+build.parametric.regressor
+</pre>
+
+For normal use you probably just want to import your experimental event data
+from a CSV file into a data frame and then use `build.parametric.regressor()`
+to generate the regressors for your design matrix.
@@ -0,0 +1,5 @@
+"Stimulus Onset","Stimulus Intensity"
+0.0,1
+17.3,10
+54.9,5
+104.9,20
@@ -0,0 +1,60 @@
+9.50218345135895e-14
+8.96393728316326
+89.8344185918076
+75.8426638806616
+23.2526632910349
+4.09246340038437
+0.50702515854648
+0.0493012453536172
+9.60823460595001e-14
+-6.27822835179074e-14
+89.6393728316318
+898.344185918075
+758.426638806616
+232.526632910350
+40.9246340038445
+5.07025158546439
+0.493012453537484
+-1.04142233808421e-13
+-1.84741111297626e-13
+7.72052405422915e-14
+1.48471616427484e-14
+-5.64192142424438e-14
+-6.15096696628146e-14
+-8.73861513830332e-14
+-7.16905805035564e-14
+-3.13911417589537e-14
+4.6238303401702e-14
+2.96094823618239e-13
+44.8196864158161
+449.172092959038
+379.213319403308
+116.263316455174
+20.4623170019221
+2.53512579273216
+0.246506226768664
+-2.71490955753113e-14
+-7.4660012832106e-14
+-2.71490955753113e-14
+-1.62894573451868e-13
+-4.07236433629669e-14
+3.13911417589537e-14
+4.02994387446027e-14
+-8.05988774892054e-14
+-4.02994387446027e-14
+-1.94285715210821e-13
+-1.31927636311278e-13
+1.20898316233808e-13
+-6.5751715846457e-14
+-2.29070493916689e-14
+2.05102932979110e-13
+-1.06263256900242e-13
+1.98103556776100e-13
+3.902682488951e-14
+179.278745663264
+1796.68837183615
+1516.85327761323
+465.053265820697
+81.8492680076887
+10.1405031709284
+0.986024907074948
@@ -0,0 +1,63 @@
+# We assume that image 1 occurs at t = 0, image 2 occurs at t = tr, and so on.
+# Therefore
+# (1) real.time == (image - 1) * tr
+# (2) round(real.time) / tr + 1 == image
+real.time.to.image.number <- function(real.time, tr = 2)
+{
+ return(round(real.time / tr) + 1)
+}
+
+image.number.to.real.time <- function(image.number, tr = 2)
+{
+ return((image.number - 1) * tr)
+}
+
+# We use Mark Cohen's gamma variate model of the HRF.
+gamma.hrf <- function(t, p = 8.6, q = 0.547)
+{
+ return((t / (p * q))^p * exp(p - t / q))
+}
+
+hrf <- function(tr = 2, model = 'gamma')
+{
+ return(100 * gamma.hrf(seq(0, 14, by = tr)))
+}
+
+# 'events' must be a data frame whose rows contain events with real time
+# onsets and parametric values that will be convolved with the HRF to
+# generate a regressor.
+build.tr.representation <- function(events, tr = 2, images = 7,
+ onsets.key = 'onsets', values.key = 'values')
+{
+ tr.representation <- rep(0, images)
+
+ for (index in 1:nrow(events))
+ {
+ real.time <- events[index, onsets.key]
+ image.number <- real.time.to.image.number(real.time, tr)
+ parametric.value <- events[index, values.key]
+ tr.representation[image.number] <- parametric.value
+ }
+
+ return(tr.representation)
+}
+
+# 'events' must be a data frame whose rows contain events with real time
+# onsets and parametric values that will be convolved with the HRF to
+# generate a regressor.
+
+# If you do not have a data frame in memory, pass an empty data frame in lieu
+# of 'events' and provide an explicit filename as an argument.
+build.parametric.regressor <- function(events, tr = 2, images = 7,
+ onsets.key = 'onsets', values.key = 'values',
+ filename = NULL, separator = ',')
+{
+ if (! is.null(filename))
+ {
+ events <- read.csv(filename, header = TRUE, sep = separator)
+ }
+
+ tr.representation <- build.tr.representation(events, tr, images, onsets.key, values.key)
+
+ return(convolve(tr.representation, rev(hrf(tr)), type = 'o')[1:images])
+}
@@ -0,0 +1,15 @@
+source('fmri_utilities.R')
+
+regressor <- build.parametric.regressor(NULL,
+ tr = 2,
+ images = 60,
+ onsets.key = 'Stimulus.Onset',
+ values.key = 'Stimulus.Intensity',
+ filename = 'example_data.csv',
+ separator = ',')
+
+write.table(regressor,
+ file = 'example_output.1D',
+ sep = ' ',
+ row.names = FALSE,
+ col.names = FALSE)
34 tests.R
@@ -0,0 +1,34 @@
+source('fmri_utilities.R')
+
+assert <- function(obtained.value, desired.value)
+{
+ if (obtained.value != desired.value)
+ {
+ cat(paste('Got: [', obtained.value, ']\nExpected: [', desired.value, ']\n', sep = ''))
+ }
+}
+
+assert.within.delta <- function(obtained.value, desired.value, delta = 0.00001)
+{
+ if (abs(obtained.value - desired.value) > delta)
+ {
+ cat(paste('Got: [', obtained.value, ']\nExpected: [', desired.value, ']\n', sep = ''))
+ }
+}
+
+assert(real.time.to.image.number(0.0, tr = 2), 1)
+assert(real.time.to.image.number(0.9, tr = 2), 1)
+assert(real.time.to.image.number(1.9, tr = 2), 2)
+assert(real.time.to.image.number(2.9, tr = 2), 2)
+assert(real.time.to.image.number(2.9999999999, tr = 2), 2)
+assert(real.time.to.image.number(3.0, tr = 2), 3)
+assert(real.time.to.image.number(3.9, tr = 2), 3)
+
+events <- data.frame(onsets = c(0.0, 10.0),
+ values = c(1, 5))
+
+expected.tr.representation <- c(1, 0, 0, 0, 0, 5, 0)
+assert.within.delta(build.tr.representation(events), expected.tr.representation)
+
+expected.parametric.regressor <- c(-8.120488e-15, 8.963937e+00, 8.983442e+01, 7.584266e+01, 2.325266e+01, 4.092463e+00, 4.532671e+01)
+assert.within.delta(build.parametric.regressor(events), expected.parametric.regressor)
@@ -0,0 +1,36 @@
+source('fmri_utilities.R')
+
+events <- data.frame(onsets = c(0.0, 10.0, 25.0),
+ values = c(2, 3, 4))
+
+power.events <- events
+power.events$values <- power.events^1.01
+
+log.events <- events
+log.events$values <- log(log.events$values)
+
+sqrt.events <- events
+sqrt.events$values <- sqrt(sqrt.events$values)
+
+plot(1:100,
+ build.parametric.regressor(power.events, images = 100),
+ main = 'Expected Hemodynamic Response',
+ xlab = 'Time (s)',
+ ylab = 'Arbitrary Units',
+ type = 'l',
+ col = 'purple')
+
+points(1:100,
+ build.parametric.regressor(events, images = 100),
+ type = 'l',
+ col = 'red')
+
+points(1:100,
+ build.parametric.regressor(log.events, images = 100),
+ type = 'l',
+ col = 'blue')
+
+points(1:100,
+ build.parametric.regressor(sqrt.events, images = 100),
+ type = 'l',
+ col = 'green')

0 comments on commit c178245

Please sign in to comment.