fMRI Analysis with emphasis on the General Linear Model

1. fMRI Analysiswith emphasis on the General Linear Model Jody Culham Brain and Mind Institute Department of Psychology University of Western Ontario http://www.fmri4newbies.com/ Last Update: January 18, 2012 Last Course: Psychology 9223, W2010, University of Western Ontario

2. Part 1 Statistical Intuitions

3. What data do we start with • 12 slices * 64 voxels x 64 voxels = 49,152 voxels • Each voxel has 136 time points (volumes) • Therefore, for each run, we have 6.7 million data points • We often have several runs for each experiment … These #s are from an obsolete scanner. With a modern 3T, we can get 3X the slices

4. The signal is much higher where there is brain, but there’s still noise Slice 9, Voxel 0, 0 Slice 9, Voxel 1, 0 Slice 9, Voxel 22, 7 Even where there’s no brain, there’s noise Slice 9, Voxel 18, 36 Slice 9, Voxel 9, 27 Here’s a couple that sort of show the right pattern but is it “real”? Here’s a voxel that responds well whenever there’s visual stimulation Slice 9, Voxel 13, 41 Slice 9, Voxel 14, 42 Here’s one that responds well whenever there’s intact objects Why do we need stats? • We could, in principle, analyze data by voxel surfing: move the cursor over different areas and see if any of the time courses look interesting

5. Why do we need stats? • Clearly voxel surfing isn’t a viable option. We’d have to do it 49,152 times in this data set and it would require a lot of subjective decisions about whether activation was real • This is why we need statistics • Statistics: • tell us where to look for activation that is related to our paradigm • help us decide how likely it is that activation is “real” The lies and damned lies come in when you write the manuscript

6. SHIFTED CONVOLVED WITH HRF Predicted Responses • fMRI is based on the Blood Oxygenation Level Dependent (BOLD) response • It takes about 5 sec for the blood to catch up with the brain • We can model the predicted activation in one of two ways: • shift the boxcar by approximately 5 seconds (2 images x 2 seconds/image = 4 sec, close enough) • convolve the boxcar with the hemodynamic responseto model the shape of the true function as well as the delay PREDICTED ACTIVATION IN OBJECT AREA PREDICTED ACTIVATION IN VISUAL AREA BOXCAR

7. Is the region truly active? Yes No Type I Error HIT Yes Does our stat test indicate that the region is active? Type II Error Correct Rejection No Types of Errors p value: probability of a Type I error e.g., p <.05 “There is less than a 5% probability that a voxel our stats have declared as “active” is in reality NOT active Slide modified from Duke course

8. Statistical Approaches in a Nutshell t-tests • compare activation levels between two conditions • use a time-shift to account for hemodynamic lag • correlations • model activation and see whether any areas show a similar pattern • Fourier analysis • Do a Fourier analysis to see if there is energy at your paradigm frequency Fourier analysis images from Huettel, Song & McCarthy, 2004, Functional Magnetic Resonance Imaging

9. r = .80 64% of variance p < 10-33 r = .50 25% of variance p < .000001 r = .40 16% of variance p < .000001 r = .24 6% of variance p < .05 Effect of Thresholds r = 0 0% of variance p < 1

10. Complications • Not only is it hard to determine what’s real, but there are all sorts of statistical problems • Potential problems • data may be contaminated by artifacts (e.g., head motion, breathing artifacts) • .05 * 49,152 = 2457 “significant” voxels by chance alone • many assumptions of statistics (adjacent voxels uncorrelated with each other; adjacent time points uncorrelated with one another) are false What’s wrong with these data? r = .24 6% of variance p < .05

11. The General Linear Model (GLM) GLM definition from Huettel et al.: • a class of statistical tests that assume that the experimental data are composed of the linear combination of different model factors, along with uncorrelated noise • Model • statistical model • Linear • things add up sensibly (1+1 = 2) • note that linearity refers to the predictors in the model and not necessarily the BOLD signal • General • many simpler statistical procedures such as correlations, t-tests and ANOVAs are subsumed by the GLM

12. Benefits of the GLM • GLM is an overarching tool that can do anything that the simpler tests do • allows any combination of contrasts (e.g., intact - scrambled, scrambled - baseline), unlike simpler methods (correlations, t-tests, Fourier analyses) • allows more complex designs (e.g., factorial designs) • allows much greater flexibility for combining data within subjects and between subjects • allows comparisons between groups • allows counterbalancing orders within and between subjects • allows modelling of known sources of noise in the data (e.g., error trials, head motion)

13. Part 2 Composition of a Voxel Time Course

14. A Simple Experiment • Lateral Occipital Complex • responds when subject views objects Blank Screen Intact Objects Scrambled Objects TIME One volume (12 slices) every 2 seconds for 272 seconds (4 minutes, 32 seconds) Condition changes every 16 seconds (8 volumes)

15. What’s real? A. C. B. D.

16. What’s real? • I created each of those time courses based by taking the predictor function and adding a variable amount of random noise signal = + noise

17. What’s real? Which of the data sets below is more convincing?

18. Formal Statistics • Formal statistics are just doing what your eyeball test of significance did • Estimate how likely it is that the signal is real given how noisy the data is • confidence: how likely is it that the results could occur purely due to chance? • “p value” = probability value • If “p = .03”, that means there is a .03/1 or 3% chance that the results are bogus • By convention, if the probability that a result could be due to chance is less than 5% (p < .05), we say that result is statistically significant • Significance depends on • signal (differences between conditions) • noise (other variability) • sample size (more time points are more convincing)

