mihaic · Dec 13, 2017 · Dec 11, 2017 · Dec 11, 2017 · Dec 12, 2017 · Dec 12, 2017
diff --git a/brainiak/isfc.py b/brainiak/isfc.py
@@ -14,8 +14,7 @@
 """Intersubject analyses (ISC/ISFC)
 
 Functions for computing intersubject correlation (ISC) and intersubject
-functional correlation (ISFC), and utility functions for loading subject data
-files and masks
+functional correlation (ISFC)
 
 Paper references:
 ISC: Hasson, U., Nir, Y., Levy, I., Fuhrmann, G. & Malach, R. Intersubject
@@ -33,14 +32,18 @@
 from brainiak.fcma.util import compute_correlation
 import numpy as np
 from scipy import stats
+from .utils.utils import phase_randomize, p_from_null
 
 
-def isc(D, collapse_subj=True):
+def isc(D, collapse_subj=True, return_p=False, num_perm=1000,
+        two_sided=False, random_state=0):
     """Intersubject correlation
 
     For each voxel, computes the correlation of each subject's timecourse with
     the mean of all other subjects' timecourses. By default the result is
-    averaged across subjects, unless collapse_subj is set to False.
+    averaged across subjects, unless collapse_subj is set to False. A null
+    distribution can optionally be computed using phase randomization, to
+    compute a p value for each voxel.
 
     Parameters
     ----------
@@ -50,35 +53,72 @@ def isc(D, collapse_subj=True):
     collapse_subj : bool, default:True
         Whether to average across subjects before returning result
 
+    return_p : bool, default:False
+        Whether to use phase randomization to compute a p value for each voxel
+
+    num_perm : int, default:1000
+        Number of null samples to use for computing p values
+
+    two_sided : bool, default:False
+        Whether the p value should be one-sided (testing only for being
+        above the null) or two-sided (testing for both significantly positive
+        and significantly negative values)
+
+    random_state : RandomState or an int seed (0 by default)
+        A random number generator instance to define the state of the
+        random permutations generator.
+
     Returns
     -------
     ISC : voxel ndarray (or voxel by subject ndarray, if collapse_subj=False)
         pearson correlation for each voxel, across subjects
+
+    p : ndarray the same shape as ISC (if return_p = True)
+        p values for each ISC value under the null distribution
     """
 
     n_vox = D.shape[0]
     n_subj = D.shape[2]
-    ISC = np.zeros((n_vox, n_subj))
 
-    # Loop across choice of leave-one-out subject
-    for loo_subj in range(n_subj):
-        group = np.mean(D[:, :, np.arange(n_subj) != loo_subj], axis=2)
-        subj = D[:, :, loo_subj]
-        for v in range(n_vox):
-            ISC[v, loo_subj] = stats.pearsonr(group[v, :], subj[v, :])[0]
+    if return_p:
+        n_perm = num_perm
+    else:
+        n_perm = 0
+
+    ISC = np.zeros((n_vox, n_subj, n_perm + 1))
+
+    for p in range(n_perm + 1):
+        # Loop across choice of leave-one-out subject
+        for loo_subj in range(n_subj):
+            group = np.mean(D[:, :, np.arange(n_subj) != loo_subj], axis=2)
+            subj = D[:, :, loo_subj]
+            for v in range(n_vox):
+                ISC[v, loo_subj, p] = stats.pearsonr(group[v, :],
+                                                     subj[v, :])[0]
+
+        # Randomize phases of D to create next null dataset
+        D = phase_randomize(D, random_state)
 
     if collapse_subj:
         ISC = np.mean(ISC, axis=1)
-    return ISC
+
+    if return_p:
+        p = p_from_null(ISC, two_sided)
+        return ISC[..., 0], p
+    else:
+        return ISC[..., 0]
 
 
-def isfc(D, collapse_subj=True):
+def isfc(D, collapse_subj=True, return_p=False,
+         num_perm=1000, two_sided=False, random_state=0):
     """Intersubject functional correlation
 
     Computes the correlation between the timecoure of each voxel in each
     subject with the average of all other subjects' timecourses in *all*
     voxels. By default the result is averaged across subjects, unless
-    collapse_subj is set to False.
+    collapse_subj is set to False. A null distribution can optionally be
+    computed using phase randomization, to compute a p value for each voxel-to-
+    voxel correlation.
 
     Uses the high performance compute_correlation routine from fcma.util
 
@@ -90,27 +130,60 @@ def isfc(D, collapse_subj=True):
     collapse_subj : bool, default:True
         Whether to average across subjects before returning result
 
+    return_p : bool, default:False
+        Whether to use phase randomization to compute a p value for each voxel
+
+    num_perm : int, default:1000
+        Number of null samples to use for computing p values
+
+    two_sided : bool, default:False
+        Whether the p value should be one-sided (testing only for being
+        above the null) or two-sided (testing for both significantly positive
+        and significantly negative values)
+
+    random_state : RandomState or an int seed (0 by default)
+        A random number generator instance to define the state of the
+        random permutations generator.
+
     Returns
     -------
     ISFC : voxel by voxel ndarray
         (or voxel by voxel by subject ndarray, if collapse_subj=False)
         pearson correlation between all pairs of voxels, across subjects
