Advanced decoding using scikit-learn

This tutorial opens the box of decoding pipelines, beyond the functionalities provided by the Decoder object. First, we reproduce basic functionalities of the Decoder object via direct calls to the underlying scikit-learn functions. Next, we give pointers towards integrating other scikit-learn estimators directly.

If some concepts seem unclear, please refer to the documentation on decoding and in particular to the advanced section. As in many other examples, we decode the visual category of stimuli in the Haxby et al.[1] dataset, focusing on distinguishing two categories: “face” and “cat” images.

Retrieve and load the fMRI data from the Haxby study

Download the data

The fetch_haxby function will download the Haxby dataset object, whose attributes include the fMRI images as Niimg objects (func), a spatial mask (mask_vt), and a CSV with the visual category label for each image (session_target).

from nilearn import datasets

haxby_dataset = datasets.fetch_haxby()
mask_filename = haxby_dataset.mask_vt[0]
fmri_filename = haxby_dataset.func[0]

# Loading the behavioral labels
import pandas as pd

behavioral = pd.read_csv(haxby_dataset.session_target[0], delimiter=" ")
behavioral

[fetch_haxby] Dataset directory found: /home/runner/work/nilearn/nilearn/nilearn_data/haxby2001

	labels	chunks
0	rest	0
1	rest	0
2	rest	0
3	rest	0
4	rest	0
...	...	...
1447	rest	11
1448	rest	11
1449	rest	11
1450	rest	11
1451	rest	11

1452 rows × 2 columns

We keep only a images from the conditions of interest (“cat” and “face”).

from nilearn.image import index_img

conditions = behavioral["labels"]
condition_mask = conditions.isin(["face", "cat"])
fmri_niimgs = index_img(fmri_filename, condition_mask)
conditions = conditions[condition_mask]
conditions = conditions.to_numpy()
run_label = behavioral["chunks"][condition_mask]

Performing decoding with scikit-learn

Importing a classifier

We can import many predictive models from scikit-learn that can be used in a decoding pipelines. They all support a .fit() method. Let’s define a Support Vector Classifier (or SVC).

from sklearn.svm import SVC

svc = SVC()

Masking the data

To use a scikit-learn estimator on brain images, you should first mask the data using a NiftiMasker to extract only the voxels inside the mask of interest, and transform 4D input fMRI data to 2D arrays of shape (n_samples, n_features) that scikit-learn estimators can work on. In our case, this means extracting arrays of shape (n_timepoints, n_voxels).

from nilearn.maskers import NiftiMasker

masker = NiftiMasker(
    mask_img=mask_filename,
    runs=run_label,
    smoothing_fwhm=4,
    standardize="zscore_sample",
    memory="nilearn_cache",
    memory_level=1,
    verbose=1,
)
fmri_masked = masker.fit_transform(fmri_niimgs)

[NiftiMasker.wrapped] Loading data from <nibabel.nifti1.Nifti1Image object at 0x7f2d98604790>
[NiftiMasker.wrapped] Loading mask from .../mask4_vt.nii.gz
[NiftiMasker.wrapped] Resampling mask
________________________________________________________________________________
[Memory] Calling nilearn.image.resampling.resample_img...
resample_img(<nibabel.nifti1.Nifti1Image object at 0x7f2d98605f60>, target_affine=None, target_shape=None, copy=False, interpolation='nearest')
_____________________________________________________resample_img - 0.0s, 0.0min
[NiftiMasker.wrapped] Finished fit
________________________________________________________________________________
[Memory] Calling nilearn.maskers.nifti_masker.filter_and_mask...
filter_and_mask(<nibabel.nifti1.Nifti1Image object at 0x7f2d98604790>, <nibabel.nifti1.Nifti1Image object at 0x7f2d98605f60>, { 'clean_args': None,
  'clean_kwargs': {},
  'cmap': 'gray',
  'detrend': False,
  'dtype': None,
  'high_pass': None,
  'high_variance_confounds': False,
  'low_pass': None,
  'reports': True,
  'runs': 21       0
22       0
23       0
24       0
25       0
        ..
1427    11
1428    11
1429    11
1430    11
1431    11
Name: chunks, Length: 216, dtype: int64,
  'smoothing_fwhm': 4,
  'standardize': 'zscore_sample',
  'standardize_confounds': True,
  't_r': None,
  'target_affine': None,
  'target_shape': None}, memory_level=1, memory=Memory(location=nilearn_cache/joblib), verbose=1, confounds=None, sample_mask=None, copy=True, sklearn_output_config=None)
