.. note:: :class: sphx-glr-download-link-note Click :ref:`here ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_03_second_level_models_plot_second_level_one_sample_test.py: Second-level fMRI model: one sample test ======================================== Full step-by-step example of fitting a GLM to perform a second-level analysis (one-sample test) and visualizing the results. More specifically: 1. A sequence of subject fMRI button press contrasts is downloaded. 2. A mask of the useful brain volume is computed. 3. A one-sample t-test is applied to the brain maps. We focus on a given contrast of the localizer dataset: the motor response to left versus right button press. Both at the individual and group level, this is expected to elicit activity in the motor cortex (positive in the right hemisphere, negative in the left hemisphere). Fetch dataset -------------- We download a list of left vs right button press contrasts from a localizer dataset. Note that we fetch individual t-maps that represent the Bold activity estimate divided by the uncertainty about this estimate. .. code-block:: default from nilearn.datasets import fetch_localizer_contrasts n_subjects = 16 data = fetch_localizer_contrasts(["left vs right button press"], n_subjects, get_tmaps=True) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none /home/kshitij/miniconda3/envs/nistats-py36-latest/lib/python3.6/site-packages/numpy/lib/npyio.py:2372: VisibleDeprecationWarning: Reading unicode strings without specifying the encoding argument is deprecated. Set the encoding, use None for the system default. output = genfromtxt(fname, **kwargs) Display subject t_maps ---------------------- We plot a grid with all the subjects t-maps thresholded at t = 2 for simple visualization purposes. The button press effect is visible among all subjects. .. code-block:: default from nilearn import plotting import matplotlib.pyplot as plt subjects = [subject_data[0] for subject_data in data['ext_vars']] fig, axes = plt.subplots(nrows=4, ncols=4) for cidx, tmap in enumerate(data['tmaps']): plotting.plot_glass_brain(tmap, colorbar=False, threshold=2.0, title=subjects[cidx], axes=axes[int(cidx / 4), int(cidx % 4)], plot_abs=False, display_mode='z') fig.suptitle('subjects t_map left-right button press') plt.show() .. image:: /auto_examples/03_second_level_models/images/sphx_glr_plot_second_level_one_sample_test_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out Out: .. code-block:: none /home/kshitij/workspace/nistats-org/nistats-repo/nistats-kchawla-pi/examples/03_second_level_models/plot_second_level_one_sample_test.py:49: UserWarning: Matplotlib is currently using agg, which is a non-GUI backend, so cannot show the figure. plt.show() 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 second_level_input = data['cmaps'] design_matrix = pd.DataFrame([1] * len(second_level_input), columns=['intercept']) Model specification and fit. .. code-block:: default from nistats.second_level_model import SecondLevelModel second_level_model = SecondLevelModel(smoothing_fwhm=8.0) second_level_model = second_level_model.fit(second_level_input, design_matrix=design_matrix) To estimate the contrast is very simple. We can just provide the column name of the design matrix. .. code-block:: default z_map = second_level_model.compute_contrast(output_type='z_score') We threshold the second level contrast at uncorrected p < 0.001 and plot it. .. code-block:: default from scipy.stats import norm p_val = 0.001 p001_unc = norm.isf(p_val) display = plotting.plot_glass_brain( z_map, threshold=p001_unc, colorbar=True, display_mode='z', plot_abs=False, title='group left-right button press (unc p<0.001)') plotting.show() .. image:: /auto_examples/03_second_level_models/images/sphx_glr_plot_second_level_one_sample_test_002.png :class: sphx-glr-single-img As expected, we find the motor cortex. Next, we compute the (corrected) p-values with a parametric test to compare them with the results from a nonparametric test. .. code-block:: default import numpy as np from nilearn.image import math_img from nilearn.input_data import NiftiMasker from nistats.utils import get_data p_val = second_level_model.compute_contrast(output_type='p_value') n_voxels = np.sum(get_data(second_level_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 = [0] # 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 discovery at all). # This threshold is much more conservative than the previous one. threshold = 1 title = ('Group left-right button press: \n' 'parametric test (FWER < 10%)') display = plotting.plot_glass_brain( neg_log_pval, colorbar=True, display_mode='z', plot_abs=False, vmax=3, cut_coords=cut_coords, threshold=threshold, title=title) plotting.show() .. image:: /auto_examples/03_second_level_models/images/sphx_glr_plot_second_level_one_sample_test_003.png :class: sphx-glr-single-img Now, we compute the (corrected) p-values with a permutation test. .. code-block:: default from nistats.second_level_model import non_parametric_inference neg_log_pvals_permuted_ols_unmasked = \ non_parametric_inference(second_level_input, design_matrix=design_matrix, model_intercept=True, n_perm=1000, two_sided_test=False, smoothing_fwhm=8.0, n_jobs=1) Let us plot the (corrected) negative log p-values for the nonparametric test. .. code-block:: default title = ('Group left-right button press: \n' 'permutation test (FWER < 10%)') display = plotting.plot_glass_brain( neg_log_pvals_permuted_ols_unmasked, colorbar=True, vmax=3, display_mode='z', plot_abs=False, cut_coords=cut_coords, threshold=threshold, title=title) plotting.show() # The neg-log p-values obtained with nonparametric testing are capped at 3 # since the number of permutations is 1e3. # The nonparametric test yields many more discoveries # and is more powerful than the usual parametric procedure. .. image:: /auto_examples/03_second_level_models/images/sphx_glr_plot_second_level_one_sample_test_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 16.101 seconds) .. _sphx_glr_download_auto_examples_03_second_level_models_plot_second_level_one_sample_test.py: .. only :: html .. container:: sphx-glr-footer :class: sphx-glr-footer-example .. container:: sphx-glr-download :download:`Download Python source code: plot_second_level_one_sample_test.py ` .. container:: sphx-glr-download :download:`Download Jupyter notebook: plot_second_level_one_sample_test.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_