.. only:: html
.. note::
:class: sphx-glr-download-link-note
Click :ref:`here ` 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_05_glm_second_level_plot_second_level_two_sample_test.py:
Second-level fMRI model: two-sample test, unpaired and paired
==============================================================
Full step-by-step example of fitting a GLM to perform a second level analysis
in experimental data and visualizing the results
More specifically:
1. A sample of n=16 visual activity fMRIs are downloaded.
2. An unpaired, two-sample t-test is applied to the brain maps in order to
see the effect of the contrast difference across subjects.
3. A paired, two-sample t-test is applied to the brain maps in order to see
the effect of the contrast difference across subjects, considering subject intercepts
The contrast is between responses to retinotopically distinct
vertical versus horizontal checkerboards. At the individual level,
these stimuli are sometimes used to map the borders of primary visual areas.
At the group level, such a mapping is not possible. Yet, we may
observe some significant effects in these areas.
.. code-block:: default
import pandas as pd
from nilearn import plotting
from nilearn.datasets import fetch_localizer_contrasts
import matplotlib.pyplot as plt
Fetch dataset
--------------
We download a list of left vs right button press contrasts from a
localizer dataset.
.. code-block:: default
n_subjects = 16
sample_vertical = fetch_localizer_contrasts(
["vertical checkerboard"], n_subjects, get_tmaps=True)
sample_horizontal = fetch_localizer_contrasts(
["horizontal checkerboard"], n_subjects, get_tmaps=True)
# Implicitly, there is a one-to-one correspondence between the two samples:
# the first image of both samples comes from subject S1, the second from subject S2 etc.
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
/home/nicolas/anaconda3/envs/nilearn/lib/python3.8/site-packages/numpy/lib/npyio.py:2405: VisibleDeprecationWarning: Reading unicode strings without specifying the encoding argument is deprecated. Set the encoding, use None for the system default.
output = genfromtxt(fname, **kwargs)
Estimate second level models
------------------------------
We define the input maps and the design matrix for the second level model
and fit it.
.. code-block:: default
second_level_input = sample_vertical['cmaps'] + sample_horizontal['cmaps']
Next, we model the effect of conditions (sample 1 vs sample 2).
.. code-block:: default
import numpy as np
condition_effect = np.hstack(([1] * n_subjects, [- 1] * n_subjects))
The design matrix for the unpaired test doesn't need any more columns
For the paired test, we include an intercept for each subject.
.. code-block:: default
subject_effect = np.vstack((np.eye(n_subjects), np.eye(n_subjects)))
subjects = [f'S{i:02d}' for i in range(1, n_subjects + 1)]
We then assemble those into design matrices
.. code-block:: default
unpaired_design_matrix = pd.DataFrame(
condition_effect[:, np.newaxis],
columns=['vertical vs horizontal'])
paired_design_matrix = pd.DataFrame(
np.hstack((condition_effect[:, np.newaxis], subject_effect)),
columns=['vertical vs horizontal'] + subjects)
and plot the designs.
.. code-block:: default
from nilearn.plotting import plot_design_matrix
_, (ax_unpaired, ax_paired) = plt.subplots(1,2, gridspec_kw={'width_ratios': [1, 17]})
plot_design_matrix(unpaired_design_matrix, rescale=False, ax=ax_unpaired)
plot_design_matrix(paired_design_matrix, rescale=False, ax=ax_paired)
ax_unpaired.set_title('unpaired design', fontsize=12)
ax_paired.set_title('paired design', fontsize=12)
plt.tight_layout()
plotting.show()
.. image:: /auto_examples/05_glm_second_level/images/sphx_glr_plot_second_level_two_sample_test_001.png
:alt: unpaired design, paired design
:class: sphx-glr-single-img
We specify the analysis models and fit them.
.. code-block:: default
from nilearn.glm.second_level import SecondLevelModel
second_level_model_unpaired = SecondLevelModel().fit(
second_level_input, design_matrix=unpaired_design_matrix)
second_level_model_paired = SecondLevelModel().fit(
second_level_input, design_matrix=paired_design_matrix)
Estimating the contrast is simple. To do so, we provide the column
name of the design matrix. The argument 'output_type' is set to return all
available outputs so that we can compare differences in the effect size,
variance, and z-score.
.. code-block:: default
stat_maps_unpaired = second_level_model_unpaired.compute_contrast(
'vertical vs horizontal',
output_type='all')
stat_maps_paired = second_level_model_paired.compute_contrast(
'vertical vs horizontal',
output_type='all')
Plot the results
---------------------------
The two effect_size images are essentially identical
.. code-block:: default
(stat_maps_unpaired['effect_size'].get_fdata()
- stat_maps_paired['effect_size'].get_fdata()).max()
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
8.881784197001252e-16
But the variance in the unpaired image is larger.
.. code-block:: default
plotting.plot_glass_brain(
stat_maps_unpaired['effect_variance'], colorbar=True, vmin=0, vmax=6,
title='vertical vs horizontal effect variance, unpaired')
plotting.plot_glass_brain(
stat_maps_paired['effect_variance'], colorbar=True, vmin=0, vmax=6,
title='vertical vs horizontal effect variance, paired')
plotting.show()
.. rst-class:: sphx-glr-horizontal
*
.. image:: /auto_examples/05_glm_second_level/images/sphx_glr_plot_second_level_two_sample_test_002.png
:alt: plot second level two sample test
:class: sphx-glr-multi-img
*
.. image:: /auto_examples/05_glm_second_level/images/sphx_glr_plot_second_level_two_sample_test_003.png
:alt: plot second level two sample test
:class: sphx-glr-multi-img
Together, this makes the z_scores from the paired test larger.
We threshold the second level contrast and plot it.
.. code-block:: default
threshold = 3.1 # corresponds to p < .001, uncorrected
display = plotting.plot_glass_brain(
stat_maps_unpaired['z_score'], threshold=threshold, colorbar=True, plot_abs=False,
title='vertical vs horizontal (unc p<0.001)', vmin=0, vmax=6)
display = plotting.plot_glass_brain(
stat_maps_paired['z_score'], threshold=threshold, colorbar=True, plot_abs=False,
title='vertical vs horizontal (unc p<0.001)', vmin=0, vmax=6)
plotting.show()
.. rst-class:: sphx-glr-horizontal
*
.. image:: /auto_examples/05_glm_second_level/images/sphx_glr_plot_second_level_two_sample_test_004.png
:alt: plot second level two sample test
:class: sphx-glr-multi-img
*
.. image:: /auto_examples/05_glm_second_level/images/sphx_glr_plot_second_level_two_sample_test_005.png
:alt: plot second level two sample test
:class: sphx-glr-multi-img
Unsurprisingly, we see activity in the primary visual cortex, both positive
and negative.
.. rst-class:: sphx-glr-timing
**Total running time of the script:** ( 0 minutes 23.411 seconds)
.. _sphx_glr_download_auto_examples_05_glm_second_level_plot_second_level_two_sample_test.py:
.. only :: html
.. container:: sphx-glr-footer
:class: sphx-glr-footer-example
.. container:: binder-badge
.. image:: images/binder_badge_logo.svg
:target: https://mybinder.org/v2/gh/nilearn/nilearn.github.io/main?filepath=examples/auto_examples/05_glm_second_level/plot_second_level_two_sample_test.ipynb
:alt: Launch binder
:width: 150 px
.. container:: sphx-glr-download sphx-glr-download-python
:download:`Download Python source code: plot_second_level_two_sample_test.py `
.. container:: sphx-glr-download sphx-glr-download-jupyter
:download:`Download Jupyter notebook: plot_second_level_two_sample_test.ipynb `
.. only:: html
.. rst-class:: sphx-glr-signature
`Gallery generated by Sphinx-Gallery `_