Giving credit

Table Of Contents

Previous topic

2. Decoding and MVPA: predicting from brain images

Next topic

2.2. Choosing the right predictive model

2.1. A decoding tutorial

This page is a decoding tutorial articulated on the analysis of the Haxby 2001 dataset. It shows how to predict from brain activity images the stimuli that the subject is viewing.

2.1.1. Data loading and preparation

2.1.1.1. The Haxby 2001 experiment

Subjects are presented visual stimuli from different categories. We are going to predict which category the subject is seeing from the fMRI activity recorded in masks of the ventral stream. Significant prediction shows that the signal in the region contains information on the corresponding category.

../_images/sphx_glr_plot_haxby_stimuli_0041.png

Face stimuli

../_images/sphx_glr_plot_haxby_stimuli_0021.png

Cat stimuli

../_images/sphx_glr_plot_haxby_masks_0011.png

Masks

../_images/sphx_glr_plot_haxby_full_analysis_0011.png

Decoding scores per mask

2.1.1.2. Loading the data into Python

Full code example

A full example, with explanation, can be found on A introduction tutorial to fMRI decoding

  • Starting an environment: Launch IPython via “ipython –matplotlib” in a terminal, or use the Jupyter notebook.

  • Fetching the data: In the tutorial, we load the data using nilearn data downloading function, nilearn.datasets.fetch_haxby (if the data is not present on the disk, it can take a while to download about 310 Mo of data from the Internet):

    It return an object with several entries that contain paths to the files downloaded on the disk.

  • Loading the behavioral labels: Behavioral information is often stored in a text file such as a CSV, and must be load with numpy.recfromcsv or pandas

  • Preparing the fMRI data: we use the nilearn.input_data.NiftiMasker to apply the mask_vt mask to the 4D fMRI data, so that its shape becomes (n_samples, n_features) (see Masking data: from 4D Nifti images to 2D data arrays for a discussion on using masks).

  • Sample mask: Masking some of the data points may be useful to restrict to a specific pair of conditions (eg cats versus faces).

Note

Seemingly minor data preparation can matter a lot on the final score, for instance standardizing the data.

2.1.2. Performing the decoding analysis

2.1.2.1. The prediction engine

2.1.2.1.1. An estimator object

To perform decoding we construct an estimator, predicting a condition label y given a set X of images.

In the tutorial, we use a simple Support Vector Classification (or SVC) with a linear kernel. The corresponding class, sklearn.svm.SVC, needs to be imported from the scikit-learn.

Note that the documentation of the object details all parameters. In IPython, it can be displayed as follows:

In [10]: svc?
Type:             SVC
Base Class:       <class 'sklearn.svm.libsvm.SVC'>
String Form:
SVC(kernel=linear, C=1.0, probability=False, degree=3, coef0=0.0, tol=0.001,
cache_size=200, shrinking=True, gamma=0.0)
Namespace:        Interactive
Docstring:
    C-Support Vector Classification.
    Parameters
    ----------
    C : float, optional (default=1.0)
        penalty parameter C of the error term.
...

2.1.2.1.2. Applying it to data: fit (train) and predict (test)

In scikit-learn, the prediction objects have two important methods:

  • a fit function that “learns” the parameters of the model from the data. Thus, we need to give some training data to fit.
  • a predict function that “predicts” a target from new data. Here, we just have to give the new set of images (as the target should be unknown):

Warning

Do not predict on data used by the fit: the prediction that we obtain here is to good to be true (see next paragraph). Here we were just doing a sanity check.

2.1.2.2. Measuring prediction performance

2.1.2.2.1. Cross-validation

However, the last analysis is wrong, as we have learned and tested on the same set of data. We need to use a cross-validation to split the data into different sets, called “folds”, in a K-Fold strategy.

There is a specific function, sklearn.cross_validation.cross_val_score that computes for you the score for the different folds of cross-validation:

>>> from sklearn.cross_validation import cross_val_score
>>> cv_scores = cross_val_score(svc, fmri_masked, target, cv=5)  

