|
37 | 37 | ############################################################################### |
38 | 38 | # First, we load the fMRI responses on the test set, which contains 10 repeats. |
39 | 39 | file_name = os.path.join(directory, 'responses', f'{subject}_responses.hdf') |
40 | | -Y_test = load_hdf5_array(file_name, "Y_test") |
| 40 | +Y_test = load_hdf5_array(file_name, key="Y_test") |
41 | 41 |
|
42 | 42 | ############################################################################### |
43 | 43 | # Then, we compute the explainable variance per voxel. |
44 | 44 | # The variance of the signal is estimated by taking the average variance over |
45 | 45 | # repeats. The variance of the component shared across repeats is estimated by |
46 | 46 | # taking the variance of the average response. Then, we compute the |
47 | 47 | # explainable variance by dividing these two quantities. |
48 | | -# Finally, an correction can be applied to account for small numbers of repeat. |
| 48 | +# Finally, a correction can be applied to account for small numbers of repeat |
| 49 | +# (parameter ``bias_correction``). |
49 | 50 |
|
50 | 51 | from voxelwise.utils import explainable_variance |
51 | | - |
52 | 52 | ev = explainable_variance(Y_test, bias_correction=False) |
53 | 53 |
|
54 | 54 | ############################################################################### |
55 | | -# Plot the distribution of explainable variance over voxels. |
| 55 | +# We can plot the distribution of explainable variance over voxels. |
| 56 | + |
56 | 57 | import matplotlib.pyplot as plt |
57 | 58 |
|
58 | | -plt.hist(ev, bins=np.linspace(0, 1, 100), log=True, |
59 | | - histtype='step') |
| 59 | +plt.hist(ev, bins=np.linspace(0, 1, 100), log=True, histtype='step') |
60 | 60 | plt.xlabel("Explainable variance") |
61 | 61 | plt.ylabel("Number of voxels") |
62 | 62 | plt.title('Histogram of explainable variance') |
63 | 63 | plt.grid('on') |
64 | 64 | plt.show() |
65 | 65 |
|
| 66 | +############################################################################### |
| 67 | +# We see that most voxels have a rather low explainable variance, around 0.1 |
| 68 | +# (when not using the bias correction). This is expected, since most voxels are |
| 69 | +# not directly driven by a visual stimulus. |
| 70 | +# We also see that some voxels reach an explainable variance of 0.7, which is |
| 71 | +# quite high. It means that these voxels consistently record the same activity |
| 72 | +# across a repeated stimulus, and thus are good targets for encoding models. |
| 73 | + |
66 | 74 | ############################################################################### |
67 | 75 | # Map to subject flatmap |
68 | 76 | # ---------------------- |
| 77 | +# |
| 78 | +# To better understand the distribution of explainable variance, we map the |
| 79 | +# values to the subject brain. This can be done with |
| 80 | +# `pycortex <https://gallantlab.github.io/pycortex/>`_, which can create |
| 81 | +# interactive 3D viewers displayed in any modern browser. |
| 82 | +# ``Pycortex`` can also display flatten maps of the cortical surface, to |
| 83 | +# visualize the entire cortical surface at once. |
| 84 | +# |
| 85 | +# Here, we do not share the anatomical information of the subjects for privacy |
| 86 | +# concerns. Instead, we provide two mappers, (i) to map the voxels to a |
| 87 | +# subject-specific flatmap, or (ii) to map the voxels to the Freesurfer average |
| 88 | +# cortical surface ("fsaverage"). |
| 89 | +# |
| 90 | +# The first mapper is a sparse CSR matrix that map each voxel to a set of pixel |
| 91 | +# in a flatmap. To ease its use, we provide here an example function |
| 92 | +# ``plot_flatmap_from_mapper``. |
| 93 | + |
69 | 94 | from voxelwise.viz import plot_flatmap_from_mapper |
70 | 95 |
|
71 | 96 | mapper_file = os.path.join(directory, 'mappers', f'{subject}_mappers.hdf') |
|
74 | 99 |
|
75 | 100 | ############################################################################### |
76 | 101 | # We can see that the explainable variance is mainly located in the visual |
77 | | -# cortex, which was expected since this is a purely visual experiment. |
| 102 | +# cortex, in early regions like V1, V2, V3, or in higher-level regions like |
| 103 | +# EBA, FFA or IPS. This was expected since this is a purely visual experiment. |
| 104 | + |
| 105 | +############################################################################### |
| 106 | +# Map to fsaverage |
| 107 | +# ---------------- |
| 108 | +# |
| 109 | +# The second mapper we provide maps the voxel data to a Freesurfer |
| 110 | +# average surface ("fsaverage"), that can be used in ``pycortex``. |
| 111 | +# First, let's download the fsaverage surface if it does not exist |
| 112 | + |
| 113 | +import cortex |
| 114 | + |
| 115 | +surface = "fsaverage_pycortex" # ("fsaverage" outside the Gallant lab) |
| 116 | + |
| 117 | +if not hasattr(cortex.db, surface): |
| 118 | + cortex.utils.download_subject(subject_id=surface) |
| 119 | + |
| 120 | +############################################################################### |
| 121 | +# Then, we load the fsaverage mapper. The mapper is a sparse CSR matrix, which |
| 122 | +# map each voxel to some vertices in the fsaverage surface. |
| 123 | +# The mapper is applied with a dot product ``@``. |
| 124 | +from voxelwise.io import load_hdf5_sparse_array |
| 125 | +voxel_to_fsaverage = load_hdf5_sparse_array(mapper_file, |
| 126 | + key='voxel_to_fsaverage') |
| 127 | +ev_projected = voxel_to_fsaverage @ ev |
| 128 | + |
| 129 | +############################################################################### |
| 130 | +# We can then create a ``Vertex`` object with the projected data. |
| 131 | +# This object can be used either in a ``pycortex`` interactive 3D viewer, or |
| 132 | +# in a ``matplotlib`` figure showing directly the flatmap. |
| 133 | + |
| 134 | +vertex = cortex.Vertex(ev_projected, surface, vmin=0, vmax=0.7, cmap='inferno') |
| 135 | + |
| 136 | +############################################################################### |
| 137 | +# To start an interactive 3D viewer in the browser, use the following function: |
| 138 | +if False: |
| 139 | + cortex.webshow(vertex, open_browser=True) |
| 140 | + |
| 141 | +############################################################################### |
| 142 | +# Alternatively, to plot a flatmap in a ``matplotlib`` figure, use the |
| 143 | +# following function: |
| 144 | + |
| 145 | +fig = cortex.quickshow(vertex, colorbar_location='right') |
| 146 | +plt.show() |
0 commit comments