[NiftiMasker.wrapped] Loading data from <nibabel.nifti1.Nifti1Image object at 0x7f2d98604790>
[NiftiMasker.wrapped] Smoothing images
[NiftiMasker.wrapped] Extracting region signals
[NiftiMasker.wrapped] Cleaning extracted signals
__________________________________________________filter_and_mask - 1.2s, 0.0min

Cross-validation with scikit-learn

To train and test the model in a meaningful way we use cross-validation with the function sklearn.model_selection.cross_val_score that computes the score for each of the different cross-validation folds.

from sklearn.model_selection import cross_val_score

# Here `cv=5` stipulates a 5-fold cross-validation
cv_scores = cross_val_score(svc, fmri_masked, conditions, cv=5)
print(f"SVC accuracy: {cv_scores.mean():.3f}")

SVC accuracy: 0.823

Tuning cross-validation parameters

You can change many parameters of the cross_validation, such as:

• using a different cross-validation scheme.

• speeding up the computation by using n_jobs = -1, which will spread the computation equally across all processors.

• use a different scoring function, as a keyword or imported from SVC; for example, sklearn.metrics.roc_auc_score.

from sklearn.model_selection import LeaveOneGroupOut

cv = LeaveOneGroupOut()
cv_scores = cross_val_score(
    svc,
    fmri_masked,
    conditions,
    cv=cv,
    scoring="roc_auc",
    groups=run_label,
    n_jobs=2,
)
print(f"SVC accuracy (tuned parameters): {cv_scores.mean():.3f}")

SVC accuracy (tuned parameters): 0.858

Measuring the chance level

sklearn.dummy.DummyClassifier (purely random) estimators are the simplest way to measure prediction performance at chance. A more controlled way, but slower, is to do permutation testing on the labels, with sklearn.model_selection.permutation_test_score.

Dummy estimator

from sklearn.dummy import DummyClassifier

null_cv_scores = cross_val_score(
    DummyClassifier(), fmri_masked, conditions, cv=cv, groups=run_label
)

print(f"Dummy accuracy: {null_cv_scores.mean():.3f}")

Dummy accuracy: 0.500

Permutation test

from sklearn.model_selection import permutation_test_score

null_cv_scores = permutation_test_score(
    svc, fmri_masked, conditions, cv=cv, groups=run_label
)[1]
print(f"Permutation test score: {null_cv_scores.mean():.3f}")

Permutation test score: 0.502

Decoding without a mask: ANOVA-SVM in scikit-learn

We can also implement feature selection before decoding. To perform the feature selection, 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).

We can then chain both steps (feature selection and decoding) into one composite estimator using a Pipeline object. Pipeline objects have several useful properties, as described in the scikit-learn documentation.

from sklearn.feature_selection import SelectPercentile, f_classif
from sklearn.pipeline import Pipeline
from sklearn.svm import LinearSVC

feature_selection = SelectPercentile(f_classif, percentile=10)
linear_svc = LinearSVC(dual=True, random_state=0)
anova_svc = Pipeline([("anova", feature_selection), ("svc", linear_svc)])

We can now use our Pipeline anova_svc object exactly as we were using our svc estimator previously. Previously, we used sklearn.model_selection.cross_val_score to return the cross-validated decoding scores. However, we now want to investigate our model’s feature selection via its weights. We can use sklearn.model_selection.cross_validate function with return_estimator = True to save the estimator.

from sklearn.model_selection import cross_validate