+
+    p : ndarray the same shape as ISC (if return_p = True)
+        p values for each ISC value under the null distribution
     """
 
     n_vox = D.shape[0]
     n_subj = D.shape[2]
-    ISFC = np.zeros((n_vox, n_vox, n_subj))
 
-    # Loop across choice of leave-one-out subject
-    for loo_subj in range(D.shape[2]):
-        group = np.mean(D[:, :, np.arange(n_subj) != loo_subj], axis=2)
-        subj = D[:, :, loo_subj]
-        ISFC[:, :, loo_subj] = compute_correlation(group, subj)
+    if return_p:
+        n_perm = num_perm
+    else:
+        n_perm = 0
 
-        # Symmetrize matrix
-        ISFC[:, :, loo_subj] = (ISFC[:, :, loo_subj] +
-                                ISFC[:, :, loo_subj].T) / 2
+    ISFC = np.zeros((n_vox, n_vox, n_subj, n_perm + 1))
+
+    for p in range(n_perm + 1):
+        # Loop across choice of leave-one-out subject
+        for loo_subj in range(D.shape[2]):
+            group = np.mean(D[:, :, np.arange(n_subj) != loo_subj], axis=2)
+            subj = D[:, :, loo_subj]
+            ISFC[:, :, loo_subj, p] = compute_correlation(group, subj)
+
+            # Symmetrize matrix
+            ISFC[:, :, loo_subj, p] = (ISFC[:, :, loo_subj, p] +
+                                       ISFC[:, :, loo_subj, p].T) / 2
+
+        # Randomize phases of D to create next null dataset
+        D = phase_randomize(D, random_state)
 
     if collapse_subj:
         ISFC = np.mean(ISFC, axis=2)
-    return ISFC
+
+    if return_p:
+        p = p_from_null(ISFC, two_sided)
+        return ISFC[..., 0], p
+    else:
+        return ISFC[..., 0]
diff --git a/brainiak/utils/utils.py b/brainiak/utils/utils.py
@@ -17,6 +17,10 @@
 import os.path
 import psutil
 from .fmrisim import generate_stimfunction, _double_gamma_hrf, convolve_hrf
+from sklearn.utils import check_random_state
+from scipy.fftpack import fft, ifft
+import math
+
 
 """
 Some utility functions that can be used by different algorithms
@@ -667,3 +671,125 @@ def usable_cpu_count():
         except AttributeError:
             result = os.cpu_count()
     return result
