# Presentation: 529735

## Outline: 6+ Hours of Edification

-  Philosophy (e.g., theory without equations)
-  Sample FMRI data
-  Theory underlying FMRI analyses: the **HRF**
-  “Simple” or “Fixed Shape” regression analysis
-  Theory and Hands-on examples
-  “Deconvolution” or “Variable Shape” analysis
-  Theory and Hands-on examples
-  Advanced Topics (followed by brain meltdown)

**Goals****: Conceptual **_**Understanding**_** ****+**** Prepare to Try It Yourself**

##  Signal = Measurable response to stimulus
 Noise = Components of measurement that interfere with detection of signal
 Statistical detection theory:
 Understand relationship between stimulus & signal
 Characterize noise statistically
 Can then devise methods to distinguish noise-only measurements from signal+noise measurements, and assess the methods’ reliability
 Methods and usefulness depend strongly on the assumptions
 Some methods are more “robust” against erroneous assumptions than others, but may be less sensitive

-  _**Noise**_ = Components of measurement that interfere with detection of signal
-  Statistical detection theory:
-  _Understand_ relationship between stimulus & signal
-  Characterize noise statistically
-  Can then devise methods to distinguish noise-only measurements from signal+noise measurements, and assess the methods’ reliability
-  Methods and usefulness depend strongly on the assumptions
-  Some methods are more “robust” against erroneous assumptions than others, but may be less sensitive

- Data Analysis Philosophy

## FMRI Philosopy: Signals and Noise

-  FMRI _Stimulus__****__Signal_ connection and _noise statistics_ are both complex and poorly characterized
-  Result: there is no “**best**” way to analyze FMRI time series data: there are only “**reasonable**” analysis methods
-  To deal with data, must make some assumptions about the signal and noise
-  Assumptions will be wrong, but must do ***something***
-  Different kinds of experiments require different kinds of analyses
-  Since signal models and questions you ask about the signal will vary
-  It is important to _**understand**_ what is going on, so you can select and evaluate “reasonable” analyses

## Meta-method for creating analysis methods

-  Write down a mathematical model connecting stimulus (or “activation”) to signal
-  Write down a statistical model for the noise
-  Combine them to produce an equation for measurements given signal+noise
-  Equation will have unknown parameters, which are to be estimated from the data
-  N.B.: signal may have zero strength (no “activation”)
-  Use statistical detection theory to produce an algorithm for processing the measurements to assess signal presence and characteristics
-  e.g., least squares fit of model parameters to data

## Time Series Analysis on Voxel Data

-  Most common forms of FMRI analysis involve fitting an activation+BOLD model to each voxel’s time series ***separately****** ***(AKA “univariate” analysis)
-  Some pre-processing steps do include inter-voxel computations; e.g.,
-  spatial smoothing to reduce noise
-  spatial registration to correct for subject motion
-  Result of model fits is a set of parameters at each voxel, estimated from that voxel’s data
-  e.g., activation amplitude (********* ***), delay, shape
-  “SPM” = statistical parametric map; e.g., ****** or *t* or *F*
-  Further analysis steps operate on individual SPMs
-  e.g., combining/contrasting data among subjects
-  sometimes called “second level analysis”

## Some Features of FMRI Voxel Time Series

-  FMRI only measures _changes_ due to neural “activity”
-  Baseline level of signal in a voxel means little or nothing about neural activity
-  Also, baseline level tends to drift around slowly (100 s time scale or so; mostly from small subject motions)
-  Therefore, an FMRI experiment must have at least 2 different neural conditions (“tasks” and/or “stimuli”)
-  Then statistically test for differences in the MRI signal level between conditions
-  Many experiments: one condition is “rest”
-  Baseline is modeled separately from activation signals, and _baseline model includes “rest” periods_
-  In AFNI, that is; in SPM, “rest” is modeled explicitly

## Why FMRI Analysis Is Hard

- Don’t know true relation between neural “activity” and BOLD signal:
- What *is* neural “activity”, anyway?
- What is connection between “activity” and hemodynamics and MRI signal?
- Noise in data is poorly characterized
- In space and in time, and in its origin
- Noise amplitude  BOLD signal
- Can some of this noise be removed by software?
- Makes both signal detection and statistical assessment hard
- Especially with 20,000+ voxels in the brain **=** 20,000+ activation decisions

## Why So Many Methods of Analysis?

