The%20statistical%20analysis%20of%20fMRI%20data

1
The statistical analysis of fMRI data

• Keith Worsley12, Chuanhong Liao1, John Aston123,
  Jean-Baptiste Poline4, Gary Duncan5, Vali Petre2,
  Frank Morales6, Alan Evans2, Tom Nichols7, Satoru
  Hayasaki8, Jonathan Taylor9
• 1Department of Mathematics and Statistics, McGill
  University,
• 2Brain Imaging Centre, Montreal Neurological
  Institute,
• 3Statistica Sinica, Taipei,
• 4Neurospin, CEA, Orsay,
• 5Centre de Recherche en Sciences Neurologiques,
  Université de Montréal,
• 6Cuban Neuroscience Centre
• 7GlaxoSmithKline and FMRIB, Oxford
• 8Wake Forest University
• 9Université de Montréal and Stanford

2
(No Transcript)
3
Before you start PCA of time ? space
1 exclude first frames
2 drift
3 long-range correlation or anatomical effect
remove by converting to of brain
4 signal?
4
Bad design 2 mins rest 2 mins Mozart 2 mins
Eminem 2 mins James Brown
5
Rest Mozart Eminem J. Brown
Temporal components (sd, variance
explained)
Period 5.2 16.1
15.6 11.6 seconds
1
0.41, 17
2
0.31, 9.5
Component
3
0.24, 5.6
0
50
100
150
200
Frame
Spatial components

1
1
0.5
Component
2
0
-0.5
3

-1
0
2
4
6
8
10
12
14
16
18
Slice (0 based)
6
Effect of stimulus on brain response
Modeled by convolving the stimulus with the
hemodynamic response function
Stimulus is delayed and dispersed by 6s
7
fMRI data, pain experiment, one slice
T (hot warm effect) / S.d. t110 if no
effect
8
How fMRI differs from other repeated measures data

• Many reps (200 time points)
• Few subjects (15)
• Df within subjects is high, so not worth pooling
  sd across subjects
• Df between subjects low, so use spatial smoothing
  to boost df
• Data sets are huge 4GB, not easy to use
  statistics packages such as R

9
FMRISTAT (Matlab) / BRAINSTAT (Python)statistica
l analysis strategy

• Analyse each voxel separately
• Borrow strength from neighbours when needed
• Break up analysis into stages
• 1st level analyse each time series separately
• 2nd level combine 1st level results over runs
• 3rd level combine 2nd level results over
  subjects
• Cut corners do a reasonable analysis in a
  reasonable time (or else no one will use it!)

10
(No Transcript)
11
1st level Linear model with AR(p) errors

• Data
• Yt fMRI data at time t
• xt (responses,1, t, t2, t3, ) to allow for
  drift
• Model
• Yt xtß et
• et a1et-1 apet-p sF?t, ?t N(0,1)
  i.i.d.
• Fit in 2 stages
• 1st pass fit by least squares, find residuals,
  estimate AR parameters a1 ap
• 2nd pass whiten data, re-fit by least squares

12
Higher levelsMixed effects model

• Data
• Ei effect (contrast in ß) from previous level
• Si sd of effect from previous level
• zi (1, treatment, group, gender, )
• Model
• Ei zi? SieiF sReiR (Si high df, so
  assumed fixed)
• eiF N(0,1) i.i.d. fixed effects error
• eiR N(0,1) i.i.d. random effects error
• Fit by ReML
• Use EM for stability, 10 iterations

13
Where we use spatial information

• 1st level smooth AR parameters to lower
  variability and increase df
• Higher levels smooth Random / Fixed effects sd
  ratio to lower variability and increase df
• Final level use random field theory to correct
  for multiple comparisons

14
1st level Autocorrelation

• AR(1) model et a1 et-1 sF?t
• Fit the linear model using least squares
• et Yt Yt
• â1 Correlation (et , et-1)
• Estimating et changes their correlation structure
  slightly, so â1 is slightly biased
• Raw autocorrelation Smoothed 12.4mm Bias
  corrected â1
• 
  -0.05 0

15
How much smoothing?

