Document Sample

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 |

OTHER DOCS BY xiaopangnv

How are you planning on using Docstoc?
BUSINESS
PERSONAL

By registering with docstoc.com you agree to our
privacy policy and
terms of service, and to receive content and offer notifications.

Docstoc is the premier online destination to start and grow small businesses. It hosts the best quality and widest selection of professional documents (over 20 million) and resources including expert videos, articles and productivity tools to make every small business better.

Search or Browse for any specific document or resource you need for your business. Or explore our curated resources for Starting a Business, Growing a Business or for Professional Development.

Feel free to Contact Us with any questions you might have.