.. 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_thresholding.py: Statistical testing of a second-level analysis ============================================== Perform a one-sample t-test on a bunch of images (a.k.a. second-level analyis in fMRI) and threshold the resulting statistical map. This example is based on the so-called localizer dataset. It shows activation related to a mental computation task, as opposed to narrative sentence reading/listening. Prepare some images for a simple t test ---------------------------------------- This is a simple manually performed second level analysis. .. code-block:: default from nilearn import datasets n_samples = 20 localizer_dataset = datasets.fetch_localizer_calculation_task( n_subjects=n_samples) .. 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) Get the set of individual statstical maps (contrast estimates) .. code-block:: default cmap_filenames = localizer_dataset.cmaps Perform the second level analysis ---------------------------------- First, we define a design matrix for the model. As the model is trivial (one-sample test), the design matrix is just one column with ones. .. code-block:: default import pandas as pd design_matrix = pd.DataFrame([1] * n_samples, columns=['intercept']) Next, we specify and estimate the model. .. code-block:: default from nistats.second_level_model import SecondLevelModel second_level_model = SecondLevelModel().fit( cmap_filenames, design_matrix=design_matrix) Compute the only possible contrast: the one-sample test. Since there is only one possible contrast, we don't need to specify it in detail. .. code-block:: default z_map = second_level_model.compute_contrast(output_type='z_score') Threshold the resulting map: false positive rate < .001, cluster size > 10 voxels. .. code-block:: default from nistats.thresholding import map_threshold thresholded_map1, threshold1 = map_threshold( z_map, alpha=.001, height_control='fpr', cluster_threshold=10) Now use FDR <.05 (False Discovery Rate) and no cluster-level threshold. .. code-block:: default thresholded_map2, threshold2 = map_threshold( z_map, alpha=.05, height_control='fdr') print('The FDR=.05 threshold is %.3g' % threshold2) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none The FDR=.05 threshold is 2.37 Now use FWER <.05 (Family-Wise Error Rate) and no cluster-level threshold. As the data has not been intensively smoothed, we can use a simple Bonferroni correction. .. code-block:: default thresholded_map3, threshold3 = map_threshold( z_map, alpha=.05, height_control='bonferroni') print('The p<.05 Bonferroni-corrected threshold is %.3g' % threshold3) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none The p<.05 Bonferroni-corrected threshold is 4.88 Visualize the results --------------------- First, the unthresholded map. .. code-block:: default from nilearn import plotting display = plotting.plot_stat_map(z_map, title='Raw z map') .. image:: /auto_examples/03_second_level_models/images/sphx_glr_plot_thresholding_001.png :class: sphx-glr-single-img Second, the p<.001 uncorrected-thresholded map (with only clusters > 10 voxels). .. code-block:: default plotting.plot_stat_map( thresholded_map1, cut_coords=display.cut_coords, threshold=threshold1, title='Thresholded z map, fpr <.001, clusters > 10 voxels') .. image:: /auto_examples/03_second_level_models/images/sphx_glr_plot_thresholding_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out Out: .. code-block:: none Third, the fdr-thresholded map. .. code-block:: default plotting.plot_stat_map(thresholded_map2, cut_coords=display.cut_coords, title='Thresholded z map, expected fdr = .05', threshold=threshold2) .. image:: /auto_examples/03_second_level_models/images/sphx_glr_plot_thresholding_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out Out: .. code-block:: none Fourth, the Bonferroni-thresholded map. .. code-block:: default plotting.plot_stat_map(thresholded_map3, cut_coords=display.cut_coords, title='Thresholded z map, expected fwer < .05', threshold=threshold3) .. image:: /auto_examples/03_second_level_models/images/sphx_glr_plot_thresholding_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out Out: .. code-block:: none These different thresholds correspond to different statistical guarantees: in the FWER-corrected image there is only a probability smaller than .05 of observing any false positive voxel. In the FDR-corrected image, 5% of the voxels found are likely to be false positive. In the uncorrected image, one expects a few tens of false positive voxels. .. code-block:: default plotting.show() .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 3.315 seconds) .. _sphx_glr_download_auto_examples_03_second_level_models_plot_thresholding.py: .. only :: html .. container:: sphx-glr-footer :class: sphx-glr-footer-example .. container:: sphx-glr-download :download:`Download Python source code: plot_thresholding.py ` .. container:: sphx-glr-download :download:`Download Jupyter notebook: plot_thresholding.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_