fitted_pipeline = cross_validate(
    anova_svc,
    fmri_masked,
    conditions,
    cv=cv,
    groups=run_label,
    return_estimator=True,
)
print(f"ANOVA+SVC test score: {fitted_pipeline['test_score'].mean():.3f}")

ANOVA+SVC test score: 0.801

Visualize the ANOVA + SVC’s discriminating weights

First, we retrieve the Pipeline object fitted on the first cross-validation fold and its SVC coefficients.

first_pipeline = fitted_pipeline["estimator"][0]
svc_coef = first_pipeline.named_steps["svc"].coef_
print(
    "After feature selection, "
    f"the SVC is trained only on {svc_coef.shape[1]} features"
)

After feature selection, the SVC is trained only on 47 features

Next, we use the inverse_transform function to invert the feature selection step and put these coefficients in the right place in our (n_timepoints, n_voxels) 2D array.

full_coef = first_pipeline.named_steps["anova"].inverse_transform(svc_coef)

print(
    "After inverting feature selection, "
    f"we have {full_coef.shape[1]} features back"
)

After inverting feature selection, we have 464 features back

Finally, we apply the inverse_transform function of our NiftiMasker object to re-create a 4D Niimg that we can visualize.

from nilearn.plotting import plot_stat_map, show

weight_img = masker.inverse_transform(full_coef)
plot_stat_map(weight_img, title="ANOVA+SVC weights", draw_cross=False)

show()

[NiftiMasker.inverse_transform] Computing image from signals
________________________________________________________________________________
[Memory] Calling nilearn.masking.unmask...
unmask(array([[-0.341414, ...,  0.      ]]), <nibabel.nifti1.Nifti1Image object at 0x7f2d98605f60>)
___________________________________________________________unmask - 0.0s, 0.0min

Going further with scikit-learn

While the above analysis mirrored what occurs in the Decoder object, we can go still further with scikit-learn. Two examples are given below, but many more are possible.

Changing the prediction engine

To change the prediction engine, we just need to import it and use in our pipeline instead of the SVC. For example, we can try Fisher’s Linear Discriminant Analysis (LDA).

# Construct the new estimator object and use it in a new Pipeline
# after feature-selection with ANOVA, as before
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis

feature_selection = SelectPercentile(f_classif, percentile=10)
lda = LinearDiscriminantAnalysis()
anova_lda = Pipeline([("anova", feature_selection), ("LDA", lda)])

# Recompute the cross-validation score:
import numpy as np

cv_scores = cross_val_score(
    anova_lda, fmri_masked, conditions, cv=cv, groups=run_label
)
classification_accuracy = np.mean(cv_scores)
n_conditions = len(set(conditions))  # number of target classes
print(
    f"ANOVA + LDA classification accuracy: {classification_accuracy:.4f} "
    f"/ Chance Level: {1.0 / n_conditions:.4f}"
)

ANOVA + LDA classification accuracy: 0.8009 / Chance Level: 0.5000

Changing the feature selection

Let’s say that you want a more sophisticated feature selection; for example, a Recursive Feature Elimination (RFE) before a SVC. We can simply follow the same principle as we did in changing the prediction engine.

from sklearn.feature_selection import RFE

svc = SVC()
rfe = RFE(SVC(kernel="linear", C=1.0), n_features_to_select=50, step=0.25)

# Create a new pipeline, composing the two classifiers `rfe` and `svc`.

rfe_svc = Pipeline([("rfe", rfe), ("svc", svc)])

# Recompute the cross-validation score
# cv_scores = cross_val_score(rfe_svc,
#                             fmri_masked,
#                             target,
#                             cv=cv,
#                             n_jobs=2,
#                             verbose=1)
# But, be aware that this can take some time....

References

[1]

James V. Haxby, M. Ida Gobbini, Maura L. Furey, Alumit Ishai, Jennifer L. Schouten, and Pietro Pietrini. Distributed and overlapping representations of faces and objects in ventral temporal cortex. Science, 293(5539):2425–2430, 2001. doi:10.1126/science.1063736.

Estimated memory usage: 1036 MB