Andy's Brain Blog Needs Your Help!

To be more specific, I (Andy) need your help; but what's good for Andy's Brain Blog is good for America - and you're all patriots, right?

As I mentioned before, I am currently applying for jobs and putting together my research and teaching portfolios, playing up all the sexy studies currently in the works, and what I plan to do for the next few years; how I can attract students to the university, students to the department, secure money, funding, recognition, and all that good stuff necessary for the vitality of the university.

However, as all of you know, this right here is one of my dearest, most ambitious projects - to make statistics, neuroimaging, and computational modeling available to everyone in straightforward, simple terms. To use online repositories to get at data unavailable to the majority of smaller, liberal arts institutions, so that students from all parts of the globe, researchers anywhere, and quite literally anyone with a computer can get a piece of the action. To make the information in dense, unreadable technical manuals accessible and easy to understand through hands-on, no-nonsense tutorials. And - perhaps most importantly - I wear a suit when I do it.

I want to continue doing this. I want to continue building, expanding, and teaching, both here and in the classroom. I will not rest: I will drink life to the lees. Heck, maybe I'll even drink the lees too.

But to do that I need your help.

Through both the comments here, on the YouTube channel, and in private correspondence, I've talked with many researchers, students, and professors around the country and around the world. Most of you I've never seen, but I've had the privilege to help out professors and scholars all the way from Italy to China; I've freelanced for PET researchers at Michigan State, schizophrenia experimenters at Indiana, designed experiments for primates in New York. The AFNI tutorials created here have been used as class material at the University of Pittsburgh, and my code for Hodgkin-Huxley simulations have been used for demonstrations at Claremont McKenna College in California. My recipe for homemade granola is used by several hungry researchers to keep them going throughout those long afternoons. The list could go on.

What I ask for is if you have ever used the materials here in an educational setting, be it for the researchers in your lab or the students in your classroom, please let me know by sending an email to ajahn [at] indiana [dot] edu. I am trying to compile a list of where it is used, to demonstrate its use and effectiveness.

Lastly - allow me to get real here for a moment - I've thought of all of you, each person that I've conversed with or replied to or Skyped with, as my children. And not as "children" in the sense of putting you all down on my tax returns as dependents; although, if the IRS isn't too strict about metaphors, I think I could get away with that. No, I've thought of you as real children: Shorter than I am, a little incoherent at times maybe, often trying to get my attention, and playing with Legos and Gak.

Regardless, it's been one of my great pleasures to help you all out. Now get away from that electrical socket.

Bayesian Inference, Step 2: Running Bayesian Inference with Your Data

If you've had the opportunity to install R Studio and JAGS, and to download the associated programs needed to run Bayesian inference, it is a small step to actually getting your parameter estimates. Two of the programs from the BEST folder - BESTExample.R, and BEST1G.R - allow you to run independent-samples and one-sample t-tests, respectively. All that's required is changing the input in the "y1" and "y2" strings with your data, separate each observation with a comma, and then run the program. Other options can be changed, such as which value you are comparing against in the one-sample t-test, but everything else can essentially remain the same.

I realize it's been a while between posts, but right now I'm currently in the process of applying for jobs; this should start to pick up again in mid-October once the deadlines pass, but in the meantime, wish me luck!

Bayesian Inference, Step 1: Installing JAGS On Your Machine

Common complaint: "Bayesian analysis is too hard! Also, I have kidney stones."
Solution: Make Bayesian analysis accessible and efficient through freeware that anyone can use!

These days, advances in technology, computers, and lithotripsy have made Bayesian analysis easy to implement on any personal computer. All it requires is a couple of programs and a library of scripts to run the actual process of Bayesian inference; all that needs to be supplied by you, the user, is the data you have collected. Conceptually, this is no more difficult then entering in data into SAS or SPSS, and, I would argue, is easier in practice.

This can be done in R, statistical software that can interface with a variety of user-created packages. You can download one such package, JAGS, to do the MCMC sampling for building up distributions of parameter estimates, and then use those parameter estimates to brag to your friends about how you've "Gone Bayes."