• Variability in
• â lowers df
• Df depends
• on contrast
• Smoothing â brings df back up

dfâ dfresidual(2
1) 1 1 2 acor(contrast
of data)2 dfeff dfresidual
dfâ
FWHMâ2 3/2 FWHMdata2

Hot-warm stimulus
Hot stimulus
FWHMdata 8.79
Residual df 110
Residual df 110
100
100
Target 100 df
Target 100 df
Contrast of data, acor 0.79
50
50
Contrast of data, acor 0.61
dfeff
dfeff
0
0
0
10
20
30
0
10
20
30
FWHM 10.3mm
FWHM 12.4mm
FWHMâ
FWHMâ
16
Higher order AR model? Try AR(3)
â
â
â
1
2
3
0.3
0.2
AR(1) seems to be adequate
0.1
0
has little effect on the T statistics
-0.1
AR(1), df100
AR(2), df99
AR(3), df98
5
0
-5
17
2nd level 4 runs, 3 df for random effects sd
very noisy sd
and Tgt15.96 for Plt0.05 (corrected)
so no response is detected
18
Solution Spatial smoothing of the sd ratio

• Basic idea increase df by spatial smoothing
  (local pooling) of the sd.
• Cant smooth the random effects sd directly, -
  too much anatomical structure.
• Instead,
• random effects sd
• fixed effects sd
• which removes the anatomical structure before
  smoothing.

? )
sd smooth ?
fixed effects sd
19

Average Si
?
Random effects sd, 3 df
Fixed effects sd, 440 df
Mixed effects sd, 100 df
0.2
0.15
0.1
0.05
0
divide
multiply
Smoothed sd ratio
Random sd / fixed sd
1.5
random effect, sd ratio 1.3
1
0.5
20
How much smoothing?

• dfratio dfrandom(2 1)
• 1 1 1
• dfeff dfratio dffixed

FWHMratio2 3/2 FWHMdata2
dfrandom 3, dffixed 4 ? 110
440, FWHMdata 8mm

fixed effects analysis, dfeff 440
dfeff
FWHM 19mm
Target 100 df
random effects analysis, dfeff 3
FWHMratio
21
Final result 19mm smoothing, 100 df
less noisy sd
and Tgt4.93 for Plt0.05 (corrected)
and now we can detect a response!
22
Final level Multiple comparisons correction
Bonferroni
4.7
4.6
4.5
True
T, 10 df
4.4
Random Field Theory
4.3
T, 20 df
Gaussianized threshold
Discrete Local Maxima (DLM)
4.2

4.1
Gaussian
4
3.9
3.8
3.7
0
1
2
3
4
5
6
7
8
9
10
FWHM of smoothing kernel (voxels)
In between use Discrete Local Maxima (DLM)
High FWHM use Random Field Theory
Low FWHM use Bonferroni
23
0.12
Gaussian
T, 20 df
T, 10 df
0.1
Random Field Theory
Bonferroni
0.08
0.06
P-value
0.04
Discrete Local Maxima
True
DLM can ½ P-value when FWHM 3 voxels
0.02
0
0
1
2
3
4
5
6
7
8
9
10
FWHM of smoothing kernel (voxels)
In between use Discrete Local Maxima (DLM)
High FWHM use Random Field Theory
Low FWHM use Bonferroni
24
(No Transcript)
25
Example single run, hot-warm
Detected by BON and DLM but not by RFT
Detected by DLM, but not by BON or RFT
26
FWHM the local smoothness of the noise
voxel size (1 correlation)1/2
FWHM (2 log
2)1/2
(If the noise is modeled as white noise smoothed
with a Gaussian kernel, this would be its FWHM)
Volume FWHM3
P-values depend on Resels
Resels
Local maximum T 4.5
Clusters above t 3.0, search volume resels 500
0.1
0.1
0.08
0.08
0.06
0.06
P value of local max
P value of cluster
0.04
0.04
0.02
0.02
0
0
0
500
1000
0
0.5
1
1.5
2
Resels of search volume
Resels of cluster
27
Non-isotropic data (spatially varying FWHM)
FWHM (mm) of scans (110 df)
FWHM (mm) of effects (3 df)
20
20
Cluster Resels1.90 P0.007
15
15
10
10
Cluster Resels0.57 P0.387
5
5
0
0

