.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/05_glm_second_level/plot_second_level_association_test.py" .. LINE NUMBERS ARE GIVEN BELOW. .. 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. .. GENERATED FROM PYTHON SOURCE LINES 14-17 .. code-block:: default # Author: Virgile Fritsch, Bertrand Thirion, 2014 -- 2018 # Jerome-Alexis Chevalier, 2019 .. GENERATED FROM PYTHON SOURCE LINES 18-19 At first, we need to load the Localizer contrasts. .. GENERATED FROM PYTHON SOURCE LINES 19-26 .. 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, legacy_format=False ) .. GENERATED FROM PYTHON SOURCE LINES 27-28 Let's print basic information on the dataset. .. GENERATED FROM PYTHON SOURCE LINES 28-31 .. 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 .. GENERATED FROM PYTHON SOURCE LINES 32-33 we also need to load the behavioral variable. .. GENERATED FROM PYTHON SOURCE LINES 33-36 .. code-block:: default tested_var = localizer_dataset.ext_vars['pseudo'] print(tested_var) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none 0 15.0 1 16.0 2 14.0 3 19.0 4 16.0 ... 89 12.0 90 16.0 91 13.0 92 25.0 93 21.0 Name: pseudo, Length: 94, dtype: float64 .. GENERATED FROM PYTHON SOURCE LINES 37-38 It is worth to do a auality check and remove subjects with missing values. .. GENERATED FROM PYTHON SOURCE LINES 38-48 .. code-block:: default import numpy as np mask_quality_check = np.where( np.logical_not(np.isnan(tested_var)) )[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].values.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 .. GENERATED FROM PYTHON SOURCE LINES 49-53 Estimate second level model --------------------------- We define the input maps and the design matrix for the second level model and fit it. .. GENERATED FROM PYTHON SOURCE LINES 53-58 .. code-block:: default import pandas as pd design_matrix = pd.DataFrame( np.hstack((tested_var, np.ones_like(tested_var))), columns=['fluency', 'intercept']) .. GENERATED FROM PYTHON SOURCE LINES 59-60 Fit of the second-level model .. GENERATED FROM PYTHON SOURCE LINES 60-64 .. 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) .. GENERATED FROM PYTHON SOURCE LINES 65-67 To estimate the contrast is very simple. We can just provide the column name of the design matrix. .. GENERATED FROM PYTHON SOURCE LINES 67-69 .. code-block:: default z_map = model.compute_contrast('fluency', output_type='z_score') .. GENERATED FROM PYTHON SOURCE LINES 70-71 We compute the fdr-corrected p = 0.05 threshold for these data .. GENERATED FROM PYTHON SOURCE LINES 71-74 .. code-block:: default from nilearn.glm import threshold_stats_img _, threshold = threshold_stats_img(z_map, alpha=.05, height_control='fdr') .. GENERATED FROM PYTHON SOURCE LINES 75-76 Let us plot the second level contrast at the computed thresholds .. GENERATED FROM PYTHON SOURCE LINES 76-84 .. 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-sg:: /auto_examples/05_glm_second_level/images/sphx_glr_plot_second_level_association_test_001.png :alt: plot second level association test :srcset: /auto_examples/05_glm_second_level/images/sphx_glr_plot_second_level_association_test_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 85-87 Computing the (corrected) p-values with parametric test to compare with non parametric test .. GENERATED FROM PYTHON SOURCE LINES 87-97 .. 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 .. GENERATED FROM PYTHON SOURCE LINES 98-99 Let us plot the (corrected) negative log p-values for the parametric test .. GENERATED FROM PYTHON SOURCE LINES 99-112 .. 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-sg:: /auto_examples/05_glm_second_level/images/sphx_glr_plot_second_level_association_test_002.png :alt: plot second level association test :srcset: /auto_examples/05_glm_second_level/images/sphx_glr_plot_second_level_association_test_002.png :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. .. GENERATED FROM PYTHON SOURCE LINES 113-114 Computing the (corrected) negative log p-values with permutation test .. GENERATED FROM PYTHON SOURCE LINES 114-123 .. 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) .. GENERATED FROM PYTHON SOURCE LINES 124-125 Let us plot the (corrected) negative log p-values .. GENERATED FROM PYTHON SOURCE LINES 125-136 .. 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-sg:: /auto_examples/05_glm_second_level/images/sphx_glr_plot_second_level_association_test_003.png :alt: plot second level association test :srcset: /auto_examples/05_glm_second_level/images/sphx_glr_plot_second_level_association_test_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 1 minutes 4.896 seconds) **Estimated memory usage:** 113 MB .. _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 `_