.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/04_glm_first_level/plot_predictions_residuals.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` 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_04_glm_first_level_plot_predictions_residuals.py: Predicted time series and residuals =================================== Here we fit a First Level :term:`GLM` with the `minimize_memory`-argument set to `False`. By doing so, the `FirstLevelModel`-object stores the residuals, which we can then inspect. Also, the predicted time series can be extracted, which is useful to assess the quality of the model fit. .. include:: ../../../examples/masker_note.rst .. GENERATED FROM PYTHON SOURCE LINES 17-19 Import modules -------------- .. GENERATED FROM PYTHON SOURCE LINES 19-41 .. code-block:: Python import pandas as pd from nilearn import image, masking from nilearn.datasets import fetch_spm_auditory from nilearn.plotting import plot_stat_map, show # load fMRI data subject_data = fetch_spm_auditory() fmri_img = subject_data.func[0] # Make an average mean_img = image.mean_img(fmri_img) mask = masking.compute_epi_mask(mean_img) # Clean and smooth data fmri_img = image.clean_img(fmri_img, standardize=None) fmri_img = image.smooth_img(fmri_img, 5.0) # load events events = pd.read_csv(subject_data.events, sep="\t") .. rst-class:: sphx-glr-script-out .. code-block:: none [fetch_spm_auditory] Dataset found in /home/runner/nilearn_data/spm_auditory .. GENERATED FROM PYTHON SOURCE LINES 42-48 Fit model --------- Note that `minimize_memory` is set to `False` so that `FirstLevelModel` stores the residuals. `signal_scaling` is set to False, so we keep the same scaling as the original data in `fmri_img`. .. GENERATED FROM PYTHON SOURCE LINES 48-62 .. code-block:: Python from nilearn.glm.first_level import FirstLevelModel fmri_glm = FirstLevelModel( t_r=subject_data.t_r, drift_model="cosine", signal_scaling=False, mask_img=mask, minimize_memory=False, verbose=1, ) fmri_glm = fmri_glm.fit(fmri_img, events) .. rst-class:: sphx-glr-script-out .. code-block:: none [FirstLevelModel.fit] Loading mask from [FirstLevelModel.fit] Loading data from [FirstLevelModel.fit] Resampling mask [FirstLevelModel.fit] Finished fit [FirstLevelModel.fit] Computing run 1 out of 1 runs (go take a coffee, a big one). [FirstLevelModel.fit] Performing mask computation. [FirstLevelModel.fit] Loading data from [FirstLevelModel.fit] Extracting region signals [FirstLevelModel.fit] Cleaning extracted signals [FirstLevelModel.fit] Masking took 0 seconds. [FirstLevelModel.fit] Performing GLM computation. [FirstLevelModel.fit] GLM took 1 seconds. [FirstLevelModel.fit] Computation of 1 runs done in 1 seconds. .. GENERATED FROM PYTHON SOURCE LINES 63-65 Calculate and plot contrast --------------------------- .. GENERATED FROM PYTHON SOURCE LINES 65-78 .. code-block:: Python z_map = fmri_glm.compute_contrast("listening") threshold = 3.1 plot_stat_map( z_map, bg_img=mean_img, threshold=threshold, title=f"listening > rest (t-test; |Z|>{threshold})", ) show() .. image-sg:: /auto_examples/04_glm_first_level/images/sphx_glr_plot_predictions_residuals_001.png :alt: plot predictions residuals :srcset: /auto_examples/04_glm_first_level/images/sphx_glr_plot_predictions_residuals_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [FirstLevelModel.compute_contrast] Computing image from signals /home/runner/work/nilearn/nilearn/examples/04_glm_first_level/plot_predictions_residuals.py:75: UserWarning: You are using the 'agg' matplotlib backend that is non-interactive. No figure will be plotted when calling matplotlib.pyplot.show() or nilearn.plotting.show(). You can fix this by installing a different backend: for example via pip install PyQt6 .. GENERATED FROM PYTHON SOURCE LINES 79-85 Extract the largest clusters ---------------------------- We can extract the 6 largest clusters surviving our threshold. and get the x, y, and z coordinates of their peaks. We then extract the time series from a sphere around each coordinate. .. GENERATED FROM PYTHON SOURCE LINES 85-101 .. code-block:: Python from nilearn.maskers import NiftiSpheresMasker from nilearn.reporting import get_clusters_table table = get_clusters_table( z_map, stat_threshold=threshold, cluster_threshold=20 ) table.set_index("Cluster ID", drop=True) print(table) coords = table.loc[range(1, 7), ["X", "Y", "Z"]].to_numpy() print(coords) masker = NiftiSpheresMasker(coords, verbose=1) real_timeseries = masker.fit_transform(fmri_img) predicted_timeseries = masker.fit_transform(fmri_glm.predicted[0]) .. rst-class:: sphx-glr-script-out .. code-block:: none Cluster ID X Y Z Peak Stat Cluster Size (mm3) 0 1 -60.0 -6.0 42.0 11.543910 20196 1 1a -51.0 -12.0 39.0 10.497958 2 1b -39.0 -12.0 39.0 10.118978 3 1c -63.0 12.0 33.0 9.754835 4 2 60.0 0.0 36.0 11.088600 32157 .. ... ... ... ... ... ... 58 21 27.0 -72.0 27.0 4.422543 729 59 21a 24.0 -66.0 18.0 4.346647 60 21b 24.0 -78.0 21.0 3.428392 61 22 -30.0 -24.0 90.0 3.877083 648 62 22a -33.0 -18.0 72.0 3.791863 [63 rows x 6 columns] [[-51. -12. 39.] [-39. -12. 39.] [-63. 12. 33.] [ 60. 0. 36.] [ 66. 12. 27.] [ 60. -18. 30.]] [NiftiSpheresMasker.wrapped] Finished fit [NiftiSpheresMasker.wrapped] Loading data from [NiftiSpheresMasker.wrapped] Extracting region signals [NiftiSpheresMasker.wrapped] Cleaning extracted signals [FirstLevelModel.predicted] Computing image from signals [NiftiSpheresMasker.wrapped] Finished fit [NiftiSpheresMasker.wrapped] Loading data from [NiftiSpheresMasker.wrapped] Extracting region signals [NiftiSpheresMasker.wrapped] Cleaning extracted signals .. GENERATED FROM PYTHON SOURCE LINES 102-104 Let's have a look at the report to make sure the spheres are well placed. .. GENERATED FROM PYTHON SOURCE LINES 104-107 .. code-block:: Python report = masker.generate_report() report .. raw:: html

NiftiSpheresMasker Class for masking of Niimg-like objects using seeds. NiftiSpheresMasker is useful when data from given seeds should be extracted. Use case: summarize brain signals from seeds that were obtained from prior knowledge.

No image to display.

This report shows the regions defined by the spheres of the masker.

The masker has 6 spheres in total (displayed together on the first image).

They can be individually browsed using Previous and Next buttons.

Regions summary
seed number coordinates position radius size (in mm^3) size (in voxels) relative size (in %)
0 [-51.0, -12.0, 39.0] [48, 27, 35] 1.0 4.19 not implemented not implemented
1 [-39.0, -12.0, 39.0] [44, 27, 35] 1.0 4.19 not implemented not implemented
2 [-63.0, 12.0, 33.0] [52, 35, 33] 1.0 4.19 not implemented not implemented
3 [60.0, 0.0, 36.0] [11, 31, 34] 1.0 4.19 not implemented not implemented
4 [66.0, 12.0, 27.0] [9, 35, 31] 1.0 4.19 not implemented not implemented
5 [60.0, -18.0, 30.0] [11, 25, 32] 1.0 4.19 not implemented not implemented 
NiftiSpheresMasker(seeds=array([[-51., -12.,  39.],
           [-39., -12.,  39.],
           [-63.,  12.,  33.],
           [ 60.,   0.,  36.],
           [ 66.,  12.,  27.],
           [ 60., -18.,  30.]]),
                       verbose=1)
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.


.. GENERATED FROM PYTHON SOURCE LINES 108-110 Plot predicted and actual time series for 6 most significant clusters --------------------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 110-140 .. code-block:: Python import matplotlib.pyplot as plt # colors for each of the clusters colors = ["blue", "navy", "purple", "magenta", "olive", "teal"] # plot the time series and corresponding locations fig1, axs1 = plt.subplots(2, 6) for i in range(6): # plotting time series axs1[0, i].set_title(f"Cluster peak {coords[i]}\n") axs1[0, i].plot(real_timeseries[:, i], c=colors[i], lw=2) axs1[0, i].plot(predicted_timeseries[:, i], c="r", ls="--", lw=2) axs1[0, i].set_xlabel("Time") axs1[0, i].set_ylabel("Signal intensity", labelpad=0) # plotting image below the time series roi_img = plot_stat_map( z_map, cut_coords=[coords[i][2]], threshold=3.1, figure=fig1, axes=axs1[1, i], display_mode="z", colorbar=False, bg_img=mean_img, ) roi_img.add_markers([coords[i]], colors[i], 300) fig1.set_size_inches(24, 14) show() .. image-sg:: /auto_examples/04_glm_first_level/images/sphx_glr_plot_predictions_residuals_002.png :alt: Cluster peak [-51. -12. 39.] , Cluster peak [-39. -12. 39.] , Cluster peak [-63. 12. 33.] , Cluster peak [60. 0. 36.] , Cluster peak [66. 12. 27.] , Cluster peak [ 60. -18. 30.] :srcset: /auto_examples/04_glm_first_level/images/sphx_glr_plot_predictions_residuals_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none /home/runner/work/nilearn/nilearn/examples/04_glm_first_level/plot_predictions_residuals.py:137: UserWarning: You are using the 'agg' matplotlib backend that is non-interactive. No figure will be plotted when calling matplotlib.pyplot.show() or nilearn.plotting.show(). You can fix this by installing a different backend: for example via pip install PyQt6 .. GENERATED FROM PYTHON SOURCE LINES 141-143 Get residuals ------------- .. GENERATED FROM PYTHON SOURCE LINES 143-147 .. code-block:: Python resid = masker.fit_transform(fmri_glm.residuals[0]) .. rst-class:: sphx-glr-script-out .. code-block:: none [FirstLevelModel.residuals] Computing image from signals [NiftiSpheresMasker.wrapped] Finished fit [NiftiSpheresMasker.wrapped] Loading data from [NiftiSpheresMasker.wrapped] Extracting region signals [NiftiSpheresMasker.wrapped] Cleaning extracted signals .. GENERATED FROM PYTHON SOURCE LINES 148-151 Plot distribution of residuals ------------------------------ Note that residuals are not really distributed normally. .. GENERATED FROM PYTHON SOURCE LINES 151-164 .. code-block:: Python fig2, axs2 = plt.subplots(2, 3, constrained_layout=True) axs2 = axs2.flatten() for i in range(6): axs2[i].set_title(f"Cluster peak {coords[i]}\n") axs2[i].hist(resid[:, i], color=colors[i]) print(f"Mean residuals: {resid[:, i].mean()}") fig2.set_size_inches(12, 7) show() .. image-sg:: /auto_examples/04_glm_first_level/images/sphx_glr_plot_predictions_residuals_003.png :alt: Cluster peak [-51. -12. 39.] , Cluster peak [-39. -12. 39.] , Cluster peak [-63. 12. 33.] , Cluster peak [60. 0. 36.] , Cluster peak [66. 12. 27.] , Cluster peak [ 60. -18. 30.] :srcset: /auto_examples/04_glm_first_level/images/sphx_glr_plot_predictions_residuals_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Mean residuals: 0.018869405610299532 Mean residuals: 0.045756026122034346 Mean residuals: -1.5670004976137923e-14 Mean residuals: -0.13042393047259077 Mean residuals: -0.01499899481931678 Mean residuals: 0.049556260499625526 /home/runner/work/nilearn/nilearn/examples/04_glm_first_level/plot_predictions_residuals.py:161: UserWarning: You are using the 'agg' matplotlib backend that is non-interactive. No figure will be plotted when calling matplotlib.pyplot.show() or nilearn.plotting.show(). You can fix this by installing a different backend: for example via pip install PyQt6 .. GENERATED FROM PYTHON SOURCE LINES 165-181 Plot R-squared -------------- Because we stored the residuals, we can plot the R-squared: the proportion of explained variance of the :term:`GLM` as a whole. Note that the R-squared is markedly lower deep down the brain, where there is more physiological noise and we are further away from the receive coils. However, R-Squared should be interpreted with a grain of salt. The R-squared value will necessarily increase with the addition of more factors (such as rest, active, drift, motion) into the GLM. Additionally, we are looking at the overall fit of the model, so we are unable to say whether a voxel/region has a large R-squared value because the voxel/region is responsive to the experiment (such as active or rest) or because the voxel/region fits the noise factors (such as drift or motion) that could be present in the :term:`GLM`. To isolate the influence of the experiment, we can use an F-test as shown in the next section. .. GENERATED FROM PYTHON SOURCE LINES 181-195 .. code-block:: Python plot_stat_map( fmri_glm.r_square[0], bg_img=mean_img, threshold=0.1, display_mode="z", cut_coords=7, cmap="inferno", title="R-squared", vmin=0, symmetric_cbar=False, ) .. image-sg:: /auto_examples/04_glm_first_level/images/sphx_glr_plot_predictions_residuals_004.png :alt: plot predictions residuals :srcset: /auto_examples/04_glm_first_level/images/sphx_glr_plot_predictions_residuals_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [FirstLevelModel.r_square] Computing image from signals .. GENERATED FROM PYTHON SOURCE LINES 196-204 Calculate and Plot F-test ------------------------- The F-test tells you how well the :term:`GLM` fits effects of interest such as the active and rest conditions together. This is different from R-squared, which tells you how well the overall :term:`GLM` fits the data, including active, rest and all the other columns in the design matrix such as drift and motion. .. GENERATED FROM PYTHON SOURCE LINES 204-224 .. code-block:: Python # f-test for 'listening' z_map_ftest = fmri_glm.compute_contrast( "listening", stat_type="F", output_type="z_score" ) plot_stat_map( z_map_ftest, bg_img=mean_img, threshold=threshold, display_mode="z", cut_coords=7, cmap="inferno", title=f"listening > rest (F-test; Z>{threshold})", symmetric_cbar=False, vmin=0, ) show() .. image-sg:: /auto_examples/04_glm_first_level/images/sphx_glr_plot_predictions_residuals_005.png :alt: plot predictions residuals :srcset: /auto_examples/04_glm_first_level/images/sphx_glr_plot_predictions_residuals_005.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [FirstLevelModel.compute_contrast] Computing image from signals /home/runner/work/nilearn/nilearn/examples/04_glm_first_level/plot_predictions_residuals.py:222: UserWarning: You are using the 'agg' matplotlib backend that is non-interactive. No figure will be plotted when calling matplotlib.pyplot.show() or nilearn.plotting.show(). You can fix this by installing a different backend: for example via pip install PyQt6 .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 33.764 seconds) **Estimated memory usage:** 778 MB .. _sphx_glr_download_auto_examples_04_glm_first_level_plot_predictions_residuals.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: binder-badge .. image:: images/binder_badge_logo.svg :target: https://mybinder.org/v2/gh/nilearn/nilearn/0.13.1?urlpath=lab/tree/notebooks/auto_examples/04_glm_first_level/plot_predictions_residuals.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_predictions_residuals.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_predictions_residuals.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_predictions_residuals.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_