All of the software and steps you need to install R, JAGS, and rjags (a program allowing JAGS to talk to R) can be found on John Kruschke's website here. Once you have that, it's simply a matter of entering in your own data, and letting the program do the nitty-gritty for you.

Gone for a Few Days

Going to be out of the state for a few days; a few of you have written questions to me about technical issues or methods, and I'll try to get back to you as soon as I can!

In the meantime, as is custom whenever I go on vacation, here is another late masterwork by Brahms. Always a mysterious figure, one of the greatest enigmas surrounding Brahms's life was his relationship with Clara Schumann, piano virtuoso and widow of the famous composer Robert Schumann. Scholars have speculated for decades what happened between them - because, obviously, knowing whether they hooked up or not is essential for how we listen to and appreciate the music of Brahms, and such a line of inquiry is a valid funding opportunity for government grants.

Fortunately, however, I recently came across one of Brahms's lost journals in a second-hand bookstore; and the secrets disclosed within will finally give us an answer. Let's take a look:


November 4th, 1876: Cloudy today. Composed first symphony.

August 9th, 1879: Went clubbing with Liszt and Wagner. Hooked up with Clara Schumann. No doubt about it.

October 24th, 1884: Note to self: Find out why I'm writing my journal in English.

April 3rd, 1897: Died from cancer.

Bayesian Inference Web App That Explains (Nearly) Everything

For those still struggling to understand the concepts of Bayesian inference - in other words, all of you - there is a web app developed by Rasmus Baath which allows one to see the process unfolding in real time. Similar to an independent-samples t-test, we are trying to estimate a parameter of the mean difference between the populations that the samples were drawn from; however, the Bayesian approach offers much richer information about a distribution of parameter values, and, more importantly, which ones are more credible than others, given the data that has been collected.

The web app is simple to use: You input data values for your two groups, specify the number of samples and burn-in samples (although the defaults for this are fine), and then hit the Start button. The MCMC chain begins sampling the posterior distribution, which builds up a distribution of credible parameter values, and 95% of the distribution containing the parameter values with the most credibility is labeled as the highest density interval, or HDI. This same procedure is applied to all of the other parameters informed by the data, including normality, effect size, and individual group means and standard deviations, which provides a much richer set of information than null hypothesis significance testing (NHST).

Because of this I plan to start incorporating more Bayesian statistics into my posts, and also because I believe it will overtake, replace, and destroy NHST as the dominant statistical method in the next ten years, burning its crops, sowing salt in its fields, looting its stores, stampeding its women, and ravishing its cattle. All of this, coming from someone slated to teach the traditional NHST approach next semester; which, understandably, has me conflicted. On the one hand, I feel pressured to do to whatever is most popular at the moment; but on the other, I am an opportunistic coward.

In any case, Bayesian inference is becoming a more tractable technique, thanks to programs that interface with the statistics package R, such such JAGS. Using this to estimate parameters for region of interest data, I think, will be a good first step for introducing Bayesian methods to neuroimagers.

Parametric Modulation with SPM: Why Order Matters

Those who are the youngest children in the family are keenly aware of the disadvantages of being last in the pecking order; the quality of toys, clothes, and nutrition is often best for the eldest children in the family, and gradually tapers off as you reach the last kids, who are usually left with what are called, in parent euphemisms, "hand-me-downs."

In reality, these are things that should have been discarded several years ago due to being broken, unsafe, containing asbestos, or clearly being overused, such as pairs of underwear that have several holes in them so large that you're unsure of where exactly your legs go.

In SPM land, a similar thing happens to parametric modulators, with later modulators getting the dregs of the variance not wanted by the first modulators. Think of variance as toys and cake, and parametric modulators as children; everybody wants the most variance, but only the first modulators get the lion's share of it; the last modulators get whatever wasn't given to the first modulators.

The reason this happens is because when GLMs are set up, SPM uses a command, spm_orth, to orthogonalize the modulators. The first modulator is essentially untouched, but all other modulators coming after it are orthogonalized with respect to the first; and this proceeds in a stepwise function, with the second modulator getting whatever was unassigned to the first, the third getting whatever was unassigned to the second, and so on.