cv=5 stipulates a 5-fold cross-validation.

You can speed up the computation by using n_jobs=-1, which will spread the computation equally across all processors (but will probably not work under Windows):

>>> cv_scores = cross_val_score(svc, fmri_masked, target, cv=5, n_jobs=-1, verbose=10) 

Prediction accuracy: We can take a look at the results of the cross_val_score function:

>>> print(cv_scores)  
[0.72727272727272729, 0.46511627906976744, 0.72093023255813948, 0.58139534883720934, 0.7441860465116279]

This is simply the prediction score for each fold, i.e. the fraction of correct predictions on the left-out data.

2.1.2.2.2. Choosing a good cross-validation strategy

There are many cross-validation strategies possible, including K-Fold or leave-one-out. When choosing a strategy, keep in mind that:

  • The test set should be as litte correlated as possible with the train set
  • The test set needs to have enough samples to enable a good measure of the prediction error (a rule of thumb is to use 10 to 20% of the data).

In these regards, leave one out is often one of the worst options.

Here, in the Haxby example, we are going to leave a session out, in order to have a test set independent from the train set. For this, we are going to use the session label, present in the behavioral data file, and sklearn.cross_validation.LeaveOneLabelOut.

Note

Full code for the above can be found on A introduction tutorial to fMRI decoding


Exercise

Compute the mean prediction accuracy using cv_scores.

Solution

>>> classification_accuracy = np.mean(cv_scores)  
>>> classification_accuracy  
0.76851...

We have a total prediction accuracy of 77% across the different sessions.

2.1.2.2.3. Choice of the prediction accuracy measure

The default metric used for measuring errors is the accuracy score, i.e. the number of total errors. It is not always a sensible metric, especially in the case of very imbalanced classes, as in such situations choosing the dominant class can achieve a low number of errors.

Other metrics, such as the f1-score, can be used:

>>> cv_scores = cross_val_score(svc, fmri_masked, target, cv=cv,  scoring='f1')  

2.1.2.2.4. Measuring the chance level

Dummy estimators: The simplest way to measure prediction performance at chance, is to use a “dummy” classifier, sklearn.dummy.DummyClassifier (purely random):

>>> from sklearn.dummy import DummyClassifier
>>> null_cv_scores = cross_val_score(DummyClassifier(), fmri_masked, target, cv=cv)  

Permutation testing: A more controlled way, but slower, is to do permutation testing on the labels, with sklearn.cross_validation.permutation_test_score:

>>> from sklearn.cross_validation import permutation_test_score
>>> null_cv_scores = permutation_test_score(svc, fmri_masked, target, cv=cv)  

Putting it all together

The ROI-based decoding example does a decoding analysis per mask, giving the f1-score of the prediction for each object.

It uses all the notions presented above, with for loop to iterate over masks and categories and Python dictionaries to store the scores.

../_images/sphx_glr_plot_haxby_masks_0011.png

Masks

../_images/sphx_glr_plot_haxby_full_analysis_0011.png

2.1.2.3. Visualizing the decoder’s weights

We can visualize the weights of the decoder:

  • we first inverse the masking operation, to retrieve a 3D brain volume of the SVC’s weights.
  • we then create a figure and plot as a background the first EPI image
  • finally we plot the SVC’s weights after masking the zero values
../_images/sphx_glr_plot_decoding_tutorial_0022.png

Note

Full code for the above can be found on A introduction tutorial to fMRI decoding

2.1.3. Decoding without a mask: Anova-SVM

2.1.3.1. Dimension reduction with feature selection

If we do not start from a mask of the relevant regions, there is a very large number of voxels and not all are useful for face vs cat prediction. We thus add a feature selection procedure. The idea is to select the k voxels most correlated to the task.

For this, we need to import the sklearn.feature_selection module and use sklearn.feature_selection.f_classif, a simple F-score based feature selection (a.k.a. Anova), that we will put before the SVC in a pipeline (sklearn.pipeline.Pipeline):