- Different assumptions about activity-to-MRI signal connection
- Different assumptions about noise ( signal fluctuations of no interest) properties and statistics
- Different experiments and different questions about the results
**Result**: **** Many “reasonable” FMRI analysis methods
- Researchers _**must**_ understand the tools (models and software) in order to make choices and to detect glitches in the analysis!!

##  First sample: Block-trial FMRI data
 “Activation” occurs over a sustained period of time (say, 10 s or longer), usually from more than one stimulation event, in rapid succession
 BOLD (hemodynamic) response accumulates from multiple close-in-time neural activations and is large
 BOLD response is often visible in time series
 Noise magnitude about same as BOLD response
 Next 2 slides: same brain voxel in 3 (of 9) EPI runs
 black curve (noisy) = data
 red curve (above data) = ideal model response
 blue curve (within data) = model fitted to data
 somatosensory task (finger being rubbed)

-  “Activation” occurs over a sustained period of time (say, 10 s or longer), usually from more than one stimulation event, in rapid succession
-  BOLD (hemodynamic) response accumulates from multiple close-in-time neural activations and is large
-  BOLD response is often visible in time series
-  Noise magnitude about same as BOLD response
-  Next 2 slides: same brain voxel in 3 (of 9) EPI runs
-  **black curve** (noisy) = data
-  **red curve** (above data) = ideal model response
-  **blue curve** (within data) = model fitted to data
-  somatosensory task (finger being rubbed)

- Some Sample FMRI Data Time Series

## Same Voxel: Runs 1 and 2

- Block-trials: 27 s “on” / 27 s “off”; TR=2.5 s; 130 time points/run

**model fitted to data**

**data**

**model regressor**

**Noise**** ******** ****same size as ********signal**

## Same Voxel: Run 3 and Average of all 9

****** **Activation amplitude & shape vary among blocks!  Why???

## More Sample FMRI Data Time Series

-  _Second sample_: Event-Related FMRI
-  “Activation” occurs in single relatively brief intervals
-  “Events” can be randomly or regularly spaced in time
-  If events are randomly spaced in time, signal model itself _looks_ noise-like (to the pitiful human eye)
-  BOLD response to stimulus tends to be weaker, since fewer nearby-in-time “activations”
-   have overlapping signal changes
-   (hemodynamic responses)
-  Next slide: Visual stimulation experiment

- “Active” voxel shown in next slide

## Two Voxel Time Series from Same Run

**correlation with ideal = 0.56**

**correlation with ideal = –**** ****0.01**

- Lesson: ER-FMRI activation is not obvious via casual inspection

## More Event-Related Data

- _More Event-Related Data_

** ****White curve = Data (first 136 TRs)**
** ****Orange curve = Model fit ****(R****2**** ****=**** ****50%)**
** ****Green = Stimulus timing**

**Four different visual stimuli**

- Very good fit  for ER data (R2=10-20% more usual).
- Noise is as big as BOLD!

## Two Fundamental Principles Underlying Most FMRI Analyses (esp. GLM): HRF  Blobs

-  _**H**_emodynamic _**R**_esponse _**F**_unction
*Convolution* model for ***temporal*** relation between stimulus/activity and response
-  Activation _**Blobs**_
- Contiguous ***spatial*** regions whose voxel time series fit HRF model
*e.g.*, Reject isolated voxels even if HRF model fit is good there
- Not the topic of these talks on time series analysis

##  HRF is the idealization of measurable FMRI signal change responding to a single activation cycle (up and down) from a stimulus in a voxel

-  **HRF** is the idealization of measurable FMRI signal change responding to a single activation cycle (up and down) from a stimulus in a voxel

- Response to brief activation (< 1 s):
-  delay of 1-2 s
-  rise time of 4-5 s
-  fall time of 4-6 s
-  model equation:
-  *h*(*t *) is signal change *t* seconds ***after*** activation

_**1 Brief Activation (Event)**_

- Hemodynamic Response Function (HRF)

## Linearity of HRF

-  Multiple activation cycles in a voxel, closer in time than duration of HRF:
-  Assume that overlapping responses add

-  Linearity is a pretty good assumption
-  But not apparently perfect — about 90% correct
-  Nevertheless, is widely taken to be true and is the basis for the “general linear model” (GLM) in FMRI analysis

_**3 Brief Activations**_

## Linearity and Extended Activation

-  Extended activation, as in a block-trial experiment:
-  HRF accumulates over its duration ( 10 s)

_**2 Long Activations (Blocks)**_

