.. 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_04_glm_first_level_plot_predictions_residuals.py:
Predicted time series and residuals
===================================
Here we fit a First Level 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.
Import modules
--------------
.. code-block:: default
from nilearn.datasets import fetch_spm_auditory
from nilearn import image
from nilearn import masking
import pandas as pd
# load fMRI data
subject_data = fetch_spm_auditory()
fmri_img = image.concat_imgs(subject_data.func)
# 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=False)
fmri_img = image.smooth_img(fmri_img, 5.)
# load events
events = pd.read_table(subject_data['events'])
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
/home/nicolas/GitRepos/nilearn-fork/nilearn/image/image.py:1106: FutureWarning: The parameter "sessions" will be removed in 0.9.0 release of Nilearn. Please use the parameter "runs" instead.
data = signal.clean(
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`.
.. code-block:: default
from nilearn.glm.first_level import FirstLevelModel
fmri_glm = FirstLevelModel(t_r=7,
drift_model='cosine',
signal_scaling=False,
mask_img=mask,
minimize_memory=False)
fmri_glm = fmri_glm.fit(fmri_img, events)
Calculate and plot contrast
---------------------------
.. code-block:: default
from nilearn import plotting
z_map = fmri_glm.compute_contrast('active - rest')
plotting.plot_stat_map(z_map, bg_img=mean_img, threshold=3.1)
.. image:: /auto_examples/04_glm_first_level/images/sphx_glr_plot_predictions_residuals_001.png
:alt: plot predictions residuals
:class: sphx-glr-single-img
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
Extract the largest clusters
----------------------------
.. code-block:: default
from nilearn.reporting import get_clusters_table
from nilearn import input_data
table = get_clusters_table(z_map, stat_threshold=3.1,
cluster_threshold=20).set_index('Cluster ID', drop=True)
table.head()
# get the 6 largest clusters' max x, y, and z coordinates
coords = table.loc[range(1, 7), ['X', 'Y', 'Z']].values
# extract time series from each coordinate
masker = input_data.NiftiSpheresMasker(coords)
real_timeseries = masker.fit_transform(fmri_img)
predicted_timeseries = masker.fit_transform(fmri_glm.predicted[0])
Plot predicted and actual time series for 6 most significant clusters
---------------------------------------------------------------------
.. code-block:: default
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(0, 6):
# plotting time series
axs1[0, i].set_title('Cluster peak {}\n'.format(coords[i]))
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 = plotting.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)
.. image:: /auto_examples/04_glm_first_level/images/sphx_glr_plot_predictions_residuals_002.png
:alt: Cluster peak [-60. -6. 42.] , Cluster peak [60. 0. 36.] , Cluster peak [30. -9. 12.] , Cluster peak [-27. -3. 15.] , Cluster peak [57. 21. 75.] , Cluster peak [39. 33. 51.]
:class: sphx-glr-single-img
Get residuals
-------------
.. code-block:: default
resid = masker.fit_transform(fmri_glm.residuals[0])
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
/home/nicolas/GitRepos/nilearn-fork/nilearn/glm/regression.py:42: FutureWarning: 'resid' has been deprecated in version 0.7.0 and will be removed in version 0.9.0. Please use 'residuals' instead.
warnings.warn(category=FutureWarning,
Plot distribution of residuals
------------------------------
Note that residuals are not really distributed normally.
.. code-block:: default
fig2, axs2 = plt.subplots(2, 3)
axs2 = axs2.flatten()
for i in range(0, 6):
axs2[i].set_title('Cluster peak {}\n'.format(coords[i]))
axs2[i].hist(resid[:, i], color=colors[i])
print('Mean residuals: {}'.format(resid[:, i].mean()))
fig2.set_size_inches(12, 7)
fig2.tight_layout()
.. image:: /auto_examples/04_glm_first_level/images/sphx_glr_plot_predictions_residuals_003.png
:alt: Cluster peak [-60. -6. 42.] , Cluster peak [60. 0. 36.] , Cluster peak [30. -9. 12.] , Cluster peak [-27. -3. 15.] , Cluster peak [57. 21. 75.] , Cluster peak [39. 33. 51.]
:class: sphx-glr-single-img
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
Mean residuals: 8.141635513917814e-16
Mean residuals: -0.05424369099334695
Mean residuals: 0.0023122829555595934
Mean residuals: 0.011694403952955032
Mean residuals: 0.00237630017859917
Mean residuals: -0.0046459978771531225
Plot R-squared
--------------
Because we stored the residuals, we can plot the R-squared: the proportion of
explained variance of the 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 GLM. To isolate the influence of the experiment,
we can use an F-test as shown in the next section.
.. code-block:: default
plotting.plot_stat_map(fmri_glm.r_square[0], bg_img=mean_img, threshold=.1,
display_mode='z', cut_coords=7)
.. image:: /auto_examples/04_glm_first_level/images/sphx_glr_plot_predictions_residuals_004.png
:alt: plot predictions residuals
:class: sphx-glr-single-img
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
Calculate and Plot F-test
-------------------------
The F-test tells you how well the 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 GLM fits the data, including active, rest and
all the other columns in the design matrix such as drift and motion.
.. code-block:: default
import numpy as np
design_matrix = fmri_glm.design_matrices_[0]
# contrast with a one for "active" and zero everywhere else
active = np.array([1 if c == 'active' else 0 for c in design_matrix.columns])
# contrast with a one for "rest" and zero everywhere else
rest = np.array([1 if c == 'rest' else 0 for c in design_matrix.columns])
effects_of_interest = np.vstack((active, rest))
# f-test for rest and activity
z_map_ftest = fmri_glm.compute_contrast(
effects_of_interest,
stat_type='F',
output_type='z_score')
plotting.plot_stat_map(z_map_ftest, bg_img=mean_img, threshold=3.1,
display_mode='z', cut_coords=7)
.. image:: /auto_examples/04_glm_first_level/images/sphx_glr_plot_predictions_residuals_005.png
:alt: plot predictions residuals
:class: sphx-glr-single-img
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
.. rst-class:: sphx-glr-timing
**Total running time of the script:** ( 0 minutes 30.552 seconds)
.. _sphx_glr_download_auto_examples_04_glm_first_level_plot_predictions_residuals.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/04_glm_first_level/plot_predictions_residuals.ipynb
:alt: Launch binder
:width: 150 px
.. container:: sphx-glr-download sphx-glr-download-python
:download:`Download Python source code: plot_predictions_residuals.py `
.. container:: sphx-glr-download sphx-glr-download-jupyter
:download:`Download Jupyter notebook: plot_predictions_residuals.ipynb `
.. only:: html
.. rst-class:: sphx-glr-signature
`Gallery generated by Sphinx-Gallery `_