• fMRI data is smoother in GM than WM
• Has little effect on peak P-values
• use average FWHM inside search region, but
• Has a big effect on cluster P-values
• smooth regions ? big clusters,
  
• rough regions ? small clusters, so
• Replace cluster volume by cluster resels volume
  / FWHM3

28
Estimating the delay of the response

• Delay or latency to the peak of the HRF is
  approximated by a linear combination of two
  optimally chosen basis functions

delay
basis1
basis2
HRF
shift

• HRF(t shift) basis1(t) w1(shift)
  basis2(t) w2(shift)
• Convolve bases with the stimulus, then add to
  the linear model

29

• Fit linear model,
• estimate w1 and w2
• Equate w2 / w1 to estimates, then solve for
  shift (Hensen et al., 2002)
• To reduce bias when the magnitude is small, use
• shift / (1 1/T2)
• where T w1 / Sd(w1) is the T statistic for the
  magnitude
• Shrinks shift to 0 where there is little
  evidence for a response.

w2 / w1
w1
w2
30
Shift of the hot stimulus
T stat for magnitude T stat for
shift
Shift (secs) Sd of shift
(secs)
31
Shift of the hot stimulus
T stat for magnitude T stat for
shift
Tgt4
T2
Shift (secs) Sd of shift
(secs)
1 sec
/- 0.5 sec
32
Combining shifts of the hot stimulus (Contours
are T stat for magnitude gt 4)
33
Shift of the hot stimulus
Shift (secs)
T stat for magnitude gt 4.93
34
Contrasts in the data used for effects
2
Hot, Sd 0.16
Warm, Sd 0.16
9 sec blocks, 9 sec gaps
1
0
-1
0
50
100
150
200
250
300
350
Hot-warm, Sd 0.19
Time (secs)
2
Hot, Sd 0.28
Warm, Sd 0.43
90 sec blocks, 90 sec gaps
Only using data near block transitions
1
0
Ignoring data in the middle of blocks
-1
0
50
100
150
200
250
300
350
Hot-warm, Sd 0.55
Time (secs)
35
Optimum block design
Sd of hot stimulus
Sd of hot-warm
0.5
0.5
20
20
0.4
0.4
15
15
Best design
0.3
0.3
Magnitude
10
10
Best design
0.2
0.2
X
5
5
0.1
0.1
0
X
0
0
0
5
10
15
20
5
10
15
20
Gap (secs)
(secs)
(secs)
1
1
20
20
0.8
0.8
15
15
0.6
0.6
Best design X
Best design X
Delay
10
10
0.4
0.4
5
5
0.2
0.2
(Not enough signal)
(Not enough signal)
0
0
0
5
10
15
20
5
10
15
20
Block (secs)
36
Optimum event design
0.5
(Not enough signal)
uniform . . . . . . . . .
____ magnitudes . delays
random .. . ... .. .
concentrated
0.4
Sd of effect (secs for delays)
0.3
0.2
12 secs best for magnitudes
0.1
0
7 secs best for delays
5
10
15
20
Average time between events (secs)
37
How many subjects?

• Largest portion of variance comes from the last
  stage i.e. combining over subjects
• sdrun2 sdsess2
  sdsubj2
• nrun nsess nsubj nsess nsubj
  nsubj
• If you want to optimize total scanner time, take
  more subjects.
• What you do at early stages doesnt matter very
  much!

38
Features special to FMRISTAT / BRAINSTAT

• Bias correction for AR coefficients
• Df boosting due to smoothing
• AR coefficients
• random/fixed effects variance
• P-value adjustment for
• peaks due to small FWHM (DLM)
• clusters due to spatially varying FWHM
• Delays analysed the same way as magnitudes
• Sd of effects before collecting data

39
Our entry in the Functional Imaging Analysis
contest

• Jonathan Taylor
• Stanford
• Keith Worsley
• McGill

40
Why a Functional Imaging Analysis Contest (FIAC)?