-  _**Black**_ curve = response to a single brief stimulus
-  _**Red**_ curve = activation intervals
-  _**Green**_ curve = summed up HRFs from activations
-  Block-trials have larger BOLD signal changes than event-related experiments

## Convolution Signal Model

-  FMRI signal model (in each voxel) is taken as sum of the individual trial HRFs (assumed equal)
-  Stimulus timing is assumed known (or measured)
-  Resulting time series (in **blue**) are called the ***convolution*** of the HRF with stimulus timing
-  Finding HRF = “deconvolution”
-  AFNI code = _**3dDeconvolve**_
-  Convolution models only the FMRI signal **changes**

**22 s**

**120 s**

-  Real data starts at and
-   returns to a nonzero,
-   slowly drifting baseline

##  Assume a fixed shape h(t ) for the HRF
 e.g., h(t ) = t 8.6 exp(-t /0.547)   [MS Cohen, 1997]
 Convolve with stimulus timing to get ideal response (temporal pattern)
 Assume a form for the baseline (data without activation)
 e.g.,  a + bt   for a constant plus a linear trend
 In each voxel, fit data Z(t ) to a curve of the form
         Z(t )  a + b  t +   r (t )
 a, b,  are unknown values to be found in each voxel
 a, b are “nuisance” parameters
  is amplitude of r (t ) in data = “how much” BOLD
In this model, each stimulus assumed to get same BOLD response — in shape and in amplitude

-  e.g., *h*(*t** *) = *t** *8.6 exp(-*t** */0.547)   [MS Cohen, 1997]
-  Convolve with stimulus timing to get ideal response (temporal pattern)
-  Assume a form for the baseline (data without activation)
-  e.g.,  *a *+* b****t*   for a constant plus a linear trend
-  In each voxel, fit data *Z*(*t** *) to a curve of the form
*         **Z*(*t** *)* **** a *+* b **** t* + ******* **** r** *(*t** *)
-  *a*, *b*, ****** are unknown values to be found in each voxel
-  *a*, *b* are “nuisance” parameters
-  ****** is amplitude of *r** *(*t** *) in data = “how much” BOLD
- In this model, each stimulus assumed to get same BOLD response — in shape and in amplitude

- Simple Regression Models

- The signal model*!*

## Simple Regression: Sample Fits

**Constant baseline**: ***a***

**Quadratic baseline: *****a****** *****+*****b************t****** *****+*****c************t****** *****2**

-  Necessary baseline model complexity depends on duration of **continuous** imaging — e.g., 1 parameter per 150 seconds

## Duration of Stimuli - Important Caveats

-  Slow baseline drift (time scale 100 s and longer) makes doing FMRI with _long duration_ stimuli difficult
-  Learning experiment: where the task is done continuously for 15 minutes and the subject is scanned to find parts of the brain that adapt during this time interval
-  Pharmaceutical challenge: where the subject is given some psychoactive drug whose action plays out over 10+ minutes (e.g., cocaine, ethanol)
-  Multiple very _short duration_ stimuli that are also very close in time to each other are very hard to tell apart, since their HRFs will have 90-95% overlap
-  Binocular rivalry, where percept switches  0.5 s

## 900 s

- _900 s_

- Is it Baseline Drift?  Or Activation?

- Is it one extended activation?
- Or four overlapping activations?

- 4 stimulus times (_**waver**_ + _**1dplot**_)

- _19 s_

- Individual HRFs

- Sum of HRFs

- not real data!

## Multiple Stimuli = Multiple Regressors

-  Usually have more than one class of stimulus or activation in an experiment
-  e.g., want to see size of “**face activation**” vis-à-vis “**house activation**”; or, “**what**” vs. “**where**” activity
-  Need to model each separate class of stimulus with a separate response function *r*1(*t** *), *r*2(*t** *), *r*3(*t** *), ....
-  Each *r**j*(*t** *) is based on the stimulus timing for activity in class number *j*
* *Calculate a ***j* amplitude = amount of *r**j** *(*t** *) in voxel data time series *Z*(*t** *) = average BOLD for stim class #*j*
-  **Contrast** *** *s to see which voxels have differential activation levels under different stimulus conditions
-  e.g., statistical test on the question ********1****–**********2**** = 0 ?**

## Multiple Stimuli - Important Caveat

