Unbiased FMRI Analysis: Leave One Subject Out

Neuroimaging researchers are incessantly bedeviled by the problem of biased region of interest (ROI) analysis. One is constantly lured by the siren song of significant results and large effect sizes radiating from the stygian depths of a non-independent ROI; and while one can at times point toward their use of independent ROIs from other studies, there is always the lurking suspicion that the researcher already knew where the activation was before the ROI was chosen. I have witnessed men, otherwise Samsons in the field and Solomons in counsel, who have had their heads shorn by the harlot of biased analysis.

The most straightforward and appropriate way to do this, of course, is with a region defined on a priori assumptions about where your quarry might lie, based on theory or based on the results of other studies. This ensures that any results extracted from that region are uninfluenced by the model used to generate the statistical maps, therefore circumventing the issue of "double-dipping", or circular analyses (see Kriegeskorte et al, 2009). Another method is to use anatomical regions based on atlases, which again should be motivated by theory.

However, there is yet another option that I was unaware of until recently: Leaving one subject out (LOSO). According to this procedure, non-independence can be mitigated by constructing a general linear model (GLM) with every subject in the study except for one; statistics such as beta weights, time courses, etc., can then be extracted from the resulting parametric map for the subject that was left out, as this subject is no longer contributing to the signal observed in the given region. This process is then repeated and the appropriate parameter extracted for each subject. It is unlikely that there will be perfect overlap between all of the subjects included in each LOSO analysis, but if the effect is real and robust, then it should survive each of these non-overlapping regions.

One consideration with this procedure is what threshold to use for each LOSO analysis. One approach is to hold the p-value constant, in which case a higher t-threshold is used for each analysis due to a reduction in the degrees of freedom. The other approach is to hold the t-value constant, leading to a slightly increased p-value. Both approaches are defensible, although if there are wide variations in the ROI results with each approach, one may want to reconsider the reliability of their finding.

More details can be found in the paper by Esterman et al (2010); I hope this provides the necessary edification and enlightenment for those benighted souls wading about in the filth of their own squalor.

AFNI Command of the Week: cdf

Not necessarily a neuroimaging-specific tool, cdf simply converts between p-values and t-statistics (or F-statistics) using the cumulative distribution function. Supply the test that you did, followed by the t-statistic (or p-value) and degrees of freedom, e.g.:

cdf -t2p fitt 3.4 15
p = 0.00396 #A t-statstic of 3.4 with 15 degrees of freedom yields a p-value of 0.00396

cdf -p2t fitt 0.001 30
t = 3.65 #We would need a t-statistic of 3.65 or greater to reach a p-value of 0.001

Degrees of freedom can be found by using 3dinfo on a statistical dataset, and then looking at the value of "statpar" associated with your statistical test of interest. Degrees of freedom can also be calculated as the number of time points minus the number of regressors in your model; for example, if you have 1200 time points and 40 regressors, then the degrees of freedom will be 1200-40 = 1160. In the following X-matrix (generated by the command "aiv X.jpg"), the first 25 columns represent regressors accounting for any drift during that run; the next nine columns are the regressors of interest; and the last six columns are motion regressors.

 

The degrees of freedom in neuroimaging data can be a little tricky to interpret, as FMRI time series are temporally autocorrelated; in other words, the value of one time point can be predicted, to a degree, by neighboring timepoints. Therefore, using ordinary GLM estimation techniques can lead to an inflated degrees of freedom. To rectify this, instead of using 3dDeconvolve, use 3dREMLfit, which will attempt to account for this autocorrelation.

Central Limit Theorem Part 2: Retaliation

Discrete population, with different probabilities associated with different numbers

Sampling distribution of means from the discrete parent population, with a sample size of n=30

In the last post, we left the central limit theorem defined as a normally-distributed sampling distribution of means reflecting the shape of the normally-distributed parent population, but with a smaller spread and less variance. However, what happens when we sample from a non-normal distribution, such as an exponential distribution or a discrete distribution?

As it turns out, the sampling distribution of means is also normal, regardless of the shape of the parent population. This holds for sample sizes of about 30 or more, which is why the central limit theorem is also sometimes referred to as the law of large numbers.

This is shown in the following video, and can be modified with this R script.

The Central Limit Theorem: Part 1

