1 / 39

A general statistical analysis for fMRI data

A general statistical analysis for fMRI data. Keith Worsley 12 , Chuanhong Liao 1 , John Aston 123 , Jean-Baptiste Poline 4 , Gary Duncan 5 , Vali Petre 2 , Alan Evans 2 1 Department of Mathematics and Statistics, McGill University, 2 Brain Imaging Centre, Montreal Neurological Institute,

fadey
Download Presentation

A general statistical analysis for fMRI data

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  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. (a) Stimulus, s(t): alternating hot and warm stimuli on forearm, separated by rest (9 seconds each). 2 warm warm hot hot 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 2 1 0 -1 0 50 100 150 200 250 300 350 400 Example: Pain perception (c) Response, x(t): sampled at the slice acquisition times every 3 seconds Time, t (seconds)

  7. First step: estimate the autocorrelation ? AR(1) model: errort = a1errort-1 + s WNt • Fit the linear model using least squares • errort = Yt – fitted Yt • â1 = Correlation ( errort , errort-1) • Estimating errort’schanges their correlation structure slightly, so â1 is slightly biased: Raw autocorrelation Smoothed 15mm Bias corrected â1 ~ -0.05 ~ 0

  8. T > 4.86 (P < 0.05, corrected) 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

  9. Higher order AR model? Try AR(4): â1 â2 ~ 0 â3 â4 ~ 0 ~ 0 AR(1) seems to be adequate

  10. But ignoring correlation … biases T up ~12%  more false positives … has no effect on the T statistics: AR(1) AR(2) AR(4)

  11. Results from 4 runs on the same subject Run 1 Run 2 Run 3 Run 4 Effect Ei Sd Si T stat Ei / Si

  12. 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, ... ‘Fixed effects’ error, due to variability within the same run Random effect, due to variability from run to run

  13. REML estimation 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: Var(Ei) = Si2 + 2 = (Si2 –minj Sj2) + (2 +minj Sj2) = Si*2 + *2 • 2 = *2 –minj Sj2 (less biased estimate) ^ ^ ? ? ^ ^

  14. Problem: 4 runs, 3 df for random effects sd  ... Run 1 Run 2 Run 3 Run 4 MULTISTAT Effect Ei … very noisy sd: Sd Si … and T>15.96 for P<0.05 (corrected): T stat Ei / Si … so no response is detected …

  15. 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 • fixed effects sd • which removes the anatomical structure before smoothing.  ) sd= smooth  fixed effects sd

  16. Regularized sd (112 df) Random effects sd Fixed effects sd  Fixed effects sd Over scans Over subjects Smooth 15mm ~1 ~3 ~1.6 Random effects sd (3 df) Fixed effects sd (448 df)

  17. Effective df FWHMratio2 3/2 FWHMdata2 dfratio = dfrandom(2 + 1) 1 1 1 dfeffdfratiodffixed e.g. dfrandom = 3, dffixed = 112, FWHMdata = 6mm: FWHMratio (mm) 0 5 10 15 20 infinite dfeff 3 11 45 112 192 448 = + Random effectsFixed effects variability bias compromise!

  18. Final result: 15mm smoothing, 112 effective df … Run 1 Run 2 Run 3 Run 4 MULTISTAT Effect Ei … less noisy sd: Sd Si … and T>4.86 for P<0.05 (corrected): T stat Ei / Si … and now we can detect a response!

  19. T > 4.86 (P < 0.05, corrected) T>4.86

  20. + + Conclusion • Largest portion of variance comes from the last stage i.e. combining over subjects: sdrun2sdsess2sdsubj2 nrunnsessnsubj nsessnsubj nsubj • If you want to optimize total scanner time, take more subjects. • What you do at early stages doesn’t matter very much!

  21. P.S. Estimating the delay of the response 0.6 • 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/(1+1/T2), T is the T statistic for the magnitude. 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)

  22. Delay of the hot stimulus T stat for magnitude = 0 T stat for delay = 5.4 secs Delay (secs) Sd of delay (secs)

  23. Varying the delay of the reference HRF T stat for mag T stat for delay Delay Sd of delay Ref. delay = 4.0 ~5.4s Ref. delay = 5.4 Ref. delay = 7.3 ~5.4s >4.86 ~0 ~5.4s 0.6s

  24. T > 4.86 (P < 0.05, corrected) Delay (secs) 6.5 5 5.5 4 4.5

  25. T > 4.86 (P < 0.05, corrected) Delay (secs) 6.5 5 5.5 4 4.5

  26. Comparison: • Different slice acquisition times: • Drift removal: • Temporal correlation: • Estimation of effects: • Rationale: • Random effects: • Map of the delay: • SPM’99: • 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

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

  28. 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

  29. Signal + Gaussian white noise P < 0.05 (uncorrected), Z > 1.64 5% of volume is false + Signal True + Noise False + FDR < 0.05, Z > 2.82 5% of discoveries is false + P < 0.05 (corrected), Z > 4.22 5% probability of any false +

  30. Comparison of thresholds • FDR depends on the ordered P-values: • 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

  31. 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?

More Related