A general statistical analysis for fMRI data

1
A general statistical analysis for fMRI data

• Keith Worsley12, Chuanhong Liao1, John Aston123,
• Jean-Baptiste Poline4, Gary Duncan5, Vali Petre2,
  Alan Evans2
• 1Department of Mathematics and Statistics, McGill
  University,
• 2Brain Imaging Centre, Montreal Neurological
  Institute,
• 3Imperial College, London,
• 4Service Hospitalier Frédéric Joliot, CEA, Orsay,
  
• 5Centre de Recherche en Sciences Neurologiques,
  Université de Montréal

2
Choices

• Time domain / frequency domain?
• AR / ARMA / state space models?
• Linear / non-linear time series model?
• Fixed HRF / estimated HRF?
• Voxel / local / global parameters?
• Fixed effects / random effects?
• Frequentist / Bayesian?

3
More importantly ...

• Fast execution / slow execution?
• Matlab / C?
• Script (batch) / GUI?
• Lazy / hard working ?
• Why not just use SPM?
• Develop new ideas ...

4
Aim Simple, general, valid, robust, fast
analysis of fMRI data

• Linear model
• ?
  ?
• Yt (stimulust HRF) b driftt c errort
• AR(p) errors
• ? ?
  ?
• errort a1 errort-1 ap errort-p s WNt

unknown parameters
5
MATLAB reads MINC or analyze format
(www/math.mcgill.ca/keith/fmristat)

• FMRIDESIGN Sets up stimulus, convolves it with
  the HRF and its derivatives (for estimating
  delay).
• FMRILM Fits model, estimates effects (contrasts
  in the magnitudes, b), standard errors, T and F
  statistics.
• MULTISTAT Combines effects from separate
  runs/sessions/subjects in a hierarchical fixed /
  random effects analysis.
• TSTAT_THRESHOLD Uses random field theory /
  Bonferroni to find thresholds for corrected
  P-values for peaks and clusters of T and F maps.

6
(No Transcript)
7
Example Pain perception
(c) Response, x(t) sampled at the slice
acquisition times every 3 seconds
Time, t (seconds)
8
(No Transcript)
9
First step estimate the autocorrelation
?

• AR(1) model errort a1 errort-1 s WNt
• Fit the linear model using least squares
• errort Yt fitted Yt
• â1 Correlation ( errort , errort-1)
• Estimating errorts changes their correlation
  structure slightly, so â1 is slightly biased
• Raw autocorrelation Smoothed 15mm Bias
  corrected â1
• 
  -0.05 0

10
Second step refit the linear model
Pre-whiten Yt Yt â1 Yt-1, then fit using
least squares Effect hot warm
Sd of effect
T statistic Effect / Sd
11
Higher order AR model? Try AR(4)
â1
â2
0
â3
â4
0
0
AR(1) seems to be adequate
12
has no effect on the T statistics
AR(1) AR(2)
AR(4)
13
(No Transcript)
14
Results from 4 runs on the same subject
Run 1 Run 2 Run 3 Run 4
Effect Ei
Sd Si
T stat Ei / Si
15
MULTISTAT combines effects from different
runs/sessions/subjects

• Ei effect for run/session/subject i
• Si standard error of effect
• Mixed effects model
• Ei covariatesi c Si WNiF ? WNiR

from FMRILM
?
?
Usually 1, but could add group, treatment,
age, sex, ...
Random effect, due to variability from run to run
Fixed effects error, due to variability within
the same run
16
REML estimation using the EM algorithm

• Slow to converge (10 iterations by default).
• Stable (maintains estimate ?2 gt 0 ), but
• ?2 biased if ?2 (random effect) is small, so
• Re-parametrise the variance model
• Var(Ei) Si2 ?2
• (Si2 minj Sj2) (?2
  minj Sj2)
• Si2
  ?2
• ?2 ?2 minj Sj2 (less biased estimate)

?
?


17
Run 1 Run 2 Run 3 Run 4
MULTISTAT
Effect Ei
very noisy sd
Sd Si
and Tgt15.96 for Plt0.05 (corrected)
T stat Ei / Si
so no response is detected
18
Solution Spatial regularization of the sd

• 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
Random effects sd (3 df)
Fixed effects sd (448 df)
20
Effective df

• dfratio dfrandom(2 1)
• 1 1 1
• dfeff dfratio dffixed
• e.g. dfrandom 3, dffixed 112, FWHMdata 6mm
• FWHMratio (mm) 0 5 10 15 20 infinite
• dfeff 3 11 45 112 192
  448

