Docstoc

boston

Document Sample
boston Powered By Docstoc
					 A general statistical analysis for
           fMRI data
  Keith Worsley12, Chuanhong Liao1, John Aston123,
  Jean-Baptiste Poline4, Gary Duncan5, Vali Petre2,
           Frank Morales6, 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,
                6Cuban Neuroscience Centre
                                                                                        First scan of fMRI data
            fMRI data: 120 scans, 3 scans each of                                                                    1000
            hot, rest, warm, rest, hot, rest, …
                                                                                                                     800
                                 (a) Highly significant correlation
            960
                                                                                                    +                600
fMRI data




                                                                           hot
            940                                                            rest                     +                400
                                                                           warm                         +
            920
                  0   50   100       150       200       250       300   350      400                                200
            940
                                  (b) No significant correlation
fMRI data




                                                                           hot                                       0
            920                                                            rest
                                                                           warm
            900                                                                         Z statistic for hot - warm
                  0   50   100       150       200       250       300   350      400
            850
                                                                                                                     6
fMRI data




                                                                                                                     4
            840
                                             (c) Drift                                                               2
            830
                  0   50   100       150      200     250          300   350      400               +                0
                                       Time, t (seconds)                                            +
                                                                                                        +            -2
                                    Z = (effect hot – warm) / S.d.                                                   -4
                                      ~ N(0,1) if no effect
                                                                                                                     -6
FMRISTAT: Simple, general, valid,
 robust, fast analysis of fMRI data
 • Linear model:
                        ?          ?
 Yt = (stimulust * HRF) b + driftt c + errort
                    unknown parameters
 • AR(p) errors:
          ?                 ?             ?
 errort = a1 errort-1 + … + ap errort-p + s WNt
       FMRIDESIGN example: Pain perception
           (a) Stimulus, s(t): alternating hot and warm stimuli on forearm, separated by rest (9 seconds each).
  2


  1


  0


  -1
       0      50           100            150                200            250            300           350        400
               (b) Hemodynamic response function, h(t): difference of two gamma densities (Glover, 1999)

0.4

0.2

  0

-0.2
       0      50              100             150              200              250           300             350   400
                        (c) Response, x(t): sampled at the slice acquisition times every 3 seconds
  2


  1


  0


  -1
       0      50             100              150             200             250             300             350   400
                                                      Time, t (seconds)
FMRILM step 1: estimate temporal correlation
                            ?
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 errort’s changes their correlation structure
  slightly, so â1 is slightly biased.
• Bias correction is very quick and effective:
Raw autocorrelation   Smoothed 15mm         Bias corrected â1
                                                                  0.4
                             ~ -0.05                  ~0          0.3

                                                                  0.2

                                                                  0.1

                                                                  0

                                                                  -0.1
    FMRILM step 2: refit the linear model
Pre-whiten: Yt* = Yt – â1 Yt-1, then fit using least squares:
         Effect: hot – warm             Sd of effect
                                   6                   2.5

                                   4                   2
                                   2
                                                       1.5
                                   0
                                                       1
                                   -2

                                   -4                  0.5

                                   -6                  0


       T statistic = Effect / Sd
                                   6

                                   4

                                   2
                                        T > 4.90
                                   0    (P < 0.05,
                                   -2   corrected)
                                   -4

                                   -6
Higher order AR model? Try AR(4):
       â1      0.4
                           â2     0.4

               0.3                0.3

               0.2                0.2

               0.1
                          ~0      0.1

               0                  0

               -0.1               -0.1


       â3                  â4
               0.4                0.4

               0.3                0.3

               0.2                0.2
      ~0       0.1
                          ~0      0.1

               0                  0

               -0.1               -0.1

     AR(1) seems to be adequate
… has no effect on the T statistics:
      AR(1)                      AR(2)
                   6                          6

                   4                          4

                   2                          2

                   0                          0

                   -2                         -2

                   -4                         -4

                   -6                         -6


      AR(4)             But ignoring correlation …
                   6                          6

                   4                          4

                   2                          2

                   0                          0

                   -2                         -2

                   -4                         -4

                   -6                         -6

              biases T up ~12%    more false positives
   Results from 4 runs on the same subject
          Run 1   Run 2   Run 3   Run 4
                                          5

Effect
                                          0
  Ei
                                          -5


                                          2.5
                                          2
 Sd                                       1.5
 Si                                       1
                                          0.5
                                          0


                                          5

T stat
                                          0