19. Let’s create a time course for one LO voxel

20. We’ll begin with activation Response to Intact Objects is 4X greater than Scrambled Objects

21. Then we’ll assume that our modelled activation is off because a transient component

22. Our modelled activation could be off for other reasons All of the following could lead to inaccurate models • different shape of function • different width of function • different latency of function

23. Reminder: Variability of HRF Intersubject variability of HRF in M1 Handwerker et al., 2004, NeuroImage

24. Now let’s add some variability due to head motion

25. …though really motion is more complex • Head motion can be quantified with 6 parameters given in any motion correction algorithm • x translation • y translation • z translation • xy rotation • xz rotation • yz rotation • For simplicity, I’ve only included parameter one in our model • Head motion can lead to other problems not predictable by these parameters

26. Now let’s throw in a pinch of linear drift • linear drift could arise from magnet noise (e.g., parts warm up) or physiological noise (e.g., subject’s head sinks)

27. and then we’ll add a dash of low frequency noise • low frequency noise can arise from magnet noise or physiological noise (e.g., subject’s cycles of alertness/drowsiness) • low frequency noise would occur over a range of frequencies but for simplicity, I’ve only included one frequency (1 cycle per run) here • Linear drift is really just very low frequency noise

28. and our last ingredient… some high frequency noise • high frequency noise can arise from magnet noise or physiological noise (e.g., subject’s breathing rate and heartrate)

29. When we add these all together, we get a realistic time course

30. Part 3 General Linear Model

31. Now let’s be the experimenter • First, we take our time course and normalize it using z scores • z = (x - mean)/SD • normalization leads to data where • mean = zero • SD = 1 Alternative: You can transform the data into % BOLD signal change. This is usually a better approach because it’s not dependent on variance

32. Wake Up!!!!! If you only pay attention to one slide in this lecture, it should be the next one!!!

33. We create a GLM with 2 predictors × 1 = + + × 2 = + fMRI Signal Design Matrix x Betas Residuals “what we CAN explain” “how much of it we CAN explain” “what we CANNOT explain” = x + “our data” Statistical significance is basically a ratio of explained to unexplained variance

34. Intact Predictor Scrambled Predictor Implementation of GLM in SPM • SPM represents time as going down • SPM represents predictors within the design matrix as grayscale plots (where black = low, white = high) over time • GLM includes a constant to take care of the average activation level throughout each run • SPM shows this explicity (BV may not) Many thanks to Øystein Bech Gadmar for creating this figure in SPM  Time

35. Effect of Beta Weights • Adjustments to the beta weights have the effect of raising or lowering the height of the predictor while keeping the shape constant

36. Dynamic Example

37. The beta weight is NOT a correlation • correlations measure goodness of fit regardless of scale • beta weights are a measure of scale small ß small r small ß large r large ß small r large ß large r

38. + + fMRI Signal Residuals We create a GLM with 2 predictors when 1=2 = + when 2=0.5 = Design Matrix x Betas “what we CAN explain” “how much of it we CAN explain” “what we CANNOT explain” = x + “our data” Statistical significance is basically a ratio of explained to unexplained variance

39. Correlated Predictors • Where possible, avoid predictors that are highly correlated with one another • This is why we NEVER include a baseline predictor • baseline predictor is almost completely correlated with the sum of existing predictors + r = -.53 = r = -.53 r = -.95 Two stimulus predictors Baseline predictor

40. Which model accounts for this data? • Because the predictors are highly correlated, the model is overdetermined and you can’t tell which beta combo is best x β = 1 x β = 0 + + OR x β = 1 x β = 0 + + x β = 0 x β = -1

41. Orthogonalizing Regressors

42. We can also examine whether a single predictor is significantly greater than another predictor Contrasts in the GLM • We can examine whether a single predictor is significant (compared to the baseline)

43. Contrasts  “balanced” • Conjunction of contrasts • e.g., (+1 -1 0) AND (+1 0 -1) • (Bio motion - Nonbio motion) AND (Bio motion > control) • more rigorous than balanced contrast • hypothetical (but not actual) conjunction p value = multiple of independent p values • e.g., .01 x .01 = .001

44. A Real Voxel • Here’s the time course from a voxel that was significant in the +Intact -Scrambled comparison

45. Maximizing Your Power signal = + noise • As we saw earlier, the GLM is basically comparing the amount of signal to the amount of noise • How can we improve our stats? • increase signal • decrease noise • increase sample size (keep subject in longer)

46. How to Reduce Noise • If you can’t get rid of an artifact, you can include it as a “predictor of no interest” to soak up variance Example: Some people include predictors from the outcome of motion correction algorithms Corollary: Never leave out predictors for conditions that will affect your data (e.g., error trials) This works best when the motion is uncorrelated with your paradigm (predictors of interest)

47. Reducing Residuals

48. Part 3 Deconvolution of Event-Related Designs Using the GLM

49. Convolution of Single Trials Neuronal Activity BOLD Signal Haemodynamic Function Time Time Slide from Matt Brown

50. Fast fMRI Detection A) BOLD Signal B) Individual Haemodynamic Components C) 2 Predictor Curves for use with GLM (summation of B) Slide from Matt Brown