FWHMratio2 3/2 FWHMdata2

Random effects Fixed effects
variability bias
compromise!
21
Final result 15mm smoothing, 112 effective df
Run 1 Run 2 Run 3 Run 4
MULTISTAT
Effect Ei
less noisy sd
Sd Si
and Tgt4.86 for Plt0.05 (corrected)
T stat Ei / Si
and now we can detect a response!
22
(No Transcript)
23
T gt 4.86 (P lt 0.05, corrected)
Tgt4.86
24
Conclusion

• 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!

25
P.S. Estimating the delay of the response

• Delays or latency in the neuronal response are
  modeled as a temporal scale shift in the
  reference HRF
• Fast voxel-wise delay estimator is found by
  adding the derivative of the reference HRF with
  respect to the log scale shift as an extra term
  to the linear model.
• Bias correction using the second derivative.
• Shrunk to the reference delay by a factor of
  1/(11/T2),
• T is the T statistic for the magnitude.

0.6
Delay 4.0 seconds, log scale shift -0.3
0.4
Delay 5.4 seconds, log scale shift 0
(reference hrf, h0)
Delay 7.3 seconds, log scale shift 0.3
0.2
0
-0.2
0
5
10
15
20
25
t (seconds)
26
Delay of the hot stimulus
T stat for magnitude 0 T stat for
delay 5.4 secs
Delay (secs) Sd of
delay (secs)
27
gt4.86 0 5.4s ?0.6s
28
T gt 4.86 (P lt 0.05, corrected)
Delay (secs)
6.5 5 5.5 4 4.5
29
T gt 4.86 (P lt 0.05, corrected)
Delay (secs)
6.5 5 5.5 4 4.5
30

• Comparison
• Different slice acquisition times
• Drift removal
• Temporal correlation
• Estimation of effects
• Rationale
• Random effects
• Map of the delay

• SPM99
• Adds a temporal derivative
• Low frequency cosines (flat at the ends)
• AR(1), global parameter, bias reduction not
  necessary
• Band pass filter, then least-squares, then
  correction for temporal correlation
• More robust,
• but lower df
• No regularization,
• low df
• No

• fmristat
• Shifts the model
• Polynomials
• (free at the ends)
• AR(p), voxel parameters, bias reduction
• Pre-whiten, then least squares (no further
  corrections needed)
• More accurate, higher df
• Regularization, high df
• Yes

31
(No Transcript)
32
(No Transcript)
33
(No Transcript)
34
(No Transcript)
35
References

• http//www.math.mcgill.ca/keith/fmristat
• Worsley et al. (2000). A general statistical
  analysis for fMRI data. NeuroImage, 11S648, and
  submitted.
• Liao et al. (2001). Estimating the delay of the
  fMRI response. NeuroImage, 13S185 and submitted.

36
False Discovery Rate (FDR)

• Benjamini and Hochberg (1995), Journal of the
  Royal Statistical Society
• Benjamini and Yekutieli (2001), Annals of
  Statistics
• Genovese et al. (2001), NeuroImage
• FDR controls the expected proportion of false
  positives amongst the discoveries, whereas
• Bonferroni / random field theory controls the
  probability of any false positives
• No correction controls the proportion of false
  positives in the volume

37
P lt 0.05 (uncorrected), Z gt 1.64 5 of volume is
false
Signal Gaussian white noise
Signal
True
Noise
False
FDR lt 0.05, Z gt 2.82 5 of discoveries is false
P lt 0.05 (corrected), Z gt 4.22 5 probability of
any false
38
Comparison of thresholds

• FDR depends on the ordered P-values
• P1 lt P2 lt lt Pn. To control the FDR at a
  0.05, find
• K max i Pi lt (i/n) a, threshold the
  P-values at PK
• Proportion of true 1 0.1 0.01
  0.001 0.0001
• Threshold Z 1.64 2.56 3.28
  3.88 4.41
• Bonferroni thresholds the P-values at a/n
• Number of voxels 1 10 100 1000
  10000
• Threshold Z 1.64 2.58 3.29
  3.89 4.42
• Random field theory resels volume / FHHM3
• Number of resels 0 1 10
  100 1000
• Threshold Z 1.64 2.82 3.46
  4.09 4.65

39
P lt 0.05 (uncorrected), Z gt 1.64 5 of volume is
false
FDR lt 0.05, Z gt 2.91 5 of discoveries is false
P lt 0.05 (corrected), Z gt 4.86 5 probability of
any false
Which do you prefer?