Classifying brain regions based on neurosynth features

In this example, we will classify two regions of the brain (ROIs) based on neurosynth features (i.e. words) using a linear support vector machine. In particular, we'll classify studies that activate medial motor regions versus those that activate ventromedial PFC.

Using Nifti masks, we will select studies that activiate only each of the two regions, train a classifier to classify between the two sets of studies, and use cross-validation to estimate the classifier's performance.

First, we need to load a dataset. Let's assume you have already generated one and saved it.

In [1]:
from neurosynth.base.dataset import Dataset
dataset = Dataset.load('../data/pickled.pkl')

We will be using the function "classify_regions" in neurosynth.analysis. Let's take a look at the argument structure:

classify_regions only has two required arguments: a neurosynth dataset and a list of Nifit binary masks to use for selecting studies to classify. Let's give it our dataset and the path to two masks:

In [2]:
from neurosynth.analysis import classify
results = classify.classify_regions(dataset, ["../neurosynth/tests/data/medial_motor.nii.gz", "../neurosynth/tests/data/vmPFC.nii.gz"])

classify_regions will by default return a Python dictionary with various elements of the classification results.

First, we can look at the number of studies selected for each mask:

In [3]:
results['n']
Out[3]:
{0: 2164, 1: 388}

This means there are 2164 studies in the first class (medial_motor) and 388 in the second (vmPFC). Next, lets see how well we classifier by looking at the score, which is by default overall accuracy.

In [4]:
results['score']
Out[4]:
0.74216300940438873

Hey! Not bad, pretty high accuracy. Before we get too excited, let's look a little deeper into classify_region's default behavior by looking at the argument structure:

classify_regions(dataset, masks, method='SVM', threshold=0.001, remove_overlap=True, regularization='scale', output='summary', studies=None, features=None, class_weight=True, classifier=None, cross_val=None):

Let's focus on the first few arguments: by default, a support vector machine is used as the classifier method, studies that activate more than 0.1% of voxels within a mask are included, studies that activate both regions are removed, the feature vectors are regularized to having unit variance and a summary dictionary is returned.

Let's try changing the method. Importantly, we need a baseline to compare to. For that, we can use a Dummy classifier. Dummy classifier use very simple strategies to classify (such as picking the most frequent class) and serve as a good baseline that is specific to your data. We can try this type of classifier by specifying method = "Dummy"

In [5]:
results = classify.classify_regions(dataset, ["../neurosynth/tests/data/medial_motor.nii.gz", "../neurosynth/tests/data/vmPFC.nii.gz"], method="Dummy")
results['score']
Out[5]:
0.74020376175548597

Uh oh, well now it looks like our classifier is not doing much better than our dummy classifier. This means we're not getting much meaningful signal over and above that one class appears more often than the other.

Maybe we need to be more selective in selecting studies by raising the threshold:

In [6]:
classify.classify_regions(dataset, ["../neurosynth/tests/data/medial_motor.nii.gz", "../neurosynth/tests/data/vmPFC.nii.gz"], threshold=0.1)
Out[6]:
{'n': {0: 875, 1: 262}, 'score': 0.69570051890289097}

Our accuracy is now 69% and we only have around 1100 total observations. Let's compare to a Dummy classifier to make sense of this

In [7]:
classify.classify_regions(dataset, ["../neurosynth/tests/data/medial_motor.nii.gz", "../neurosynth/tests/data/vmPFC.nii.gz"], threshold=0.1, method="Dummy")
Out[7]:
{'n': {0: 875, 1: 262}, 'score': 0.57870953792933033}

It looks like now we're classifier a bit better with a SVM than with a Dummy classifier. Great! To really get at this we'd need to run a permutation test to see if they are significantly different but for now, let's rejoice!