Ei / Si
                                          -5
    MULTISTAT combines effects from
     different runs/sessions/subjects:
 • Ei = effect for run/session/subject i
 • Si = standard error of effect
                                           from     }
                                         FMRILM
 • Mixed effects model:
                       ?               ?
                                  F +  WN R
      Ei = covariatesi c + Si WNi           i

Usually 1, but
could add group,   ‘Fixed effects’ error,   Random effect,
treatment, age,    due to variability       due to variability
sex, ...           within the same run      from run to run
    REML estimation of the mixed effects
      model using the EM algorithm
• Slow to converge (10 iterations by default).
• Stable (maintains estimate 2 > 0 ), but
                                 ^
• 2 biased if 2 (random effect) is small, so:
  ^
• Re-parametrise the variance model:
                   ?
              2 + 2
  Var(Ei) = Si
          = (Si2 – minj Sj2) + (2 + minj Sj2)
                                     ? 2
          =       Si* 2      +      *
• 2 = *2 – minj Sj2 (less biased estimate)
  ^ ^
 Problem: 4 runs, 3 df for random effects sd             ...
          Run 1   Run 2     Run 3     Run 4    MULTISTAT
                                                             5

Effect
                                                             0
  Ei
                                                             -5

                                          … very noisy sd:
                                                             2
 Sd
 Si                                                          1


                                                             0
                     … and T>15.96 for P<0.05 (corrected):
                                                             5

T stat
                                                             0
Ei / Si
                                                             -5

                               … so no response is detected …
Solution: Spatial regularization of the sd
• Basic idea: increase df by spatial smoothing
(local pooling) of the sd.
• Can’t smooth the random effects sd directly,
- too much anatomical structure.
• Instead,
              random effects sd
sd = smooth
               fixed effects sd   ) fixed effects sd
which removes the anatomical structure
before smoothing.
Random effects sd   Fixed effects sd    Regularized sd
     (3 df)            (448 df)            (112 df)
                                                          4


                                                          3


                                                          2


                                                          1


                                                          0



Random effects sd                      Fixed effects sd
 Fixed effects sd
                     Over runs          Over subjects
                                                          3


              Smooth                                      2
               15mm            ~1                   ~3
                                                          1

                           ~1.6
                                                          0
Effective df depends on the smoothing
                        (
         dfratio = dfrandom 2
                               FWHMratio2
                               FWHMdata   2+1   )   3/2


                     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
             Random effects        Fixed effects
               variability             bias
                           compromise!
Final result: 15mm smoothing, 112 effective df …
          Run 1   Run 2      Run 3     Run 4   MULTISTAT
                                                               5

Effect
                                                               0
  Ei
                                                               -5
                                            … less noisy sd:
                                                               2
 Sd
 Si                                                            1


                                                               0
                      … and T>4.90 for P<0.05 (corrected):
                                                               5

T stat                                                         0
Ei / Si
                                                               -5
                          … and now we can detect a response!
Conjunction: All Ti > threshold = Min Ti > threshold
        ‘Minimum of Ti’          ‘Average of Ti’
                          6                        6
                          4                        4
                          2                        2
                          0                        0
                          -2                       -2
                          -4                       -4
                          -6                       -6
         For P=0.05,              For P=0.05,
       threshold = 1.82   1
                                threshold = 4.90   1

                          0.8                      0.8

                          0.6                      0.6

                          0.4                      0.4

                          0.2                      0.2
       Efficiency = 82%
                          0                        0
       If the conjunction is significant,
       does it mean that all effects > 0?
• Problem: for the conjunction of 20 effects, the threshold
  can be negative!?!?!
• Reason: significance is based on the wrong null
  hypothesis, namely: all effects = 0
• Correct null hypothesis is: at least one effect = 0.
  Unfortunately the P-value depends on the unknown > 0
  effects …
• If the effects are random, all effects > 0 is meaningless.
  The only parameter is the (single) population effect, so
  that the conjunction just tests if population effect > 0.
• P-values now depend on the random effects sd, not the
  fixed effects sd. But the minimum (i.e. the conjunction) is
  less efficient (sensitive) than the average (the usual test).
 FWHM – the local smoothness of the noise
• Used by STAT_THRESHOLD to find the P-value of local maxima and
  the spatial extent of clusters of voxels above a threshold.
• u = normalised residuals from linear model = residuals / sd
  ·
• u = vector of spatial derivatives of u
            ·