(You can see this for yourself by doing a couple of quick simulations in Matlab. Type x1=rand(10,1), and x2=rand(10,1), and then combine the two into a matrix with xMat = [x1 x2]. Type spm_orth(xMat), and notice that the first column is unchanged, while the second is greatly altered. You can then flip it around by creating another matrix, xMatRev = [x2 x1], and using spm_orth again to see what effect this change has.)

Be aware of this, as this can dramatically alter your results depending on how your modulators were ordered. (If they were completely orthogonal by design, then the spm_orth command should have virtually no impact; but this scenario is rare.) This came to my attention when looking at a contrast of a parametric modulator, which was one of three assigned to a main regressor. In one GLM, the modulator was last, and correspondingly explained almost no variance in the data; while in another GLM, the modulator was first, and loaded on much more of the brain. The following two figures show the difference, using exactly the same contrast and correction threshold, but simply reversing the order of the modulators:

Results window when parametric modulator was the last one in order.

Results window when parametric modulator was the first one in order.

This can be dealt with in a couple of ways:

1. Comment out the spm_orth call in two scripts used in setting up the GLM: line 229 in spm_get_ons.m, and lines 285-287 in spm_fMRI_design.m. However, this might affect people who actually do want the orthogonalization on, and also this increases the probability that you will screw something up. Also note that doing this means that the modulators will not get de-meaned when entered into the GLM.
2. Create your own modulators by de-meaning the regressor, inserting your value at the TR where the modulator occurred, then convolving this with an HRF (usually the one in SPM.xBF.bf). If needed, downsample the output to match the original TR.

It is well to be aware of all this, not only for your own contrasts, but to see how this is done in other papers. I recall reading one a long time ago - I can't remember the name - where several different models were used, in particular different models where the modulators were entered in different orders. Not surprisingly, the results were vastly different, although why this was the case wasn't clearly explained. Just something to keep in mind, as both a reader and a reviewer.

I should also point out that this orthogonalization does not occur in AFNI; the order of modulators has no influence on how they are estimated. A couple of simulations should help illustrate this:

Design matrices from AFNI with simulated data (one regressor, two parametric modulators.) Left: PM1 entered first, PM2 second. Right: Order reversed. The values are the same in both windows, just transposed to illustrate which one was first in the model.

More information on the topic can be found here, and detailed instructions on how to construct your own modulators can be found on Tor Wager's lab website.

DTI Analysis: Soup to Nuts Playlist

Instead of going through each DTI analysis step individually, I've collated everything into a Youtube playlist down below. Just remember that we are using data from the FSL practical course here, and also remember that suing somebody for giving out bad advice, although it is admittedly an easy way to become fantastically wealthy, won't necessarily make you happier.

In any case, just to briefly go over the rest of the steps: After correcting for magnetic field distortions and eddy currents, tensors are fitted using the dtifit command (or simply going through the FDT menu in the FSL interface). Once this has been done for each subject, a series of TBSS tools are used, each one prefixed by "tbss"; for example, tbss_1_preproc, tbss_2_reg, and so on. (You can find all of these in the $FSLDIR/bin directory, and if you have a good handle on Unix programming, you can inspect the code yourself.) After you have run all of those for your dataset, you set up the appropriate experimental design and contrasts, and use the randomise tool to perform statistics in your tractography mask.

Keep in mind that this is just beginner's material; and that if you were to compare your DTI competence to dating, if would be like you were still in that awkward teenager phase, unable to talk to anybody or make eye contact.

However, much of this material has already been covered in other blogs and online resources, provided by several other highly talented scientists and researchers, and - as much as I constantly fantasize about establishing a monopoly over neuroimaging information - there is no practical way to destroy them all.

Therefore, instead of posting redundant information, I highly recommend checking out an ongoing step-by-step series on TORTOISE by Peter Molfese, which you can compare to your results with FSL, and another blog dedicated to diffusion imaging, diffusion-imaging.com. The latter site covers virtually every piece of diffusion imaging software and pipeline out there, and is a good place to start.