.. only:: html

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

        Click :ref:`here <sphx_glr_download_auto_examples_07_advanced_plot_haxby_mass_univariate.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_haxby_mass_univariate.py:


Massively univariate analysis of face vs house recognition
==========================================================

A permuted Ordinary Least Squares algorithm is run at each voxel in
order to detemine whether or not it behaves differently under a "face
viewing" condition and a "house viewing" condition.
We consider the mean image per session and per condition.
Otherwise, the observations cannot be exchanged at random because
a time dependance exists between observations within a same session
(see [1] for more detailed explanations).

The example shows the small differences that exist between
Bonferroni-corrected p-values and family-wise corrected p-values obtained
from a permutation test combined with a max-type procedure [2].
Bonferroni correction is a bit conservative, as revealed by the presence of
a few false negative.

References
----------
[1] Winkler, A. M. et al. (2014).
    Permutation inference for the general linear model. Neuroimage.

[2] Anderson, M. J. & Robinson, J. (2001).
    Permutation tests for linear models.
    Australian & New Zealand Journal of Statistics, 43(1), 75-88.
    (http://avesbiodiv.mncn.csic.es/estadistica/permut2.pdf)


.. code-block:: default

    # Author: Virgile Fritsch, <virgile.fritsch@inria.fr>, Feb. 2014








Load Haxby dataset


.. code-block:: default

    from nilearn import datasets
    haxby_dataset = datasets.fetch_haxby(subjects=[2])

    # print basic information on the dataset
    print('Mask nifti image (3D) is located at: %s' % haxby_dataset.mask)
    print('Functional nifti image (4D) is located at: %s' % haxby_dataset.func[0])





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

 Out:

 .. code-block:: none

    Mask nifti image (3D) is located at: /home/varoquau/nilearn_data/haxby2001/mask.nii.gz
    Functional nifti image (4D) is located at: /home/varoquau/nilearn_data/haxby2001/subj2/bold.nii.gz




Mask data


.. code-block:: default

    mask_filename = haxby_dataset.mask
    from nilearn.input_data import NiftiMasker
    nifti_masker = NiftiMasker(
        smoothing_fwhm=8,
        mask_img=mask_filename,
        memory='nilearn_cache', memory_level=1)  # cache options
    func_filename = haxby_dataset.func[0]
    fmri_masked = nifti_masker.fit_transform(func_filename)








Restrict to faces and houses


.. code-block:: default

    import numpy as np
    import pandas as pd
    labels = pd.read_csv(haxby_dataset.session_target[0], sep=" ")
    conditions = labels['labels']
    categories = conditions.unique()
    conditions_encoded = np.zeros_like(conditions)
    for c, category in enumerate(categories):
        conditions_encoded[conditions == category] = c
    sessions = labels['chunks']
    condition_mask = conditions.isin(['face', 'house'])
    conditions_encoded = conditions_encoded[condition_mask]
    fmri_masked = fmri_masked[condition_mask]

    # We consider the mean image per session and per condition.
    # Otherwise, the observations cannot be exchanged at random because
    # a time dependence exists between observations within a same session.
    n_sessions = np.unique(sessions).size
    grouped_fmri_masked = np.empty((2 * n_sessions,  # two conditions per session
                                    fmri_masked.shape[1]))
    grouped_conditions_encoded = np.empty((2 * n_sessions, 1))

    for s in range(n_sessions):
        session_mask = sessions[condition_mask] == s
        session_house_mask = np.logical_and(session_mask,
                                            conditions[condition_mask] == 'house')
        session_face_mask = np.logical_and(session_mask,
                                           conditions[condition_mask] == 'face')
        grouped_fmri_masked[2 * s] = fmri_masked[session_house_mask].mean(0)
        grouped_fmri_masked[2 * s + 1] = fmri_masked[session_face_mask].mean(0)
        grouped_conditions_encoded[2 * s] = conditions_encoded[
            session_house_mask][0]
        grouped_conditions_encoded[2 * s + 1] = conditions_encoded[
            session_face_mask][0]








Perform massively univariate analysis with permuted OLS

We use a two-sided t-test to compute p-values, but we keep trace of the
effect sign to add it back at the end and thus observe the signed effect


.. code-block:: default

    from nilearn.mass_univariate import permuted_ols
    neg_log_pvals, t_scores_original_data, _ = permuted_ols(
        grouped_conditions_encoded, grouped_fmri_masked,
        # + intercept as a covariate by default
        n_perm=10000, two_sided_test=True,
        n_jobs=1)  # can be changed to use more CPUs
    signed_neg_log_pvals = neg_log_pvals * np.sign(t_scores_original_data)
    signed_neg_log_pvals_unmasked = nifti_masker.inverse_transform(
        signed_neg_log_pvals)








scikit-learn F-scores for comparison

F-test does not allow to observe the effect sign (pure two-sided test)


.. code-block:: default

    from sklearn.feature_selection import f_regression
    _, pvals_bonferroni = f_regression(
        grouped_fmri_masked,
        grouped_conditions_encoded)  # f_regression implicitly adds intercept
    pvals_bonferroni *= fmri_masked.shape[1]
    pvals_bonferroni[np.isnan(pvals_bonferroni)] = 1
    pvals_bonferroni[pvals_bonferroni > 1] = 1
    neg_log_pvals_bonferroni = -np.log10(pvals_bonferroni)
    neg_log_pvals_bonferroni_unmasked = nifti_masker.inverse_transform(
        neg_log_pvals_bonferroni)





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

 Out:

 .. code-block:: none

    /home/varoquau/dev/scikit-learn/sklearn/utils/validation.py:73: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples, ), for example using ravel().
      return f(**kwargs)




Visualization


.. code-block:: default

    import matplotlib.pyplot as plt
    from nilearn.plotting import plot_stat_map, show

    # Use the fmri mean image as a surrogate of anatomical data
    from nilearn import image
    from nilearn.image import get_data

    mean_fmri_img = image.mean_img(func_filename)

    threshold = -np.log10(0.1)  # 10% corrected

    vmax = min(signed_neg_log_pvals.max(),
               neg_log_pvals_bonferroni.max())

    # Plot thresholded p-values map corresponding to F-scores
    display = plot_stat_map(neg_log_pvals_bonferroni_unmasked, mean_fmri_img,
                            threshold=threshold, cmap=plt.cm.RdBu_r,
                            display_mode='z', cut_coords=[-1, ],
                            vmax=vmax)

    neg_log_pvals_bonferroni_data = get_data(neg_log_pvals_bonferroni_unmasked)
    n_detections = (neg_log_pvals_bonferroni_data > threshold).sum()
    title = ('Negative $\\log_{10}$ p-values'
             '\n(Parametric two-sided F-test'
             '\n+ Bonferroni correction)'
             '\n%d detections') % n_detections

    display.title(title, y=1.1)

    # Plot permutation p-values map
    display = plot_stat_map(signed_neg_log_pvals_unmasked, mean_fmri_img,
                            threshold=threshold, cmap=plt.cm.RdBu_r,
                            display_mode='z', cut_coords=[-1, ],
                            vmax=vmax)

    n_detections = (np.abs(signed_neg_log_pvals) > threshold).sum()
    title = ('Negative $\\log_{10}$ p-values'
             '\n(Non-parametric two-sided test'
             '\n+ max-type correction)'
             '\n%d detections') % n_detections

    display.title(title, y=1.1)

    show()



.. rst-class:: sphx-glr-horizontal


    *

      .. image:: /auto_examples/07_advanced/images/sphx_glr_plot_haxby_mass_univariate_001.png
          :alt: plot haxby mass univariate
          :class: sphx-glr-multi-img

    *

      .. image:: /auto_examples/07_advanced/images/sphx_glr_plot_haxby_mass_univariate_002.png
          :alt: plot haxby mass univariate
          :class: sphx-glr-multi-img


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

 Out:

 .. code-block:: none

    /home/varoquau/dev/nilearn/nilearn/plotting/displays.py:1608: MatplotlibDeprecationWarning: Adding an axes using the same arguments as a previous axes currently reuses the earlier instance.  In a future version, a new instance will always be created and returned.  Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.
      ax = fh.add_axes([fraction * index * (x1 - x0) + x0, y0,





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

   **Total running time of the script:** ( 1 minutes  21.975 seconds)


.. _sphx_glr_download_auto_examples_07_advanced_plot_haxby_mass_univariate.py:


.. only :: html

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


  .. container:: binder-badge

    .. image:: https://mybinder.org/badge_logo.svg
      :target: https://mybinder.org/v2/gh/nilearn/nilearn.github.io/master?filepath=examples/auto_examples/07_advanced/plot_haxby_mass_univariate.ipynb
      :width: 150 px


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

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



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

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


.. only:: html

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

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