• λ = |Var(u)|1/2 (mm-3)
• FWHM = (4 log 2)1/2 λ-1/3 (mm)
  (If residuals are modeled as white noise smoothed with a Gaussian kernel,
  this would be its FWHM).
• λ and FWHM are corrected for low df and large voxel size so they are
  approximately unbiased.
• For a search region S, the number of ‘resolution elements’ is
  Resels(S) = Vol(S) AvgS(FWHM-3) = Vol(S) AvgS(λ) (4 log 2)-3/2
• For local maxima in S, P_value = Resels(S) x (function of threshold).
• For a cluster C, P-value depends on Resels(C) instead of Vol(C), so that
  clusters in smooth regions are less significant.
• Need a correction for the randomness of λ and FWHM - depends on df .
• Correction is more important for small clusters C than for large search
  regions S.
…. FWHM depends on the spatial correlation between neighbours
   FWHM (mm) over scans (448 df)      FWHM (mm) over runs (3 df)
                              20                               20

                              15                               15
     Resels=1.90
     P=0.007                  10                               10


             Resels=0.57      5                                5

             P=0.387          0                                0


   smoothed FWHM over runs         smoothed (runs FWHM / scans FWHM)
                              20                               1.5

                              15

                              10                               1

                              5

                              0                                0.5
T > 4.90
 T>4.86
(P < 0.05,
corrected)
  Smooth the data before analysis?
• Temporal smoothing or low-pass filtering is used by
  SPM’99 to validate a global AR(1) model. For our local
  AR(p) model, it is not necessary (but ~ harmless).

• Spatial smoothing is used by SPM’99 to validate
  random field theory. Can be harmful for focal signals.
  Should fix the theory! STAT_THRESHOLD uses the
  better of the Bonferroni or the random field theory.

• A better reason for spatial smoothing is greater
  detectability of extensive activation: choose the FWHM
  to match the activation (e.g. 10mm FWHM for 10mm
  activations) – or try a range of FWHM’s i.e. scale space
  – but thresholds are higher …
         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
  Signal + Gaussian            P < 0.05 (uncorrected), Z > 1.64
     white noise                   5% of volume is false +
                          4                                4


                          2                                2

        Signal            0
                                         True +            0


                          -2                               -2
        Noise                            False +
                          -4                               -4

  FDR < 0.05, Z > 2.82         P < 0.05 (corrected), Z > 4.22
5% of discoveries is false +   5% probability of any false +
                          4                                4


                          2                                2


                          0                                0


                          -2                               -2


                          -4                               -4
          Comparison of thresholds
