|
| 1 | +""" |
| 2 | +Predicted time series and residuals |
| 3 | +=================================== |
| 4 | +
|
| 5 | +Here we fit a First Level GLM with the `minimize_memory`-argument set to `False`. |
| 6 | +By doing so, the `FirstLevelModel`-object stores the residuals, which we can then inspect. |
| 7 | +Also, the predicted time series can be extracted, which is useful to assess the quality of the model fit. |
| 8 | +""" |
| 9 | + |
| 10 | + |
| 11 | +######################################################################### |
| 12 | +# Import modules |
| 13 | +# -------------- |
| 14 | +from nistats.datasets import fetch_spm_auditory |
| 15 | +from nilearn import image |
| 16 | +from nilearn import masking |
| 17 | +import pandas as pd |
| 18 | + |
| 19 | + |
| 20 | +# load fMRI data |
| 21 | +subject_data = fetch_spm_auditory() |
| 22 | +fmri_img = image.concat_imgs(subject_data.func) |
| 23 | + |
| 24 | +# Make an average |
| 25 | +mean_img = image.mean_img(fmri_img) |
| 26 | +mask = masking.compute_epi_mask(mean_img) |
| 27 | + |
| 28 | +# Clean and smooth data |
| 29 | +fmri_img = image.clean_img(fmri_img, standardize=False) |
| 30 | +fmri_img = image.smooth_img(fmri_img, 5.) |
| 31 | + |
| 32 | +# load events |
| 33 | +events = pd.read_table(subject_data['events']) |
| 34 | + |
| 35 | + |
| 36 | +######################################################################### |
| 37 | +# Fit model |
| 38 | +# --------- |
| 39 | +# Note that `minimize_memory` is set to `False` so that `FirstLevelModel` |
| 40 | +# stores the residuals. |
| 41 | +# `signal_scaling` is set to False, so we keep the same scaling as the |
| 42 | +# original data in `fmri_img`. |
| 43 | +from nistats.first_level_model import FirstLevelModel |
| 44 | + |
| 45 | +fmri_glm = FirstLevelModel(t_r=7, |
| 46 | + drift_model='cosine', |
| 47 | + signal_scaling=False, |
| 48 | + mask_img=mask, |
| 49 | + minimize_memory=False) |
| 50 | + |
| 51 | +fmri_glm = fmri_glm.fit(fmri_img, events) |
| 52 | + |
| 53 | + |
| 54 | +######################################################################### |
| 55 | +# Calculate and plot contrast |
| 56 | +# --------------------------- |
| 57 | +from nilearn import plotting |
| 58 | + |
| 59 | +z_map = fmri_glm.compute_contrast('active - rest') |
| 60 | + |
| 61 | +plotting.plot_stat_map(z_map, bg_img=mean_img, threshold=3.1) |
| 62 | + |
| 63 | +######################################################################### |
| 64 | +# Extract the largest clusters |
| 65 | +# ---------------------------- |
| 66 | +from nistats.reporting import get_clusters_table |
| 67 | +from nilearn import input_data |
| 68 | + |
| 69 | +table = get_clusters_table(z_map, stat_threshold=3.1, |
| 70 | + cluster_threshold=20).set_index('Cluster ID', drop=True) |
| 71 | +table.head() |
| 72 | + |
| 73 | +# get the 6 largest clusters' max x, y, and z coordinates |
| 74 | +coords = table.loc[range(1, 7), ['X', 'Y', 'Z']].values |
| 75 | + |
| 76 | +# extract time series from each coordinate |
| 77 | +masker = input_data.NiftiSpheresMasker(coords) |
| 78 | +real_timeseries = masker.fit_transform(fmri_img) |
| 79 | +predicted_timeseries = masker.fit_transform(fmri_glm.predicted[0]) |
| 80 | + |
| 81 | + |
| 82 | +######################################################################### |
| 83 | +# Plot predicted and actual time series for 6 most significant clusters |
| 84 | +# --------------------------------------------------------------------- |
| 85 | +import matplotlib.pyplot as plt |
| 86 | + |
| 87 | +# colors for each of the clusters |
| 88 | +colors = ['blue', 'navy', 'purple', 'magenta', 'olive', 'teal'] |
| 89 | +# plot the time series and corresponding locations |
| 90 | +fig1, axs1 = plt.subplots(2, 6) |
| 91 | +for i in range(0, 6): |
| 92 | + # plotting time series |
| 93 | + axs1[0, i].set_title('Cluster peak {}\n'.format(coords[i])) |
| 94 | + axs1[0, i].plot(real_timeseries[:, i], c=colors[i], lw=2) |
| 95 | + axs1[0, i].plot(predicted_timeseries[:, i], c='r', ls='--', lw=2) |
| 96 | + axs1[0, i].set_xlabel('Time') |
| 97 | + axs1[0, i].set_ylabel('Signal intensity', labelpad=0) |
| 98 | + # plotting image below the time series |
| 99 | + roi_img = plotting.plot_stat_map( |
| 100 | + z_map, cut_coords=[coords[i][2]], threshold=3.1, figure=fig1, |
| 101 | + axes=axs1[1, i], display_mode='z', colorbar=False, bg_img=mean_img) |
| 102 | + roi_img.add_markers([coords[i]], colors[i], 300) |
| 103 | +fig1.set_size_inches(24, 14) |
| 104 | + |
| 105 | + |
| 106 | +######################################################################### |
| 107 | +# Get residuals |
| 108 | +# ------------- |
| 109 | +resid = masker.fit_transform(fmri_glm.residuals[0]) |
| 110 | + |
| 111 | + |
| 112 | +######################################################################### |
| 113 | +# Plot distribution of residuals |
| 114 | +# ------------------------------ |
| 115 | +# Note that residuals are not really distributed normally. |
| 116 | +fig2, axs2 = plt.subplots(2, 3) |
| 117 | +axs2 = axs2.flatten() |
| 118 | +for i in range(0, 6): |
| 119 | + axs2[i].set_title('Cluster peak {}\n'.format(coords[i])) |
| 120 | + axs2[i].hist(resid[:, i], color=colors[i]) |
| 121 | + print('Mean residuals: {}'.format(resid[:, i].mean())) |
| 122 | + |
| 123 | +fig2.set_size_inches(12, 7) |
| 124 | +fig2.tight_layout() |
| 125 | + |
| 126 | + |
| 127 | +######################################################################### |
| 128 | +# Plot R-squared |
| 129 | +# -------------- |
| 130 | +# Because we stored the residuals, we can plot the R-squared: the proportion |
| 131 | +# of explained variance of the GLM as a whole. Note that the R-squared is markedly |
| 132 | +# lower deep down the brain, where there is more physiological noise and we |
| 133 | +# are further away from the receive coils. However, R-Squared should be interpreted |
| 134 | +# with a grain of salt. The R-squared value will necessarily increase with |
| 135 | +# the addition of more factors (such as rest, active, drift, motion) into the GLM. |
| 136 | +# Additionally, we are looking at the overall fit of the model, so we are |
| 137 | +# unable to say whether a voxel/region has a large R-squared value because |
| 138 | +# the voxel/region is responsive to the experiment (such as active or rest) |
| 139 | +# or because the voxel/region fits the noise factors (such as drift or motion) |
| 140 | +# that could be present in the GLM. To isolate the influence of the experiment, |
| 141 | +# we can use an F-test as shown in the next section. |
| 142 | +plotting.plot_stat_map(fmri_glm.r_square[0], |
| 143 | + bg_img=mean_img, threshold=.1, display_mode='z', cut_coords=7) |
| 144 | + |
| 145 | + |
| 146 | +######################################################################### |
| 147 | +# Calculate and Plot F-test |
| 148 | +# ------------------------- |
| 149 | +# The F-test tells you how well the GLM fits effects of interest such as |
| 150 | +# the active and rest conditions together. This is different from R-squared, |
| 151 | +# which tells you how well the overall GLM fits the data, including active, rest |
| 152 | +# and all the other columns in the design matrix such as drift and motion. |
| 153 | +import numpy as np |
| 154 | + |
| 155 | +design_matrix = fmri_glm.design_matrices_[0] |
| 156 | + |
| 157 | +# contrast with a one for "active" and zero everywhere else |
| 158 | +active = np.array([1 if c == 'active' else 0 for c in design_matrix.columns]) |
| 159 | + |
| 160 | +# contrast with a one for "rest" and zero everywhere else |
| 161 | +rest = np.array([1 if c == 'rest' else 0 for c in design_matrix.columns]) |
| 162 | + |
| 163 | +effects_of_interest = np.vstack((active, rest)) |
| 164 | +# f-test for rest and activity |
| 165 | +z_map_ftest = fmri_glm.compute_contrast( |
| 166 | + effects_of_interest, |
| 167 | + stat_type='F', |
| 168 | + output_type='z_score') |
| 169 | + |
| 170 | +plotting.plot_stat_map(z_map_ftest, |
| 171 | + bg_img=mean_img, threshold=3.1, display_mode='z', cut_coords=7) |
0 commit comments