+
+
+def phase_randomize(D, random_state=0):
+    """Randomly shift signal phases
+
+    For each timecourse (from each voxel and each subject), computes its DFT
+    and then randomly shifts the phase of each frequency before inverting
+    back into the time domain. This yields timecourses with the same power
+    spectrum (and thus the same autocorrelation) as the original timecourses,
+    but will remove any meaningful temporal relationships between the
+    timecourses.
+
+    This procedure is described in:
+    Simony E, Honey CJ, Chen J, Lositsky O, Yeshurun Y, Wiesel A, Hasson U
+    (2016) Dynamic reconfiguration of the default mode network during narrative
+    comprehension. Nat Commun 7.
+
+    Parameters
+    ----------
+    D : voxel by time by subject ndarray
+        fMRI data to be phase randomized
+
+    random_state : RandomState or an int seed (0 by default)
+        A random number generator instance to define the state of the
+        random permutations generator.
+
+    Returns
+    ----------
+    ndarray of same shape as D
+        phase randomized timecourses
+    """
+
+    random_state = check_random_state(random_state)
+
+    F = fft(D, axis=1)
+    if D.shape[1] % 2 == 0:
+        pos_freq = np.arange(1, D.shape[1] // 2)
+        neg_freq = np.arange(D.shape[1] - 1, D.shape[1] // 2, -1)
+    else:
+        pos_freq = np.arange(1, (D.shape[1] - 1) // 2 + 1)
+        neg_freq = np.arange(D.shape[1] - 1, (D.shape[1] - 1) // 2, -1)
+
+    shift = random_state.rand(D.shape[0], len(pos_freq),
+                              D.shape[2]) * 2 * math.pi
+
+    # Shift pos and neg frequencies symmetrically, to keep signal real
+    F[:, pos_freq, :] *= np.exp(1j * shift)
+    F[:, neg_freq, :] *= np.exp(-1j * shift)
+
+    return np.real(ifft(F, axis=1))
+
+
+def ecdf(x):
+    """Empirical cumulative distribution function
+
+    Given a 1D array of values, returns a function f(q) that outputs the
+    fraction of values less than or equal to q.
+
+    Parameters
+    ----------
+    x : 1D array
+        values for which to compute CDF
+
+    Returns
+    ----------
+    ecdf_fun: Callable[[float], float]
+        function that returns the value of the CDF at a given point
+    """
+    xp = np.sort(x)
+    yp = np.arange(len(xp) + 1) / len(xp)
+
+    def ecdf_fun(q):
+        return yp[np.searchsorted(xp, q, side="right")]
+
+    return ecdf_fun
+
+
+def p_from_null(X, two_sided=False):
+    """Compute p value of true result from null distribution
+
+    Given an array containing both a real result and a set of null results,
+    computes the fraction of null results larger than the real result (or,
+    if two_sided=True, the fraction of null results more extreme than the real
+    result in either the positive or negative direction).
+
+    Note that all real results are compared to a pooled null distribution,
+    which is the max/min over all null results, providing multiple
+    comparisons correction.
+
+    Parameters
+    ----------
+    X : ndarray with arbitrary number of dimensions
+        The last dimension of X should contain the real result in X[..., 0]
+        and the null results in X[..., 1:]
+
+    two_sided : bool, default:False
+        Whether the p value should be one-sided (testing only for being
+        above the null) or two-sided (testing for both significantly positive
+        and significantly negative values)
+
+    Returns
+    -------
+    p : ndarray the same shape as X, without the last dimension
+        p values for each true X value under the null distribution
+    """
+    leading_dims = tuple([int(d) for d in np.arange(X.ndim - 1)])
+
+    # Compute maximum/minimum in each null dataset
+    max_null = np.max(X[..., 1:], axis=leading_dims)
+    min_null = np.min(X[..., 1:], axis=leading_dims)
+
+    # Compute where the true values fall on the null distribution
+    max_null_ecdf = ecdf(max_null)
+    if two_sided:
+        min_null_ecdf = ecdf(min_null)
+        p = 2 * np.minimum(1 - max_null_ecdf(X[..., 0]),
+                           min_null_ecdf(X[..., 0]))
+        p = np.minimum(p, 1)
+    else:
+        p = 1 - max_null_ecdf(X[..., 0])
+
+    return p
diff --git a/docs/newsfragments/310.feature b/docs/newsfragments/310.feature
@@ -0,0 +1,2 @@
+Added option to compute p values for ISC and ISFC
+
diff --git a/tests/isfc/test_isfc.py b/tests/isfc/test_isfc.py
@@ -5,21 +5,28 @@
 
 
 def test_ISC():
+    # Create dataset in which one voxel is highly correlated across subjects
+    # and the other is not
     D = np.zeros((2, 5, 3))
     D[:, :, 0] = \
         [[-0.36225433, -0.43482456,  0.26723158,  0.16461712, -0.37991465],
          [-0.62305959, -0.46660116, -0.50037994,  1.81083754,  0.23499509]]
     D[:, :, 1] = \
-        [[-0.31484153, -0.47486988,  0.18966625, -0.12568572, -0.40535156],
-         [-0.78433298, -0.35875085, -0.6121344,  1.68267639,  0.28603493]]
+        [[-0.30484153, -0.49486988,  0.10966625, -0.19568572, -0.20535156],
+         [1.68267639, -0.78433298, -0.35875085, -0.6121344,  0.28603493]]
     D[:, :, 2] = \
-        [[-0.26593192, -0.56914734,  0.28397317,  0.30276589, -0.42637472],
-         [-0.67598379, -0.51549055, -0.64196342,  1.60686666,  0.04127293]]
+        [[-0.36593192, -0.50914734,  0.21397317,  0.30276589, -0.42637472],
+         [0.04127293, -0.67598379, -0.51549055, -0.64196342,  1.60686666]]
 
-    ISC = brainiak.isfc.isc(D)
-    assert np.isclose(ISC, [0.9540602, 0.99585304]).all(), \
+    (ISC, p) = brainiak.isfc.isc(D, return_p=True, num_perm=100,
+                                 two_sided=True, random_state=0)
+
+    assert np.isclose(ISC, [0.8909243, 0.0267954]).all(), \
         "Calculated ISC does not match ground truth"
 
+    assert np.isclose(p, [0.02, 1]).all(), \
+        "Calculated p values do not match ground truth"
+
 
 def test_ISFC():
     curr_dir = os.path.dirname(__file__)
@@ -35,13 +42,19 @@ def test_ISFC():
 
     assert D.shape == (4, 5, 2), "Loaded data has incorrect shape"
 
-    ISFC = brainiak.isfc.isfc(D)
+    (ISFC, p) = brainiak.isfc.isfc(D, return_p=True, num_perm=100,
+                                   two_sided=True, random_state=0)
 
     ground_truth = \
         [[1, 1, 0, -1],
          [1, 1, 0, -1],
          [0, 0, 1,  0],
          [-1, -1, 0, 1]]
 
+    ground_truth_p = 1 - np.abs(ground_truth)
+
     assert np.isclose(ISFC, ground_truth).all(), \
         "Calculated ISFC does not match ground truth"
+
+    assert np.isclose(p, ground_truth_p).all(), \
+        "Calculated p values do not match ground truth"
Original file line number	Diff line number	Diff line change
		@@ -0,0 +1,2 @@
		Added option to compute p values for ISC and ISFC