-  In AFNI: do _***not***_ model baseline (“control”) condition
-  e.g., “rest”, visual fixation, high-low tone discrimination, or some other simple task
-  FMRI can only measure _***changes***_ in MR signal levels between tasks
-  So you need some simple-ish task to serve as a reference point
-  The baseline model (e.g., ***a****** *****+**** *****b****** **********t****** ***) takes care of the signal level to which the MR signal returns when the “active” tasks are turned off
-  Modeling the reference task explicitly would be redundant (or “collinear”, to anticipate a forthcoming concept)

## Multiple Stimuli - Experiment Design

-  How many distinct stimuli do you need in each class?  Our rough recommendations:
-  _Short event-related designs_: at least 25 events in each stimulus class (spread across multiple imaging runs) — and more is better
-  _Block designs_: at least 5 blocks in each stimulus class — 10 would be better
-  While we’re on the subject: _**How many subjects?**_
-  Several independent studies agree that 20-25 subjects in each category are needed for highly reliable results
-  This number is more than has usually been the custom in FMRI-based studies**!!**

##  IM = Individual Modulation
Compute separate amplitude of HRF for each event
Instead of the standard computation of the average amplitude of all responses to multiple stimuli in the same class
Response amplitudes (s) for each individual block/event will be highly noisy
Can’t use individual activation maps for much
Must pool the computed s in some further statistical analysis (t-test via 3dttest? inter-voxel correlations in the s? correlate s with something?)
Further description and examples given in the Advanced Topics presentation in this series (afni07_advanced)

- Compute separate amplitude of HRF for each event
- Instead of the standard computation of the average amplitude of all responses to multiple stimuli in the same class
- Response amplitudes (**s) for each individual block/event will be highly noisy
- Can’t use individual activation maps for much
- Must pool the computed **s in some further statistical analysis (*t*-test via **3dttest**? inter-voxel correlations in the **s? correlate **s with something?)
- Further description and examples given in the *Advanced Topics* presentation in this series (afni07_advanced)

_**IM**__** Regression - an Aside**_

## Multiple Regressors: Cartoon Animation

-  _**Red**_ curve = signal model for class #1
-  _**Green**_ curve = signal model for #2
-  _**Blue**_ curve =
-     ********1****#1**+********2****#2**
- where ********1** and ********2** vary from 0.1 to 1.7 in the animation
-  Goal of regression is to find ********1** and ********2** that make the blue curve best fit the data time series
** **_**Gray**_ curve =
**  ****1.5****#1**+**0.6****#2**+**noise**
-  = simulated data

## Multiple Regressors: Collinearity!!

- _Multiple Regressors: Collinearity!!_

_**Green**_ curve = signal model for #1
_**Red**_ curve = signal model for class #2
_**Blue**_ curve = signal model for #3
**Purple** curve =
**    ****#****1**** **+ **#****2**** **+ **#****3**
- which is exactly = 1
-  We cannot — ***in principle or in practice****** ***— distinguish sum of 3 signal models from constant baseline!!

_***No***_ analysis can distinguish the cases
*   **Z*(*t** *)=10+  5**#****1**                       and
*   **Z*(*t** *)=  0+15**#****1**+10**#****2**+10**#****3**
- and an infinity of other possibilities

- Collinear designs
- are _**bad**__** **__**bad**__** **__**bad**_!

## Multiple Regressors: Near Collinearity

_**Red**_ curve = signal model for class #1
_**Green**_ curve = signal model for #2
_**Blue**_ curve =
-     ********1****#****1**+**(1–**********1****)****#****2**
- where ******1** varies randomly  from 0.0 to 1.0 in animation
_**Gray**_ curve =
** **** ****0.66****#****1**+**0.33****#****2**
-  = simulated data *with no noise*
-  Lots of different combinations of **#****1 **and **#****2** are decent fits to gray curve

- Stimuli are too close in time to distinguish
- response **#****1 **from **#****2**, considering noise

**Red** & **Green** stimuli average 2 s apart

## The Geometry of Collinearity - 1

- z1

- z2

**Basis vectors**

**r****2**

**r****1**

**z=Data value**
**=**** ****1.3********r****1****+1.1********r****2**

- }

- Non-collinear
- (well-posed)

- z1

- z2

**r****2**

**r****1**

**z=Data value**
**=**** ********1.8********r****1****+7.2********r****2**

- }

- Near-collinear
- (ill-posed)

-  Trying to fit data as a sum of basis vectors that are nearly parallel doesn’t work well: solutions can be huge
-  Exactly parallel basis vectors would be impossible:
-  Determinant of matrix to invert would be zero

## The Geometry of Collinearity - 2

- z1

- z2

