ENH add explainable_variance function

Author: TomDLT

tutorials/movies_3T/01_plot_explainable_variance.py

@@ -41,20 +41,21 @@
 
 ###############################################################################
 # Then, we compute the explainable variance per voxel.
-# The variance of the signal is estimated by taking the variance over
-# time (``var``). The variance of the component shared across repeats
-# is estimated by taking the variance of the average response (``var_mean``).
-# Then, we can compute the explainable variance by dividing the two quantities.
+# The variance of the signal is estimated by taking the average variance over
+# repeats. The variance of the component shared across repeats is estimated by
+# taking the variance of the average response. Then, we compute the
+# explainable variance by dividing these two quantities.
+# Finally, an correction can be applied to account for small numbers of repeat.
 
-var = np.var(Y_test.reshape(-1, Y_test.shape[-1]), axis=0)
-var_mean = np.var(np.mean(Y_test, axis=0), axis=0)
-explainable_variance = var_mean / var
+from voxelwise.utils import explainable_variance
+
+ev = explainable_variance(Y_test, bias_correction=False)
 
 ###############################################################################
 # Plot the distribution of explainable variance over voxels.
 import matplotlib.pyplot as plt
 
-plt.hist(explainable_variance, bins=np.linspace(0, 1, 100), log=True,
+plt.hist(ev, bins=np.linspace(0, 1, 100), log=True,
          histtype='step')
 plt.xlabel("Explainable variance")
 plt.ylabel("Number of voxels")
@@ -68,7 +69,7 @@
 from voxelwise.viz import plot_flatmap_from_mapper
 
 mapper_file = os.path.join(directory, 'mappers', f'{subject}_mappers.hdf')
-plot_flatmap_from_mapper(explainable_variance, mapper_file, vmin=0, vmax=0.7)
+plot_flatmap_from_mapper(ev, mapper_file, vmin=0, vmax=0.7)
 plt.show()
 
 ###############################################################################

tutorials/movies_3T/README.rst

@@ -23,13 +23,14 @@ This tutorial implements different voxelwise models:
 **Requirements:**
 This tutorial requires the following Python packages:
 
-- numpy  (for the data array)
-- scipy  (for motion energy extraction)
-- h5py  (for loading the data files)
-- scikit-learn  (for preprocessing and modeling)
-- himalaya  (for modeling)
-- pymoten  (for extracting motion energy)
+- numpy
+- scipy
+- h5py
+- scikit-learn
 - voxelwise  (this repository)
+- himalaya
+- pymoten  (optional, for extracting motion energy)
+- cupy/pytorch (optional, to use GPU in himalaya)
 
 
 **References:**

tutorials/movies_4T/README.rst

@@ -26,13 +26,14 @@ publication [5]_, and the CRCNS data set [6]_.
 **Requirements:**
 This tutorial requires the following Python packages:
 
-- numpy  (for the data array)
-- scipy  (for motion energy extraction)
-- h5py  (for loading the data files)
-- scikit-learn  (for preprocessing and modeling)
-- himalaya  (for modeling)
-- pymoten  (for extracting motion energy)
+- numpy
+- scipy
+- h5py
+- scikit-learn
 - voxelwise  (this repository)
+- himalaya
+- pymoten  (optional, for extracting motion energy)
+- cupy/pytorch (optional, to use GPU in himalaya)
 
 **References:**
 

voxelwise/utils.py

@@ -1,4 +1,5 @@
 import numpy as np
+import scipy.stats
 from sklearn.utils.validation import check_random_state
 
 
@@ -43,3 +44,34 @@ def generate_leave_one_run_out(n_samples, run_onsets, random_state=None,
             [runs[jj] for jj in range(n_runs) if jj not in val_runs])
         val = np.hstack([runs[jj] for jj in range(n_runs) if jj in val_runs])
         yield train, val
+
+
+def explainable_variance(data, bias_correction=True, do_zscore=True):
+    """Compute explainable variance for a set of voxels.
+
+    Parameters
+    ----------
+    data : array of shape (n_repeats, n_times, n_voxels)
+        fMRI reponses of the repeated test set.
+    bias_correction: bool
+        Perform bias correction based on the number of repetitions.
+    do_zscore: bool
+        z-score the data in time. Only set to False if your data time courses
+        are already z-scored.
+
+    Returns
+    -------
+    ev : array of shape (n_voxels, )
+        Explainable variance per voxel.
+    """
+    if do_zscore:
+        data = scipy.stats.zscore(data, axis=1)
+
+    mean_var = data.var(axis=1, dtype=np.float64, ddof=1).mean(axis=0)
+    var_mean = data.mean(axis=0).var(axis=0, dtype=np.float64, ddof=1)
+    ev = var_mean / mean_var
+
+    if bias_correction:
+        n_repeats = data.shape[0]
+        ev = ev - (1 - ev) / (n_repeats - 1)
+    return ev