.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "_auto_examples/shortclips/05_plot_motion_energy_model.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 .. rst-class:: sphx-glr-example-title .. _sphx_glr__auto_examples_shortclips_05_plot_motion_energy_model.py: ============================================= Fit a ridge model with motion energy features ============================================= In this example, we model the fMRI responses with motion-energy features extracted from the movie stimulus. The model is a regularized linear regression model. This tutorial reproduces part of the analysis described in Nishimoto et al (2011) [1]_. See this publication for more details about the experiment, the motion-energy features, along with more results and more discussions. *Motion-energy features:* Motion-energy features result from filtering a video stimulus with spatio-temporal Gabor filters. A pyramid of filters is used to compute the motion-energy features at multiple spatial and temporal scales. Motion-energy features were introduced in [1]_. *Summary:* As in the previous example, we first concatenate the features with multiple delays, to account for the slow hemodynamic response. A linear regression model then weights each delayed feature with a different weight, to build a predictive model of BOLD activity. Again, the linear regression is regularized to improve robustness to correlated features and to improve generalization. The optimal regularization hyperparameter is selected independently on each voxel over a grid-search with cross-validation. Finally, the model generalization performance is evaluated on a held-out test set, comparing the model predictions with the ground-truth fMRI responses. .. GENERATED FROM PYTHON SOURCE LINES 30-32 Path of the data directory -------------------------- .. GENERATED FROM PYTHON SOURCE LINES 32-36 .. code-block:: Python from voxelwise_tutorials.io import get_data_home directory = get_data_home(dataset="shortclips") print(directory) .. rst-class:: sphx-glr-script-out .. code-block:: none /home/jlg/mvdoc/voxelwise_tutorials_data/shortclips .. GENERATED FROM PYTHON SOURCE LINES 37-41 .. code-block:: Python # modify to use another subject subject = "S01" .. GENERATED FROM PYTHON SOURCE LINES 42-46 Load the data ------------- We first load the fMRI responses. .. GENERATED FROM PYTHON SOURCE LINES 46-57 .. code-block:: Python import os import numpy as np from voxelwise_tutorials.io import load_hdf5_array file_name = os.path.join(directory, "responses", f"{subject}_responses.hdf") Y_train = load_hdf5_array(file_name, key="Y_train") Y_test = load_hdf5_array(file_name, key="Y_test") print("(n_samples_train, n_voxels) =", Y_train.shape) print("(n_repeats, n_samples_test, n_voxels) =", Y_test.shape) .. rst-class:: sphx-glr-script-out .. code-block:: none (n_samples_train, n_voxels) = (3600, 84038) (n_repeats, n_samples_test, n_voxels) = (10, 270, 84038) .. GENERATED FROM PYTHON SOURCE LINES 58-60 We average the test repeats, to remove the non-repeatable part of fMRI responses. .. GENERATED FROM PYTHON SOURCE LINES 60-64 .. code-block:: Python Y_test = Y_test.mean(0) print("(n_samples_test, n_voxels) =", Y_test.shape) .. rst-class:: sphx-glr-script-out .. code-block:: none (n_samples_test, n_voxels) = (270, 84038) .. GENERATED FROM PYTHON SOURCE LINES 65-66 We fill potential NaN (not-a-number) values with zeros. .. GENERATED FROM PYTHON SOURCE LINES 66-69 .. code-block:: Python Y_train = np.nan_to_num(Y_train) Y_test = np.nan_to_num(Y_test) .. GENERATED FROM PYTHON SOURCE LINES 70-71 Then we load the precomputed "motion-energy" features. .. GENERATED FROM PYTHON SOURCE LINES 71-80 .. code-block:: Python feature_space = "motion_energy" file_name = os.path.join(directory, "features", f"{feature_space}.hdf") X_train = load_hdf5_array(file_name, key="X_train") X_test = load_hdf5_array(file_name, key="X_test") print("(n_samples_train, n_features) =", X_train.shape) print("(n_samples_test, n_features) =", X_test.shape) .. rst-class:: sphx-glr-script-out .. code-block:: none (n_samples_train, n_features) = (3600, 6555) (n_samples_test, n_features) = (270, 6555) .. GENERATED FROM PYTHON SOURCE LINES 81-86 Define the cross-validation scheme ---------------------------------- We define the same leave-one-run-out cross-validation split as in the previous example. .. GENERATED FROM PYTHON SOURCE LINES 86-94 .. code-block:: Python from sklearn.model_selection import check_cv from voxelwise_tutorials.utils import generate_leave_one_run_out # indice of first sample of each run run_onsets = load_hdf5_array(file_name, key="run_onsets") print(run_onsets) .. rst-class:: sphx-glr-script-out .. code-block:: none [ 0 300 600 900 1200 1500 1800 2100 2400 2700 3000 3300] .. GENERATED FROM PYTHON SOURCE LINES 95-96 We define a cross-validation splitter, compatible with ``scikit-learn`` API. .. GENERATED FROM PYTHON SOURCE LINES 96-100 .. code-block:: Python n_samples_train = X_train.shape[0] cv = generate_leave_one_run_out(n_samples_train, run_onsets) cv = check_cv(cv) # copy the cross-validation splitter into a reusable list .. GENERATED FROM PYTHON SOURCE LINES 101-106 Define the model ---------------- We define the same model as in the previous example. See the previous example for more details about the model definition. .. GENERATED FROM PYTHON SOURCE LINES 106-128 .. code-block:: Python from sklearn.pipeline import make_pipeline from sklearn.preprocessing import StandardScaler from voxelwise_tutorials.delayer import Delayer from himalaya.kernel_ridge import KernelRidgeCV from himalaya.backend import set_backend backend = set_backend("torch_cuda", on_error="warn") X_train = X_train.astype("float32") X_test = X_test.astype("float32") alphas = np.logspace(1, 20, 20) pipeline = make_pipeline( StandardScaler(with_mean=True, with_std=False), Delayer(delays=[1, 2, 3, 4]), KernelRidgeCV( alphas=alphas, cv=cv, solver_params=dict(n_targets_batch=500, n_alphas_batch=5, n_targets_batch_refit=100)), ) .. GENERATED FROM PYTHON SOURCE LINES 129-133 .. code-block:: Python from sklearn import set_config set_config(display='diagram') # requires scikit-learn 0.23 pipeline .. raw:: html 
Pipeline(steps=[('standardscaler', StandardScaler(with_std=False)),
                    ('delayer', Delayer(delays=[1, 2, 3, 4])),
                    ('kernelridgecv',
                     KernelRidgeCV(alphas=array([1.e+01, 1.e+02, 1.e+03, 1.e+04, 1.e+05, 1.e+06, 1.e+07, 1.e+08,
           1.e+09, 1.e+10, 1.e+11, 1.e+12, 1.e+13, 1.e+14, 1.e+15, 1.e+16,
           1.e+17, 1.e+18, 1.e+19, 1.e+20]),
                                   cv=_CVIterableWrapper(cv=[(array([ 300,  301, ..., 3598, 3599]), a..., (array([   0,    1, ..., 3598, 3599]), array([3000, 3001, ..., 3298, 3299])), (array([   0,    1, ..., 3598, 3599]), array([1800, 1801, ..., 2098, 2099])), (array([   0,    1, ..., 3598, 3599]), array([1200, 1201, ..., 149... 2701, ..., 2998, 2999])), (array([   0,    1, ..., 3598, 3599]), array([600, 601, ..., 898, 899]))]),
                                   solver_params={'n_alphas_batch': 5,
                                                  'n_targets_batch': 500,
                                                  'n_targets_batch_refit': 100}))])
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 134-138 Fit the model ------------- We fit on the train set, and score on the test set. .. GENERATED FROM PYTHON SOURCE LINES 138-146 .. code-block:: Python pipeline.fit(X_train, Y_train) scores_motion_energy = pipeline.score(X_test, Y_test) scores_motion_energy = backend.to_numpy(scores_motion_energy) print("(n_voxels,) =", scores_motion_energy.shape) .. rst-class:: sphx-glr-script-out .. code-block:: none (n_voxels,) = (84038,) .. GENERATED FROM PYTHON SOURCE LINES 147-150 Plot the model performances --------------------------- The performances are computed using the :math:`R^2` scores. .. GENERATED FROM PYTHON SOURCE LINES 150-159 .. code-block:: Python import matplotlib.pyplot as plt from voxelwise_tutorials.viz import plot_flatmap_from_mapper mapper_file = os.path.join(directory, "mappers", f"{subject}_mappers.hdf") ax = plot_flatmap_from_mapper(scores_motion_energy, mapper_file, vmin=0, vmax=0.5) plt.show() .. image-sg:: /_auto_examples/shortclips/images/sphx_glr_05_plot_motion_energy_model_001.png :alt: 05 plot motion energy model :srcset: /_auto_examples/shortclips/images/sphx_glr_05_plot_motion_energy_model_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 160-163 The motion-energy features lead to large generalization scores in the early visual cortex (V1, V2, V3, ...). For more discussions about these results, we refer the reader to the original publication [1]_. .. GENERATED FROM PYTHON SOURCE LINES 165-171 Compare with the wordnet model ------------------------------ Interestingly, the motion-energy model performs well in different brain regions than the semantic "wordnet" model fitted in the previous example. To compare the two models, we first need to fit again the wordnet model. .. GENERATED FROM PYTHON SOURCE LINES 171-180 .. code-block:: Python feature_space = "wordnet" file_name = os.path.join(directory, "features", f"{feature_space}.hdf") X_train = load_hdf5_array(file_name, key="X_train") X_test = load_hdf5_array(file_name, key="X_test") X_train = X_train.astype("float32") X_test = X_test.astype("float32") .. GENERATED FROM PYTHON SOURCE LINES 181-183 We can create an unfitted copy of the pipeline with the ``clone`` function, or simply call fit again if we do not need to reuse the previous model. .. GENERATED FROM PYTHON SOURCE LINES 183-189 .. code-block:: Python if False: from sklearn.base import clone pipeline_wordnet = clone(pipeline) pipeline_wordnet .. GENERATED FROM PYTHON SOURCE LINES 190-198 .. code-block:: Python pipeline.fit(X_train, Y_train) scores_wordnet = pipeline.score(X_test, Y_test) scores_wordnet = backend.to_numpy(scores_wordnet) ax = plot_flatmap_from_mapper(scores_wordnet, mapper_file, vmin=0, vmax=0.5) plt.show() .. image-sg:: /_auto_examples/shortclips/images/sphx_glr_05_plot_motion_energy_model_002.png :alt: 05 plot motion energy model :srcset: /_auto_examples/shortclips/images/sphx_glr_05_plot_motion_energy_model_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 199-204 We can also plot the comparison of model prediction accuracies with a 2D histogram. All ~70k voxels are represented in this histogram, where the diagonal corresponds to identical prediction accuracy for both models. A distribution deviating from the diagonal means that one model has better predictive performance than the other. .. GENERATED FROM PYTHON SOURCE LINES 204-212 .. code-block:: Python from voxelwise_tutorials.viz import plot_hist2d ax = plot_hist2d(scores_wordnet, scores_motion_energy) ax.set(title='Generalization R2 scores', xlabel='semantic wordnet model', ylabel='motion energy model') plt.show() .. image-sg:: /_auto_examples/shortclips/images/sphx_glr_05_plot_motion_energy_model_003.png :alt: Generalization R2 scores :srcset: /_auto_examples/shortclips/images/sphx_glr_05_plot_motion_energy_model_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 213-216 Interestingly, the well predicted voxels are different in the two models. To further describe these differences, we can plot both performances on the same flatmap, using a 2D colormap. .. GENERATED FROM PYTHON SOURCE LINES 216-226 .. code-block:: Python from voxelwise_tutorials.viz import plot_2d_flatmap_from_mapper mapper_file = os.path.join(directory, "mappers", f"{subject}_mappers.hdf") ax = plot_2d_flatmap_from_mapper(scores_wordnet, scores_motion_energy, mapper_file, vmin=0, vmax=0.25, vmin2=0, vmax2=0.5, label_1="wordnet", label_2="motion energy") plt.show() .. image-sg:: /_auto_examples/shortclips/images/sphx_glr_05_plot_motion_energy_model_004.png :alt: 05 plot motion energy model :srcset: /_auto_examples/shortclips/images/sphx_glr_05_plot_motion_energy_model_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 227-249 The blue regions are well predicted by the motion-energy features, the orange regions are well predicted by the wordnet features, and the white regions are well predicted by both feature spaces. A large part of the visual semantic areas are not only well predicted by the wordnet features, but also by the motion-energy features, as indicated by the white color. Since these two features spaces encode quite different information, two interpretations are possible. In the first interpretation, the two feature spaces encode complementary information, and could be used jointly to further increase the generalization performance. In the second interpretation, both feature spaces encode the same information, because of spurious stimulus correlations. For example, imagine that the visual stimulus contained faces that appeared consistetly in the same portion of the visual field. In this case, position in the visual field would be perfectly correlated with the "face" semantic category. Thus, motion-energy features could predict responses in face-responsive areas without encoding any semantic information. To better disentangle the two feature spaces, we developed a joint model called `banded ridge regression` [2]_, which fits multiple feature spaces simultaneously with optimal regularization for each feature space. This model is described in the next example. .. GENERATED FROM PYTHON SOURCE LINES 251-261 References ---------- .. [1] Nishimoto, S., Vu, A. T., Naselaris, T., Benjamini, Y., Yu, B., & Gallant, J. L. (2011). Reconstructing visual experiences from brain activity evoked by natural movies. Current Biology, 21(19), 1641-1646. .. [2] Nunez-Elizalde, A. O., Huth, A. G., & Gallant, J. L. (2019). Voxelwise encoding models with non-spherical multivariate normal priors. Neuroimage, 197, 482-492. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 58.974 seconds) .. _sphx_glr_download__auto_examples_shortclips_05_plot_motion_energy_model.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 05_plot_motion_energy_model.ipynb <05_plot_motion_energy_model.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 05_plot_motion_energy_model.py <05_plot_motion_energy_model.py>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_