.. 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_association_test.py:
Example of generic design in second-level models
================================================
This example shows the results obtained in a group analysis using a more
complex contrast than a one- or two-sample t test.
We use the [left button press (auditory cue)] task from the Localizer
dataset and seek association between the contrast values and a variate
that measures the speed of pseudo-word reading. No confounding variate
is included in the model.
.. code-block:: default
# Author: Virgile Fritsch, Bertrand Thirion, 2014 -- 2018
# Jerome-Alexis Chevalier, 2019
At first, we need to load the Localizer contrasts.
.. code-block:: default
from nilearn import datasets
n_samples = 94
localizer_dataset = datasets.fetch_localizer_contrasts(
['left button press (auditory cue)'], n_subjects=n_samples)
.. 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)
Let's print basic information on the dataset.
.. code-block:: default
print('First contrast nifti image (3D) is located at: %s' %
localizer_dataset.cmaps[0])
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
First contrast nifti image (3D) is located at: /home/nicolas/nilearn_data/brainomics_localizer/brainomics_data/S01/cmaps_LeftAuditoryClick.nii.gz
we also need to load the behavioral variable.
.. code-block:: default
tested_var = localizer_dataset.ext_vars['pseudo']
print(tested_var)
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
[b'15.0' b'16.0' b'14.0' b'19.0' b'16.0' b'18.0' b'22.0' b'19.0' b'17.0'
b'15.0' b'10.0' b'21.0' b'17.0' b'21.0' b'n/a' b'14.0' b'22.0' b'17.0'
b'23.0' b'15.0' b'15.0' b'18.0' b'17.0' b'18.0' b'20.0' b'27.0' b'18.0'
b'16.0' b'18.0' b'17.0' b'19.0' b'22.0' b'15.0' b'16.0' b'21.0' b'20.0'
b'12.0' b'n/a' b'19.0' b'19.0' b'16.0' b'22.0' b'23.0' b'14.0' b'24.0'
b'22.0' b'20.0' b'25.0' b'23.0' b'15.0' b'12.0' b'16.0' b'20.0' b'18.0'
b'14.0' b'14.0' b'18.0' b'20.0' b'19.0' b'14.0' b'27.0' b'n/a' b'13.0'
b'17.0' b'19.0' b'19.0' b'14.0' b'17.0' b'15.0' b'15.0' b'14.0' b'20.0'
b'16.0' b'15.0' b'15.0' b'15.0' b'19.0' b'17.0' b'14.0' b'15.0' b'n/a'
b'20.0' b'15.0' b'17.0' b'18.0' b'17.5' b'n/a' b'15.0' b'23.0' b'12.0'
b'16.0' b'13.0' b'25.0' b'21.0']
It is worth to do a auality check and remove subjects with missing values.
.. code-block:: default
import numpy as np
mask_quality_check = np.where(tested_var != b'n/a')[0]
n_samples = mask_quality_check.size
contrast_map_filenames = [localizer_dataset.cmaps[i]
for i in mask_quality_check]
tested_var = tested_var[mask_quality_check].astype(float).reshape((-1, 1))
print("Actual number of subjects after quality check: %d" % n_samples)
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
Actual number of subjects after quality check: 89
Estimate second level model
---------------------------
We define the input maps and the design matrix for the second level model
and fit it.
.. code-block:: default
import pandas as pd
design_matrix = pd.DataFrame(
np.hstack((tested_var, np.ones_like(tested_var))),
columns=['fluency', 'intercept'])
Fit of the second-level model
.. code-block:: default
from nilearn.glm.second_level import SecondLevelModel
model = SecondLevelModel(smoothing_fwhm=5.0)
model.fit(contrast_map_filenames, design_matrix=design_matrix)
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
SecondLevelModel(smoothing_fwhm=5.0)
To estimate the contrast is very simple. We can just provide the column
name of the design matrix.
.. code-block:: default
z_map = model.compute_contrast('fluency', output_type='z_score')
We compute the fdr-corrected p = 0.05 threshold for these data
.. code-block:: default
from nilearn.glm import threshold_stats_img
_, threshold = threshold_stats_img(z_map, alpha=.05, height_control='fdr')
Let us plot the second level contrast at the computed thresholds
.. code-block:: default
from nilearn import plotting
plotting.plot_stat_map(
z_map, threshold=threshold, colorbar=True,
title='Group-level association between motor activity \n'
'and reading fluency (fdr=0.05)')
plotting.show()
.. image:: /auto_examples/05_glm_second_level/images/sphx_glr_plot_second_level_association_test_001.png
:alt: plot second level association test
:class: sphx-glr-single-img
Computing the (corrected) p-values with parametric test to compare with
non parametric test
.. code-block:: default
from nilearn.image import math_img
from nilearn.image import get_data
p_val = model.compute_contrast('fluency', output_type='p_value')
n_voxels = np.sum(get_data(model.masker_.mask_img_))
# Correcting the p-values for multiple testing and taking negative logarithm
neg_log_pval = math_img("-np.log10(np.minimum(1, img * {}))"
.format(str(n_voxels)),
img=p_val)
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
:1: RuntimeWarning: divide by zero encountered in log10
Let us plot the (corrected) negative log p-values for the parametric test
.. code-block:: default
cut_coords = [38, -17, -3]
# Since we are plotting negative log p-values and using a threshold equal to 1,
# it corresponds to corrected p-values lower than 10%, meaning that there
# is less than 10% probability to make a single false discovery
# (90% chance that we make no false discoveries at all).
# This threshold is much more conservative than the previous one.
threshold = 1
title = ('Group-level association between motor activity and reading: \n'
'neg-log of parametric corrected p-values (FWER < 10%)')
plotting.plot_stat_map(neg_log_pval, colorbar=True, cut_coords=cut_coords,
threshold=threshold, title=title)
plotting.show()
.. image:: /auto_examples/05_glm_second_level/images/sphx_glr_plot_second_level_association_test_002.png
:alt: plot second level association test
:class: sphx-glr-single-img
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
/home/nicolas/GitRepos/nilearn-fork/nilearn/_utils/niimg.py:61: UserWarning: Non-finite values detected. These values will be replaced with zeros.
warn(
Computing the (corrected) negative log p-values with permutation test
.. code-block:: default
from nilearn.glm.second_level import non_parametric_inference
neg_log_pvals_permuted_ols_unmasked = \
non_parametric_inference(contrast_map_filenames,
design_matrix=design_matrix,
second_level_contrast='fluency',
model_intercept=True, n_perm=1000,
two_sided_test=False, mask=None,
smoothing_fwhm=5.0, n_jobs=1)
Let us plot the (corrected) negative log p-values
.. code-block:: default
title = ('Group-level association between motor activity and reading: \n'
'neg-log of non-parametric corrected p-values (FWER < 10%)')
plotting.plot_stat_map(neg_log_pvals_permuted_ols_unmasked, colorbar=True,
cut_coords=cut_coords, threshold=threshold,
title=title)
plotting.show()
# The neg-log p-values obtained with non parametric testing are capped at 3
# since the number of permutations is 1e3.
# The non parametric test yields a few more discoveries
# and is then more powerful than the usual parametric procedure.
.. image:: /auto_examples/05_glm_second_level/images/sphx_glr_plot_second_level_association_test_003.png
:alt: plot second level association test
:class: sphx-glr-single-img
.. rst-class:: sphx-glr-timing
**Total running time of the script:** ( 0 minutes 51.691 seconds)
.. _sphx_glr_download_auto_examples_05_glm_second_level_plot_second_level_association_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_association_test.ipynb
:alt: Launch binder
:width: 150 px
.. container:: sphx-glr-download sphx-glr-download-python
:download:`Download Python source code: plot_second_level_association_test.py `
.. container:: sphx-glr-download sphx-glr-download-jupyter
:download:`Download Jupyter notebook: plot_second_level_association_test.ipynb `
.. only:: html
.. rst-class:: sphx-glr-signature
`Gallery generated by Sphinx-Gallery `_