written by Luke Chang 8/15/2017
This tutorial will provide an introduction to the nltools toolbox, which is loosely based off of Tor Wager's Matlab Toolbox and is the analysis engine that power neuro-learn.org.
In addition, we will try to replicate training the Picture Induced Negative Emotion Signature using rating data stored in neurovault. For details about the methods see our paper.
Here we fetch the affective rating dataset used in Chang et al., 2015 from neurovault. In this dataset there are 182 subjects with approximately 5 separate beta images reflecting varying intensities of their negative affective ratings in response to arousing IAPS images. The data will be downloaded to ~/nilearn_data, and automatically loaded as a Brain_Data() instance. The image metadata will be stored in data.X.
In [4]:
%matplotlib inline
from nltools.datasets import fetch_emotion_ratings
data = fetch_emotion_ratings()
The bulk of the nltools toolbox is built around the Brain_Data() class. This class represents imaging data as a vectorized features by observations matrix. Each image is an observation and each voxel is a feature. The concept behind the class is to have a similar feel to a pandas dataframe, which means that it should feel intuitive to manipulate the data.
In [18]:
data.X.head()
Out[18]:
In [4]:
print(len(data))
Find the dimensions of the data. images x voxels
In [5]:
print(data.shape())
We can use any type of indexing to slice the data such as integers, lists of integers, or boolean.
In [6]:
print(data[[1,6,2]])
Calculate the mean for every voxel over images
In [7]:
data.mean()
Out[7]:
Calculate the standard deviation for every voxel over images
In [8]:
data.std()
Out[8]:
Methods can be chained. Here we get the shape of the mean.
In [9]:
print(data.mean().shape())
Brain_Data instances can be added and subtracted
In [10]:
new = data[1]+data[2]
Brain_Data instances can be manipulated with basic arithmetic operations Here we add 10 to every voxel and scale by 2
In [11]:
data2 = (data+10)*2
Brain_Data instances can be copied
In [12]:
new = data.copy()
Brain_Data instances can be easily converted to nibabel instances, which store the data in a 3D/4D matrix. This is useful for interfacing with other python toolboxes such as nilearn
In [13]:
data.to_nifti()
Out[13]:
Brain_Data instances can be concatenated using the append method
In [14]:
new = new.append(data[4])
Any Brain_Data object can be written out to a nifti file
In [16]:
data[1].write('Tmp_Data.nii.gz')
Images within a Brain_Data() instance are iterable. Here we use a list comprehension to calculate the overall mean across all voxels within an image.
For speed we only calculate the mean for the first 10 images.
In [17]:
[x.mean() for x in data[:10]]
Out[17]:
In [5]:
from nilearn.plotting import plot_glass_brain
plot_glass_brain(data.mean().to_nifti())
Out[5]:
There is also a fast montage plotting method. Here we plot the average image it will render a separate plot for each image. There is a 'limit' flag which allows you to specify the maximum number of images to display.
In [6]:
data.mean().plot()
We also have functions to create a more comprehensive plot from multiple orientations
In [7]:
from nltools.plotting import plotBrain
plotBrain(data.mean())
Running MVPA style analyses using multivariate regression is even easier and faster than univariate methods. All you need to do is specify the algorithm and cross-validation parameters. Currently, we have several different linear algorithms implemented from scikit-learn.
Here is a helpful blog post on different algorithms and reasonable default parameters.
Cross-validation is critical for evaluating the performance of models. We recommend splitting the data into a training and separate holdout test sets. The final holdout set should only be used to test the model once, while the training data can be explored and fit multiple times (though we do recommend cross-validation here too). As an example, we will split the data into train (2/3 of data) and test (1/3 of data) based on a stratifed sampling approach. See details here.
In [15]:
train = data[data.X['Holdout']=='Train']
train.Y = train.X['Rating']
test = data[data.X['Holdout']=='Test']
test.Y = test.X['Rating']
In [19]:
train.X['SubjectID'].head()
Out[19]:
We can now predict the output variable is a dictionary of the most
useful output from the prediction analyses. The predict function runs
the prediction multiple times. One of the iterations uses all of the
data to calculate the weight_map. The other iterations are to estimate
the cross-validated predictive accuracy. For speed, we are only going to estimate the model using all of the data. See examples below for different algorithms and cross-validation schemes.
In [12]:
stats = train.predict(algorithm='ridge',plot=False)
Display the available data in the output dictionary
In [44]:
stats.keys()
Out[44]:
Plot the multivariate weight map
In [64]:
plotBrain(stats['weight_map'])
Once we have a model that we are happy with and performs the best in cross-validation on the training data, we can test it on the holdout test dataset. It is important to mention that this should only happen once and should be the final results reported in the paper. The reason is because it is easy to introduce bias and trick yourself into thinking that the model will generalize better than it should.
This two level cross-validation scheme provides a way to explore and try lots of different algorithms and feature selection methods, without overfitting the data.
Because this is a regression we also need to remember to add the intercept to shift the predicted values.
In [13]:
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
predicted_emotion = test.similarity(stats['weight_map'],'dot_product') + stats['intercept']
dat = pd.DataFrame(data={'SubjectID':test.X['SubjectID'],'Predicted':predicted_emotion,'Rating':test.X['Rating']})
with sns.plotting_context(context='paper',font_scale=2):
sns.factorplot(data=dat,x='Rating',y='Predicted')
We can also do similarity instead of dot product to compare performance of model across data on different scales
In [66]:
predicted_emotion = test.similarity(stats['weight_map'],'correlation')
dat = pd.DataFrame(data={'SubjectID':test.X['SubjectID'],'Predicted':predicted_emotion,'Rating':test.X['Rating']})
with sns.plotting_context(context='paper',font_scale=2):
sns.factorplot(data=dat,x='Rating',y='Predicted')
Now that we have evaluated the sensitivity of the model in capturing affect in the holdout data, we might want to know how specific it is. If it is simply reflecting negative arousal or salience (whatever that is) than it should be confused by other negative arousing stimuli such as thermal pain.
Here we first download a pain data set to assess the specificity from our paper.
In [3]:
from nltools.datasets import fetch_pain
pain = fetch_pain()
pain.Y = pain.X['PainLevel']
subject_id = pain.X['SubjectID']
Now let's test how well the PINES generalizes to the pain test data.
In [67]:
predicted_pain = pain.similarity(stats['weight_map'],'correlation')
pain_dat = pd.DataFrame(data={'SubjectID':pain.X['SubjectID'],'Predicted':predicted_pain,'Rating':pain.X['PainLevel']})
with sns.plotting_context(context='paper',font_scale=2):
sns.factorplot(data=pain_dat,x='Rating',y='Predicted')
Let's put the two datasets on the same plot to show the sensitivity and specificity (at least with respect to pain)
In [68]:
dat['Data']='IAPS'
pain_dat['Data']='Pain'
d = pd.concat([dat,pain_dat])
with sns.plotting_context(context='paper',font_scale=2):
sns.factorplot(data=d,x='Rating',y='Predicted',col='Data')
Feature selection describes the process of deciding which features to include when training the model. Here it is simply, which voxels should we use to train the model?
There are several ways to perform feature selection. Searchlights are a popular approach. I personally have a preference for using parcellation schemes.
Here we download a single 50 parcel map from a forthcoming paper on conducting automated parcellations using neurosynth.
Yarkoni, T., de la Vega, A., & Chang, L.J. (In Prep). Fully automated meta-analytic clustering and decoding of human brain activity
Some of the details can be found here
In [1]:
from nltools.mask import expand_mask, collapse_mask
from nltools.data import Brain_Data
mask = Brain_Data('http://neurovault.org/media/images/2099/Neurosynth%20Parcellation_0.nii.gz')
mask.plot()
Out[1]:
We can easily expand the mask into separate binary masks, which makes it easier to use.
In [52]:
mask_x = expand_mask(mask)
mask_x.plot(limit=3)
First, let's see how well we can predict pain intensity using whole brain
In [53]:
stats = pain.predict(algorithm='ridge', cv_dict={'type': 'kfolds','n_folds': 5,'subject_id':subject_id})
Now let's see how well we can do with a subset of the voxels by only training voxels located in the ACC
You can see we are doing pretty good with only the ACC, but the whole brain is capturing more variance. This implies that information is distributed more broadly than just the ACC.
See Lieberman & Eisenberger (2016) and our rebuttal Wager et al (2016) along with Tal Yarkoni's, Tor Wager, and Alex Shackman's Blogpost and for a more detailed debate about the selectivity of the dACC in processing pain.
In [54]:
masked_stats = pain.apply_mask(mask_x[0]).predict(algorithm='ridge', cv_dict={'type': 'kfolds','n_folds': 5,'subject_id':subject_id})
In [17]:
subject_id = train.X['SubjectID']
svr_stats = train.predict(algorithm='svr',
cv_dict={'type': 'kfolds','n_folds': 5,
'subject_id':subject_id}, **{'kernel':"linear"})
Lasso Regression
In [ ]:
lasso_stats = train.predict(algorithm='lasso',
cv_dict={'type': 'kfolds','n_folds': 5,
'subject_id':subject_id}, **{'alpha':.1})
Principal Components Regression
In [ ]:
pcr_stats = train.predict(algorithm='pcr',
cv_dict={'type': 'kfolds','n_folds': 5,
'subject_id':subject_id})
There are several different ways to perform cross-validation. The standard
approach is to use k-folds, where the data is equally divided into k subsets
and each fold serves as both training and test.
Often we want to hold out the same subjects in each fold. This can be done by passing in a vector of unique subject IDs that
correspond to the images in the data frame.
In [ ]:
subject_id = data.X['SubjectID']
ridge_stats = data.predict(algorithm='ridge',
cv_dict={'type': 'kfolds','n_folds': 5,'subject_id':subject_id},
plot=False, **{'alpha':.1})
Sometimes we want to ensure that the training labels are balanced across folds. This can be done using the stratified k-folds method.
In [ ]:
ridge_stats = data.predict(algorithm='ridge',
cv_dict={'type': 'kfolds','n_folds': 5, 'stratified':data.Y},
plot=False, **{'alpha':.1})
Leave One Subject Out Cross-Validaiton (LOSO) is when k=n subjects.
This can be performed by passing in a vector indicating subject id's of
each image and using the loso flag.
In [ ]:
ridge_stats = data.predict(algorithm='ridge',
cv_dict={'type': 'loso','subject_id': subject_id},
plot=False, **{'alpha':.1})
There are also methods to estimate the shrinkage parameter for the penalized methods using nested crossvalidation with the ridgeCV and lassoCV algorithms.
In [ ]:
import numpy as np
ridgecv_stats = data.predict(algorithm='ridgeCV',
cv_dict={'type': 'kfolds','n_folds': 5, 'stratified':data.Y},
plot=False, **{'alphas':np.linspace(.1, 10, 5)})