.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples/07_advanced/plot_advanced_decoding_scikit.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_auto_examples_07_advanced_plot_advanced_decoding_scikit.py>`
        to download the full example code. or to run this example in your browser via Binder

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_auto_examples_07_advanced_plot_advanced_decoding_scikit.py:


Advanced decoding using scikit-learn
====================================

This tutorial opens the box of decoding pipelines, beyond the
functionalities provided by the
:class:`~nilearn.decoding.Decoder` object.
First, we reproduce basic functionalities of the
:class:`~nilearn.decoding.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 :ref:`documentation on decoding <decoding_intro>`
and in particular to the :ref:`advanced section <going_further>`.
As in many other examples, we decode the visual category of
stimuli in the :footcite:t:`Haxby2001` dataset,
focusing on distinguishing two categories:
"face" and "cat" images.

.. GENERATED FROM PYTHON SOURCE LINES 25-35

Retrieve and load the :term:`fMRI` data from the Haxby study
------------------------------------------------------------

Download the data
.................
The :func:`~nilearn.datasets.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``).

.. GENERATED FROM PYTHON SOURCE LINES 35-48

.. code-block:: Python


    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





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

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


.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>labels</th>
          <th>chunks</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>rest</td>
          <td>0</td>
        </tr>
        <tr>
          <th>1</th>
          <td>rest</td>
          <td>0</td>
        </tr>
        <tr>
          <th>2</th>
          <td>rest</td>
          <td>0</td>
        </tr>
        <tr>
          <th>3</th>
          <td>rest</td>
          <td>0</td>
        </tr>
        <tr>
          <th>4</th>
          <td>rest</td>
          <td>0</td>
        </tr>
        <tr>
          <th>...</th>
          <td>...</td>
          <td>...</td>
        </tr>
        <tr>
          <th>1447</th>
          <td>rest</td>
          <td>11</td>
        </tr>
        <tr>
          <th>1448</th>
          <td>rest</td>
          <td>11</td>
        </tr>
        <tr>
          <th>1449</th>
          <td>rest</td>
          <td>11</td>
        </tr>
        <tr>
          <th>1450</th>
          <td>rest</td>
          <td>11</td>
        </tr>
        <tr>
          <th>1451</th>
          <td>rest</td>
          <td>11</td>
        </tr>
      </tbody>
    </table>
    <p>1452 rows × 2 columns</p>
    </div>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 49-50

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

.. GENERATED FROM PYTHON SOURCE LINES 50-59

.. code-block:: Python

    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]








.. GENERATED FROM PYTHON SOURCE LINES 60-62

Performing decoding with scikit-learn
-------------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 64-71

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 :sklearn:`SVC <modules/svm.html>`).

.. GENERATED FROM PYTHON SOURCE LINES 71-76

.. code-block:: Python


    from sklearn.svm import SVC

    svc = SVC()








.. GENERATED FROM PYTHON SOURCE LINES 77-86

Masking the data
................
To use a scikit-learn estimator on brain images, you should first mask the
data using a :class:`~nilearn.maskers.NiftiMasker` to extract only the
voxels inside the mask of interest,
and transform 4D input :term:`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)`.

.. GENERATED FROM PYTHON SOURCE LINES 86-99

.. code-block:: Python

    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)





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    [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




.. GENERATED FROM PYTHON SOURCE LINES 100-105

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

.. GENERATED FROM PYTHON SOURCE LINES 105-111

.. code-block:: Python

    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}")





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    SVC accuracy: 0.823




.. GENERATED FROM PYTHON SOURCE LINES 112-125

Tuning cross-validation parameters
..................................
You can change many parameters of the cross_validation, such as:

* using a different
  :sklearn:`cross-validation scheme <modules/cross_validation.html>`.

* 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
  :sklearn:`SVC <modules/model_evaluation.html>`;
  for example, :func:`sklearn.metrics.roc_auc_score`.

.. GENERATED FROM PYTHON SOURCE LINES 125-139

.. code-block:: Python

    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}")





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    SVC accuracy (tuned parameters): 0.858




.. GENERATED FROM PYTHON SOURCE LINES 140-146

Measuring the chance level
--------------------------
:class:`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
:func:`sklearn.model_selection.permutation_test_score`.

.. GENERATED FROM PYTHON SOURCE LINES 148-150

