diff --git a/brainiak/isfc.py b/brainiak/isfc.py index 176042d8b..32bb96f7d 100644 --- 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 index 4c6772293..4b60497e7 100644 --- 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 new file mode 100644 index 000000000..d07a06846 --- /dev/null +++ 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 index cfe2163c6..e4d516be0 100644 --- 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,7 +42,8 @@ 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], @@ -43,5 +51,10 @@ def test_ISFC(): [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" diff --git a/tests/utils/test_utils.py b/tests/utils/test_utils.py index d6b468a73..f58a2768c 100644 --- a/tests/utils/test_utils.py +++ b/tests/utils/test_utils.py @@ -173,3 +173,62 @@ def test_center_mass_exp(): assert np.isclose(result, 1.0+1e-10), 'for a small '\ 'enough interval, the center of mass should be '\ 'close to its mid-point' + + +def test_phase_randomize(): + from brainiak.utils.utils import phase_randomize + import numpy as np + from scipy.fftpack import fft + import math + from scipy.stats import pearsonr + + # Generate auto-correlated signals + nv = 2 + T = 100 + ns = 3 + D = np.zeros((nv, T, ns)) + for v in range(nv): + for s in range(ns): + D[v, :, s] = np.sin(np.linspace(0, math.pi * 5 * (v + 1), T)) + \ + np.sin(np.linspace(0, math.pi * 6 * (s + 1), T)) + + freq = fft(D, axis=1) + D_pr = phase_randomize(D) + freq_pr = fft(D_pr, axis=1) + p_corr = pearsonr(np.angle(freq).flatten(), np.angle(freq_pr).flatten())[0] + + assert np.isclose(abs(freq), abs(freq_pr)).all(), \ + "Amplitude spectrum not preserved under phase randomization" + + assert abs(p_corr) < 0.03, \ + "Phases still correlated after randomization" + + +def test_ecdf(): + from brainiak.utils.utils import ecdf + import numpy as np + + x = np.array([1, 4, 3]) + cdf_fun = ecdf(x) + + assert np.isclose(cdf_fun(0), [0]), "Left side of cdf should be 0" + assert np.isclose(cdf_fun(5), [1]), "Right side of cdf should be 1" + assert np.isclose(cdf_fun(1.5), [1 / 3]), "CDF value incorrect" + assert np.isclose(cdf_fun(1), [1 / 3]), "CDF should be right-continuous" + + +def test_p_from_null(): + from brainiak.utils.utils import p_from_null + import numpy as np + + X = np.zeros((2, 5)) # One true value, 4 null values + X[0, 0] = 1 + X[1, 0] = -2 + X[0, 1:] = [-1.0, 0.00, 0.50, 2.00] + X[1, 1:] = [-1.5, 0.25, -0.25, 0.25] + + p_1side = p_from_null(X, two_sided=False) + assert np.isclose(p_1side, [0.25, 1]).all(), "One-sided p value incorrect" + + p_2side = p_from_null(X, two_sided=True) + assert np.isclose(p_2side, [0.5, 0]).all(), "Two-sided p value incorrect"