• Competing packages produce slightly different
  results, which is correct?
• Simulated data?
• Real data, compare analyses
• Contest session at 2005 Human Brain Map
  conference
• 9 entrants
• Results in a special issue of Human Brain Mapping
  in May, 2006

41
The main participants

• SPM (Statistical Parametric Mapping, 1993),
  University College, London, SAS, (MATLAB)
• AFNI (1995), NIH, more display and manipulation,
  not much stats (C)
• FSL (2000), Oxford, the upstart (C)
• .
• FMRISTAT (2001), McGill, stats only (MATLAB)
• BRAINSTAT (2005), Stanford/McGill, Python version
  of FMRISTAT

42
FIAC paradigm

• 16 subjects
• 4 runs per subject
• 2 runs event design
• 2 runs block design
• 4 conditions per run
• Same sentence, same speaker
• Same sentence, different speaker
• Different sentence, same speaker
• Different sentence, different speaker
• 3T, 191 frames, TR2.5s

43
Response

• Events
• Blocks

Beginning of block/run
44
Design matrix for block expt

• B1, B2 are basis functions for magnitude and
  delay

45
1st level analysis

• Motion and slice time correction (using FSL)
• 5 conditions
• Smoothing of temporal autocorrelation to control
  the effective df

3 contrasts Beginning of block/run Same sent, same speak Same sent, diff speak Diff sent, same speak Diff sent, diff speak
Sentence 0 -0.5 -0.5 0.5 0.5
Speaker 0 -0.5 0.5 -0.5 0.5
Interaction 0 1 -1 -1 1
46
Efficiency

• Sd of contrasts (lower is better) for a single
  run, assuming additivity of responses
• For the magnitudes, event and block have similar
  efficiency
• For the delays, event is much better.

47
2nd level analysis

• Analyse events and blocks separately
• Register contrasts to Talairach (using FSL)
• Bad registration on 2 subjects - dropped
• Combine 2 runs using fixed FX
• Combine remaining 14 subjects using random FX
• 3 contrasts event/block magnitude/delay 12
• Threshold using best of Bonferroni, random field
  theory, and discrete local maxima (new!)

3rd level analysis
48
Part of slice z -2 mm
49
(No Transcript)
50
(No Transcript)
51
(No Transcript)
52
(No Transcript)
53
Event
Block
Magnitude
Delay
54
Delays in different same sentence

• Events 0.140.04s Blocks 1.190.23s
• Both significant, Plt0.05 (corrected) (!?!)
• Answer take a look at blocks

Greater magnitude
Different sentence (sustained interest)
Best fitting block
Same sentence (lose interest)
Greater delay
55
FMRISTAT/
SPM BRAINSTAT
56

• Magnitude increase for
• Sentence, Event
• Sentence, Block
• Sentence, Combined
• Speaker, Combined
• at (-54,-14,-2)

57

• Magnitude decrease for
• Sentence, Block
• Sentence, Combined
• at (-54,-54,40)

58

• Delay increase for
• Sentence, Event
• at (58,-18,2)
• inside the region where all
• conditions are activated

59
Conclusions

• Greater BOLD response for
• different same sentences (1.080.16)
• different same speaker (0.470.0.8)
• Greater latency for
• different same sentences (0.1480.035 secs)

60
The main effects of sentence repetition (in red)
and of speaker repetition (in blue). 1 Meriaux
et al, Madic 2 Goebel et al, Brain voyager 3
Beckman et al, FSL 4 Dehaene-Lambertz et al,
SPM2.
Fmristat/ Brainstat combined block and event,
threshold at Tgt5.67, Plt0.05.
61
Functional connectivity

• Measured by the correlation between residuals at
  every pair of voxels (6D data!)
• Local maxima are larger than all 12 neighbours
• P-value can be calculated using random field
  theory
• Good at detecting focal connectivity, but
• PCA of residuals x voxels is better at detecting
  large regions of co-correlated voxels

Activation only
Correlation only
Voxel 2


Voxel 2







Voxel 1
Voxel 1



62
Correlations gt 0.7, Plt10-10 (corrected)
First Principal Component gt threshold