Dummy estimator
...............

.. GENERATED FROM PYTHON SOURCE LINES 150-158

.. code-block:: Python

    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}")





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Dummy accuracy: 0.500




.. GENERATED FROM PYTHON SOURCE LINES 159-161

Permutation test
................

.. GENERATED FROM PYTHON SOURCE LINES 161-168

.. code-block:: Python

    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}")





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Permutation test score: 0.502




.. GENERATED FROM PYTHON SOURCE LINES 169-183

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 :mod:`sklearn.feature_selection` module and use
:func:`sklearn.feature_selection.f_classif`, a simple F-score
based feature selection (a.k.a.
`ANOVA <https://en.wikipedia.org/wiki/Analysis_of_variance#The_F-test>`_).

We can then chain both steps (feature selection and decoding)
into one composite estimator using
a :class:`~sklearn.pipeline.Pipeline` object.
Pipeline objects have several useful properties, as described in
the :sklearn:`scikit-learn documentation <modules/compose.html>`.

.. GENERATED FROM PYTHON SOURCE LINES 183-191

.. code-block:: Python

    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)])








.. GENERATED FROM PYTHON SOURCE LINES 192-200

We can now use our Pipeline ``anova_svc`` object exactly
as we were using our ``svc`` estimator previously.
Previously, we used :func:`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 :func:`sklearn.model_selection.cross_validate` function
with ``return_estimator = True`` to save the estimator.

.. GENERATED FROM PYTHON SOURCE LINES 200-212

.. code-block:: Python

    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}")





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    ANOVA+SVC test score: 0.801




.. GENERATED FROM PYTHON SOURCE LINES 213-217

Visualize the :term:`ANOVA` + SVC's discriminating weights
..........................................................
First, we retrieve the Pipeline object fitted on the first
cross-validation fold and its SVC coefficients.

.. GENERATED FROM PYTHON SOURCE LINES 217-225

.. code-block:: Python


    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"
    )





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

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




.. GENERATED FROM PYTHON SOURCE LINES 226-230

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.

.. GENERATED FROM PYTHON SOURCE LINES 230-237

.. code-block:: Python

    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"
    )





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    After inverting feature selection, we have 464 features back




.. GENERATED FROM PYTHON SOURCE LINES 238-241

Finally, we apply the ``inverse_transform`` function
of our :class:`~nilearn.maskers.NiftiMasker` object
to re-create a 4D Niimg that we can visualize.

.. GENERATED FROM PYTHON SOURCE LINES 241-248

.. code-block:: Python

    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()




.. image-sg:: /auto_examples/07_advanced/images/sphx_glr_plot_advanced_decoding_scikit_001.png
   :alt: plot advanced decoding scikit
   :srcset: /auto_examples/07_advanced/images/sphx_glr_plot_advanced_decoding_scikit_001.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    [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




.. GENERATED FROM PYTHON SOURCE LINES 249-255

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

.. GENERATED FROM PYTHON SOURCE LINES 257-264

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
:sklearn:`Linear Discriminant Analysis (LDA)
<auto_examples/decomposition/plot_pca_vs_lda.html>`.

.. GENERATED FROM PYTHON SOURCE LINES 264-286

.. code-block:: Python


    # 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}"
    )





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

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




.. GENERATED FROM PYTHON SOURCE LINES 287-295

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

.. GENERATED FROM PYTHON SOURCE LINES 295-314

.. code-block:: Python


    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....








.. GENERATED FROM PYTHON SOURCE LINES 315-319

References
----------

.. footbibliography::


.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 20.376 seconds)

**Estimated memory usage:**  1036 MB


.. _sphx_glr_download_auto_examples_07_advanced_plot_advanced_decoding_scikit.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: binder-badge

      .. image:: images/binder_badge_logo.svg
        :target: https://mybinder.org/v2/gh/nilearn/nilearn/0.14.1?urlpath=lab/tree/notebooks/auto_examples/07_advanced/plot_advanced_decoding_scikit.ipynb
        :alt: Launch binder
        :width: 150 px

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_advanced_decoding_scikit.ipynb <plot_advanced_decoding_scikit.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_advanced_decoding_scikit.py <plot_advanced_decoding_scikit.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: plot_advanced_decoding_scikit.zip <plot_advanced_decoding_scikit.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_