• FDR depends on the ordered P-values (not smoothness):
  P1 < P2 < … < Pn. To control the FDR at a = 0.05, find
  K = max {i : Pi < (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
P < 0.05 (uncorrected), Z > 1.64
    5% of volume is false +


                                     FDR < 0.05, Z > 2.91
                                   5% of discoveries is false +



                                                                  P < 0.05 (corrected), Z > 4.86
                                                                  5% probability of any false +




                    Which do you prefer?
 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
  0.6


  0.4
             basis1                              basis2
  0.2
                        HRF
   0


 -0.2


 -0.4
    -5           0         5            10            15   20         25
                                   t (seconds)

               shift

         HRF(t + shift) ~ basis1(t) w1(shift) + basis2(t) w2(shift)

• Convolve bases with the stimulus, then add to the linear model
3                                 • Fit linear model, estimate w1 and w2

2                w2 / w1          • Equate w2 / w1 to estimates, then
                                  solve for shift (Hensen et al., 2002)
1      w1
                                  • To reduce bias when the magnitude
                                  is small, use
0

       w2                                  shift / (1 + 1/T2)
-1

                                  where T = w1 / Sd(w1) is the T statistic
-2                                for the magnitude

                                  • Shrinks shift to 0 where there is little
-3
  -5             0            5   evidence for a response.
            shift (seconds)
Delay of the hot stimulus (= shift + 5.4 sec)
      T stat for magnitude          T stat for shift
                             5                         5



                             0                         0



                             -5                        -5

         Delay (secs)             Sd of delay (secs)
                             10                        5
                             8                         4
                             6                         3
                             4                         2
                             2                         1
                             0                         0
Varying the delay and dispersion of the reference HRF
                                         T stat for magnitude                 T stat for shift
                                                                5                                 5
                                   10.4                               10.4
  Reference dispersion (seconds)




                                   5.2                          0     5.2                         0


                                   2.6                                2.6
                                                                -5                                -5
                                           2.4   5.4   8.4                    2.4    5.4    8.4

                                            Delay (secs)                    Sd of delay (secs)
                                                                10                                5
                                   10.4                         8     10.4                        4

                                                                6                                 3
                                   5.2                                5.2
                                                                4                                 2

                                   2.6                          2     2.6                         1

                                                               0                                  0
                                           2.4   5.4   8.4                   2.4     5.4    8.4
                                                        Reference delay (seconds)
T > 4.86
(P < 0.05,
corrected)

             Delay
             (secs)
                6.5

                5

                5.5

                4

                4.5
T > 4.86
(P < 0.05,
corrected)


             Delay
             (secs)
                6.5

                5

                5.5

                4

                4.5
 EFFICIENCY for optimum block design
                                               Sd of hot stimulus                     Sd of hot-warm
                                          20                               0.5   20                                 0.5

                                                                           0.4                                      0.4
                                          15                                     15
Magnitude                                               Optimum            0.3                                      0.3
                                          10            design                   10
          InterStimulus Interval (secs)




                                                       X                   0.2                Optimum               0.2
                                          5                                      5
                                                                           0.1                design                0.1
                                          0                                0     0            X                     0
                                                  5     10   15       20                5     10   15       20
                                                                       (secs)                                    (secs)
                                          20                              1      20                                 1

                                                                           0.8                                      0.8
                                          15                                     15
  Delay                                               Optimum
                                                                           0.6                                      0.6
                                          10                                     10         Optimum
                                                      design               0.4                                      0.4
                                                                                            design
                                          5       X                              5
                                                                           0.2          X                           0.2
                                                (Not enough signal)                   (Not enough signal)
                                                                           0     0                                  0
                                                  5     10   15       20                5     10   15       20
                                                                  Stimulus Duration (secs)
EFFICIENCY for optimum event design
                                    0.5
                                              (Not             ____ magnitudes          uniform . . . . . . . . .
                                   0.45       enough                                    random .. . ... .. .
                                                               ……. delays
                                               signal)                                  concentrated :
                                    0.4
  Sd of effect (secs for delays)




                                   0.35

                                    0.3

                                   0.25

                                    0.2

                                   0.15

                                    0.1

                                   0.05

                                     0
                                          5                  10                    15                           20
                                                         Average time between events (secs)
            How many subjects?
• Variance =       sdrun2          sdsess2     sdsubj2
               nrun nsess nsubj + nsess nsubj + nsubj

• The largest portion of variance comes from the
  last stage, i.e. combining over subjects.

• If you want to optimize total scanner time, take
  more subjects, rather than more scans per subject.

• What you do at early stages doesn’t matter very
  much - any reasonable design will do …
Comparison SPM’99:                            FMRISTAT:
• Different slice • Adds a temporal           •Shifts the model
acquisition times: derivative
• Drift removal: • Low frequency cosines • Polynomials
                   (flat at the ends)         (free at the ends)
• Temporal         • AR(1), global parameter, • AR(p), voxel
correlation:       bias reduction not         parameters, bias
                   necessary                  reduction
• Estimation of    • Band pass filter, then   • Pre-whiten, then least
effects:           least-squares, then        squares (no further
                   correction for temporal    corrections needed)
                   correlation
• Rationale:       • More robust, low df      • More efficient, high df
• Random effects: • No regularization, low df • Regularization, high df
• FWHM:            • Global, ~ OK for local • Local, is OK for local
                   maxima, but not clusters maxima and clusters
• Map of delay: • No                          • Yes
                 References
• Worsley et al. (2002). A general statistical
  analysis for fMRI data. NeuroImage, 15:1-15.

• Liao et al. (2002). Estimating the delay of the
  fMRI response. NeuroImage, 16:593-606.

• http://www.math.mcgill.ca/keith/fmristat
  - 200K of MATLAB code
  - fully worked example
          Functional connectivity
• Measured by the correlation between residuals at
  every pair of voxels (6D data!)
     Activation only            Correlation only
       Voxel 2    ++              Voxel 2
                   +                             +
                  +++
                                              ++
                  Voxel 1                   +  Voxel 1
                                        +
                                    +
•   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
 |Correlations| > 0.7,
 P<10-10 (corrected)




   First Principal
Component > threshold

				
DOCUMENT INFO
Shared By:
Categories:
Tags:
Stats:
views:0
posted:10/27/2012
language:English
pages:46