**Basis**
**vectors**

**r****2**

**r****1**

**z=Data value**
**=**** ****1.7********r****1****+2.8********r****2**
** ****=**** 5.1********r****2**** ******** ****3.1********r****3**
**    ****=**** an ******** of other**
**combinations**

- }

- Multi-collinear
**=** more than
- one solution
- fits the data
**=** over-determined

-  Trying to fit data with too many regressors (basis vectors) doesn’t work: no unique solution

**r****3**

## Equations: Notation

-  Will approximately follow notation of manual for the AFNI program _**3dDeconvolve**_
-  Time: continuous in reality, but in steps in the data
-  Functions of continuous time are written like *f** *(*t** *)
-  Functions of discrete time expressed like     * * where *n** *= 0,1,2,... and TR=time step
-  Usually use subscript notion *f**n* as shorthand
-  Collection of numbers assembled in a column is a 			***vector*** and is printed in boldface:

## Equations: Single Response Function

- _Equations: Single Response Function_

-  In each voxel, fit data *Z**n* to a curve of the form
*       **Z**n** **** a *+* b****t**n* + ***r**n*   for *n** *=* *0,1,...,*N** *–1 (*N** *=* *# time pts)
* **a*, *b*, ** are unknown parameters to be calculated in each voxel
-  *a*,*b* are “nuisance” baseline parameters
-  ** is amplitude of *r** *(*t** *) in data = “how much” BOLD
-  Baseline model should be more complicated for long (> 150 s) continuous imaging runs:
-  150 < *T* < 300 s: *a*+*b**t** *+*c**t** *2
-  Longer:               *a*+*b**t** *+*c**t** *2 + *T** */150 low frequency components
**3dDeconvolve** actually uses Legendre polynomials for baseline
-  Often, also include as extra baseline components the estimated subject head movement time series, in order to remove residual contamination from such artifacts (will see example of this later)

- 1 param per 150 s

## Equations: Multiple Response Functions

-  In each voxel, fit data *Z**n* to a curve of the form
*         *
-  ***j* is amplitude in data of *r**n*(*j *)=*r**j *(*t**n*) ; i.e., “how much” of the *j** *th response function is in the data time series
-  In simple regression, each *r**j*(*t *) is derived directly from stimulus timing ***and*** user-chosen HRF model
-  In terms of stimulus times:
- Where      is the *k**th* stimulus time in the *j**th* stimulus class
- These times are input using the **-stim_times** option to program **3dDeconvolve**

## Equations: Matrix-Vector Form

-  Express ***known*** data vector as a sum of ***known*** columns with ***unknown*** coefficents:

_*or*_

_*or*_

**‘********’*** *means “least squares”

- the “design” matrix; AKA **X**

**z** depends on the voxel; **R** doesn’t

-  Const  baseline
-  Linear trend
-  Response to stim#1
-  Response to stim#2

## Visualizing the R Matrix

-  Can graph columns (program **1dplot**)
-  But might have 20-50 columns
-  Can plot columns on a grayscale (program **1dgrayplot** or **3dDeconvolve**** ****-xjpeg**)
-  Easier way to show many columns
-  In this plot, darker bars means larger numbers

*constant baseline: column **#**1*

*linear trend: column **#**2*

*response to stim A: column **#**3*

*response to stim B column **#**4*

## Solving zR for 

-  Number of equations = number of time points
-  100s per run, but perhaps 1000s per subject
-  Number of unknowns usually in range 5–50
-  Least squares solution:
-      denotes an *estimate* of the true (unknown)
-  From    , calculate            as the *fitted model*
-             is the ***residual time series*** = noise (we hope)
- Statistics measure how much each regressor helps reduce residuals
* *Collinearity: when matrix         can’t be inverted
* *Near collinearity: when inverse exists but is huge

## Simple Regression: Recapitulation

-  Choose HRF model ***h*****(*****t*****)** [AKA *fixed-model regression*]
-  Build model responses ***r******n*****(*****t*****)** to each stimulus class
-  Using ***h*****(*****t*****)** and the stimulus timing
-  Choose baseline model time series
-  Constant + linear + quadratic (+ movement?)
-  Assemble model and baseline time series into the columns of the **R** matrix
-  For each voxel time series **z**, solve **z********R******** for
-  **Individual subject maps**: Test the coefficients in that you care about for statistical significance
-  **Group maps**: Transform the coefficients in    that you care about to Talairach space, and perform statistics on the collection of ** **  values across subjects