# Define the prediction function to be used.
# Here we use a Support Vector Classification, with a linear kernel
from sklearn.svm import SVC
svc = SVC(kernel='linear')

# Define the dimension reduction to be used.
# Here we use a classical univariate feature selection based on F-test,
# namely Anova. We set the number of features to be selected to 500
from sklearn.feature_selection import SelectKBest, f_classif
feature_selection = SelectKBest(f_classif, k=500)

# We have our classifier (SVC), our feature selection (SelectKBest), and now,
# we can plug them together in a *pipeline* that performs the two operations
# successively:
from sklearn.pipeline import Pipeline
anova_svc = Pipeline([('anova', feature_selection), ('svc', svc)])

#############################################################################
# Fit the decoder and predict

anova_svc.fit(X, y)
y_pred = anova_svc.predict(X)

#############################################################################

We can use our anova_svc object exactly as we were using our svc object previously.

2.1.3.2. Visualizing the results

To visualize the results, we need to:

  • first get the support vectors of the SVC and inverse the feature selection mechanism
  • then, as before, inverse the masking process to retrieve the weights and plot them.
../_images/sphx_glr_plot_haxby_anova_svm_0011.png
# Look at the SVC's discriminating weights
coef = svc.coef_
# reverse feature selection
coef = feature_selection.inverse_transform(coef)
# reverse masking
weight_img = nifti_masker.inverse_transform(coef)


# Create the figure
from nilearn import image
from nilearn.plotting import plot_stat_map, show

# Plot the mean image because we have no anatomic data
mean_img = image.mean_img(func_filename)

plot_stat_map(weight_img, mean_img, title='SVM weights')

# Saving the results as a Nifti file may also be important
weight_img.to_filename('haxby_face_vs_house.nii')

#############################################################################

Final script

The complete script to do an SVM-Anova analysis can be found as an example.

2.1.4. Going further with scikit-learn

We have seen a very simple analysis with scikit-learn, but it may be interesting to explore the wide variety of supervised learning algorithms in the scikit-learn.

2.1.4.1. Changing the prediction engine

We now see how one can easily change the prediction engine, if needed. We can try Fisher’s Linear Discriminant Analysis (LDA)

Import the module:

>>> from sklearn.lda import LDA

Construct the new estimator object and use it in a pipeline:

>>> from sklearn.pipeline import Pipeline
>>> lda = LDA()
>>> anova_lda = Pipeline([('anova', feature_selection), ('LDA', lda)])

and recompute the cross-validation score:

>>> cv_scores = cross_val_score(anova_lda, fmri_masked, target, cv=cv, verbose=1)  
>>> classification_accuracy = np.mean(cv_scores)  
>>> n_conditions = len(set(target))  # number of target classes
>>> print("Classification accuracy: %.4f / Chance Level: %.4f" % \
...    (classification_accuracy, 1. / n_conditions))  
Classification accuracy: 0.7846 / Chance level: 0.5000

2.1.4.2. Changing the feature selection

Let’s start by defining a linear SVM as a first classifier:

>>> clf = LinearSVC()

Let’s say that you want a more sophisticated feature selection, for example a Recursive Feature Elimination (RFE)

Import the module:

>>> from sklearn.feature_selection import RFE

Construct your new fancy selection:

>>> rfe = RFE(SVC(kernel='linear', C=1.), 50, step=0.25)

and create a new pipeline, composing the two classifiers rfe and clf:

>>> rfe_svc = Pipeline([('rfe', rfe), ('svc', clf)])

and recompute the cross-validation score:

>>> cv_scores = cross_val_score(rfe_svc, fmri_masked, target, cv=cv,
...     n_jobs=-1, verbose=1)  

But, be aware that this can take A WHILE...


See also

  • The scikit-learn documentation has very detailed explanations on a large variety of estimators and machine learning techniques. To become better at decoding, you need to study it.