Random sample of numbers from a normal distribution, N ~ (100, 10). Actual normal distribution is superimposed in red.

One fundamental concept for hypothesis testing is something called the Central Limit Theorem. This theorem states that, for large enough sample sizes and for enough samples, we begin to build a sampling distribution that is approximately normal. More importantly, when we build sampling distributions of the means selected from a population, the average mean is identical to the mean of the parent population.

To illustrate this in R, from the parent population we can take random samples of several different sizes - 10, 50, 300 - and plot those samples as a histogram. These samples will roughly follow the shape of the population they were drawn from - in this case, the normal distribution with a mean of 100 and a standard deviation of 10 - and the more observations we have in our sample, the more closely it reflects the actual parent population. Theoretically, if our sample were large enough, it would in essence sample the entire population and therefore be the same as the population.

However, for smaller sample sizes, we can calculate the mean of each sample and then plot that value in a histogram. If we do this enough times, the mean of the sampling distribution has less spread and more tightly clusters around the mean of the parent population. Increasing the sample size does the same thing.

The demo script can be downloaded here; I have basically copied the code from this website, but distilled it into an R script that can be used and modified by anybody.

Brief Overview of Standard Error

As I begin teaching a statistics course next semester, I've been spending the past couple of weeks hitting the books and refreshing my statistical knowledge; however, to my dismay, I remember virtually nothing of what was taught during my salad days of college, when my greatest concern was how fast I could run eight kilometers, and whether there would be enough ice cream left over in the Burton Dining Hall after a late workout. You laugh now, but during certain eras of one's lifetime, there are specific things that take on especial significance, only to be later ridiculed or belittled; as what is important to an adult may seem insignificant to a child, whereas what is a matter of life and death for the child may seem silly to the adult, even though there is, deep down, recognition of the same hopes and fears, the child the father of the man.

In any case, imagine my sphincter-tightening (and subsequent releasing) horror when I realized how little I actually knew, and with what haste I began to relearn the fundamentals; not only in statistics, but in several other related fields, such as biology, physics, eschatology, chemistry, and astrology, which are needed to have any sense about what one is doing when analyzing neuroimaging data. It is one thing to bandy about the usual formulas and fill them in as needed; it is completely another to learn enough jargon so that, even if you still do not understand it, you can use enough impressive-sounding words to allay any fears that you are hopelessly, utterly ignorant. And this, I maintain, is the end of all good education.

I leave you with one of the most famous quotes about the value of education, from George Washington's second inaugural address:

Power flows to the one who knows how. Desire alone is not enough.

More details about standard error can be found in the following video, which features a legit, squeaking chalkboard.

AFNI Command of the Week: 3dcopy

Nothing earth-shattering here: Just a straightforward conversion command to let AFNI import data from different file types, such as .nii or .img/.hdr extensions. However, when I first started out doing FMRI analysis it took me quite a while to realize that this command was out there; so if you're in the same boat, here it is.

For conversion to AFNI's BRIK/HEAD format, simply type in the name of the dataset to be converted, followed by an output name, e.g.:

3dcopy r01.nii r01
----> r01+orig.HEAD r01+orig.BRIK

Datasets can also be converted the other way around, using something like 3dAFNItoNIFTI, 3dAFNItoANALYZE, etc. Using 3dcopy on a dataset already in .BRIK/.HEAD format will simply rename the dataset without affecting the extension, similar to 3drename.

Note that virtually all of AFNI's commands will output data in BRIK/HEAD format no matter what the input dataset, so you probably will not use this command that often. However, if you absolutely, positively, need to have a dataset in AFNI format - and you need it NOW - then 3dcopy is your new best friend. Move over, thirty pack of Keystone Light!

Resting State Functional Connectivity Database

For those of you not in the know, there is a gimundo FMRI database comprising resting state datasets from research laboratories all over the world. (Meaning the United States, Belgium, Ireland, and the Netherlands.) Scanning parameters and sample features are included with each dataset, making it straightforward to download and analyze entire datasets with relative ease. This is my first time doing resting state functional connectivity analysis, as I have never before collected a dataset using this technique, and I will be sure to document my progress as I go along.


Thanks to Omar Maximo, who can dunk a basketball literally AND figuratively.