From c4dbb9595446cc9533e3cca214e0d1b0367da66b Mon Sep 17 00:00:00 2001 From: Christopher Baldassano Date: Mon, 11 Dec 2017 16:01:20 -0500 Subject: [PATCH 1/5] Adding statistical tests to ISC and ISFC --- brainiak/isfc.py | 209 +++++++++++++++++++++++++++++++++++----- tests/isfc/test_isfc.py | 22 ++++- 2 files changed, 202 insertions(+), 29 deletions(-) diff --git a/brainiak/isfc.py b/brainiak/isfc.py index 176042d8b..6da6a04a1 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 scipy.fftpack import fft, ifft +import math -def isc(D, collapse_subj=True): +def isc(D, collapse_subj=True, return_p=False, num_perm=1000, two_sided=False): """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,85 @@ 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) + 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) if collapse_subj: ISC = np.mean(ISC, axis=1) - return ISC + + if not return_p: + return np.squeeze(ISC) + + # Compute maximum/minimum ISC in each null dataset + if collapse_subj: + max_null = np.max(ISC[:, 1:], axis=0) + min_null = np.min(ISC[:, 1:], axis=0) + ISC = ISC[:, 0] + else: + max_null = np.max(ISC[:, :, 1:], axis=(0, 1)) + min_null = np.min(ISC[:, :, 1:], axis=(0, 1)) + ISC = ISC[:, :, 0] + + # 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(ISC), min_null_ecdf(ISC)) + p = np.minimum(p, 1) + else: + p = 1 - max_null_ecdf(ISC) + return ISC, p -def isfc(D, collapse_subj=True): +def isfc(D, collapse_subj=True, return_p=False, + num_perm=1000, two_sided=False): """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 +143,135 @@ 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) + 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 + + 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 - # Symmetrize matrix - ISFC[:, :, loo_subj] = (ISFC[:, :, loo_subj] + - ISFC[:, :, loo_subj].T) / 2 + # Randomize phases of D to create next null dataset + D = phase_randomize(D) if collapse_subj: ISFC = np.mean(ISFC, axis=2) - return ISFC + + if not return_p: + return np.squeeze(ISFC) + + # Compute maximum/minimum ISFC in each null dataset + if collapse_subj: + max_null = np.max(ISFC[:, :, 1:], axis=(0, 1)) + min_null = np.min(ISFC[:, :, 1:], axis=(0, 1)) + ISFC = ISFC[:, :, 0] + else: + max_null = np.max(ISFC[:, :, :, 1:], axis=(0, 1, 2)) + min_null = np.min(ISFC[:, :, 1:], axis=(0, 1)) + ISFC = ISFC[:, :, :, 0] + + # 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(ISFC), min_null_ecdf(ISFC)) + p = np.minimum(p, 1) + else: + p = 1 - max_null_ecdf(ISFC) + return ISFC, p + + +def phase_randomize(D): + """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. + + Parameters + ---------- + D : voxel by time by subject ndarray + fMRI data to be phase randomized + + Returns + ---------- + ndarray of same shape as D + phase randomized timecourses + """ + + 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 = np.random.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 diff --git a/tests/isfc/test_isfc.py b/tests/isfc/test_isfc.py index cfe2163c6..533bb028a 100644 --- a/tests/isfc/test_isfc.py +++ b/tests/isfc/test_isfc.py @@ -11,15 +11,20 @@ def test_ISC(): [-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]] + [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.04127293, -0.67598379, -0.51549055, -0.64196342, 1.60686666]] - ISC = brainiak.isfc.isc(D) - assert np.isclose(ISC, [0.9540602, 0.99585304]).all(), \ + np.random.seed(0) + (ISC, p) = brainiak.isfc.isc(D, return_p=True, + num_perm=100, two_sided=True) + assert np.isclose(ISC, [0.9540602, 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 +40,9 @@ def test_ISFC(): assert D.shape == (4, 5, 2), "Loaded data has incorrect shape" - ISFC = brainiak.isfc.isfc(D) + np.random.seed(0) + (ISFC, p) = brainiak.isfc.isfc(D, return_p=True, + num_perm=100, two_sided=True) ground_truth = \ [[1, 1, 0, -1], @@ -43,5 +50,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" From 178d7df0c5e5facb7920899f90bc863d53074e1e Mon Sep 17 00:00:00 2001 From: Christopher Baldassano Date: Mon, 11 Dec 2017 16:07:14 -0500 Subject: [PATCH 2/5] Adding news snippet --- docs/newsfragments/310.feature | 2 ++ 1 file changed, 2 insertions(+) create mode 100644 docs/newsfragments/310.feature 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 + From 726d4aa99756c6b3595f65bb577736c641a85665 Mon Sep 17 00:00:00 2001 From: Christopher Baldassano Date: Tue, 12 Dec 2017 17:38:27 -0500 Subject: [PATCH 3/5] Code review tweaks, and refactoring p_from_null --- brainiak/isfc.py | 151 ++++++++++++++-------------------------- brainiak/utils/utils.py | 79 +++++++++++++++++++++ 2 files changed, 132 insertions(+), 98 deletions(-) diff --git a/brainiak/isfc.py b/brainiak/isfc.py index 6da6a04a1..820bca7b8 100644 --- a/brainiak/isfc.py +++ b/brainiak/isfc.py @@ -32,11 +32,11 @@ from brainiak.fcma.util import compute_correlation import numpy as np from scipy import stats -from scipy.fftpack import fft, ifft -import math +from .utils.utils import phase_randomize, ecdf -def isc(D, collapse_subj=True, return_p=False, num_perm=1000, two_sided=False): +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 @@ -64,6 +64,10 @@ def isc(D, collapse_subj=True, return_p=False, num_perm=1000, two_sided=False): 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) @@ -93,37 +97,20 @@ def isc(D, collapse_subj=True, return_p=False, num_perm=1000, two_sided=False): subj[v, :])[0] # Randomize phases of D to create next null dataset - D = phase_randomize(D) + D = phase_randomize(D, random_state) if collapse_subj: ISC = np.mean(ISC, axis=1) - if not return_p: - return np.squeeze(ISC) - - # Compute maximum/minimum ISC in each null dataset - if collapse_subj: - max_null = np.max(ISC[:, 1:], axis=0) - min_null = np.min(ISC[:, 1:], axis=0) - ISC = ISC[:, 0] - else: - max_null = np.max(ISC[:, :, 1:], axis=(0, 1)) - min_null = np.min(ISC[:, :, 1:], axis=(0, 1)) - ISC = ISC[:, :, 0] - - # 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(ISC), min_null_ecdf(ISC)) - p = np.minimum(p, 1) + if return_p: + p = p_from_null(ISC, two_sided) + return ISC[..., 0], p else: - p = 1 - max_null_ecdf(ISC) - return ISC, p + return ISC[..., 0] def isfc(D, collapse_subj=True, return_p=False, - num_perm=1000, two_sided=False): + num_perm=1000, two_sided=False, random_state=0): """Intersubject functional correlation Computes the correlation between the timecoure of each voxel in each @@ -154,6 +141,10 @@ def isfc(D, collapse_subj=True, return_p=False, 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 @@ -186,92 +177,56 @@ def isfc(D, collapse_subj=True, return_p=False, ISFC[:, :, loo_subj, p].T) / 2 # Randomize phases of D to create next null dataset - D = phase_randomize(D) + D = phase_randomize(D, random_state) if collapse_subj: ISFC = np.mean(ISFC, axis=2) - if not return_p: - return np.squeeze(ISFC) - - # Compute maximum/minimum ISFC in each null dataset - if collapse_subj: - max_null = np.max(ISFC[:, :, 1:], axis=(0, 1)) - min_null = np.min(ISFC[:, :, 1:], axis=(0, 1)) - ISFC = ISFC[:, :, 0] - else: - max_null = np.max(ISFC[:, :, :, 1:], axis=(0, 1, 2)) - min_null = np.min(ISFC[:, :, 1:], axis=(0, 1)) - ISFC = ISFC[:, :, :, 0] - - # 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(ISFC), min_null_ecdf(ISFC)) - p = np.minimum(p, 1) + if return_p: + p = p_from_null(ISFC, two_sided) + return ISFC[..., 0], p else: - p = 1 - max_null_ecdf(ISFC) - return ISFC, p + return ISFC[..., 0] -def phase_randomize(D): - """Randomly shift signal phases +def p_from_null(X, two_sided=False): + """Compute p value of true result from null distribution - 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. + 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). Parameters ---------- - D : voxel by time by subject ndarray - fMRI data to be phase randomized + 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 - ---------- - ndarray of same shape as D - phase randomized timecourses + ------- + 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)]) - 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 = np.random.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. + # 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) - 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) + # 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]) - def ecdf_fun(q): - return yp[np.searchsorted(xp, q, side="right")] - return ecdf_fun + return p diff --git a/brainiak/utils/utils.py b/brainiak/utils/utils.py index 4c6772293..ee5032c52 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,78 @@ 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 From 6a12d8cb00c6aca26a7c87eec0e45cb02844b0e5 Mon Sep 17 00:00:00 2001 From: Christopher Baldassano Date: Tue, 12 Dec 2017 17:42:58 -0500 Subject: [PATCH 4/5] Added updated tests --- tests/isfc/test_isfc.py | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/tests/isfc/test_isfc.py b/tests/isfc/test_isfc.py index 533bb028a..7eb4b5391 100644 --- a/tests/isfc/test_isfc.py +++ b/tests/isfc/test_isfc.py @@ -5,21 +5,22 @@ 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.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.36593192, -0.50914734, 0.21397317, 0.30276589, -0.42637472], [0.04127293, -0.67598379, -0.51549055, -0.64196342, 1.60686666]] - np.random.seed(0) - (ISC, p) = brainiak.isfc.isc(D, return_p=True, - num_perm=100, two_sided=True) - assert np.isclose(ISC, [0.9540602, 0.0267954]).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(), \ @@ -41,8 +42,8 @@ def test_ISFC(): assert D.shape == (4, 5, 2), "Loaded data has incorrect shape" np.random.seed(0) - (ISFC, p) = brainiak.isfc.isfc(D, return_p=True, - num_perm=100, two_sided=True) + (ISFC, p) = brainiak.isfc.isfc(D, return_p=True, num_perm=100, + two_sided=True, random_state=0) ground_truth = \ [[1, 1, 0, -1], From 7556639bfbab02fc7f0a3dc41cc7f6bec49e15e8 Mon Sep 17 00:00:00 2001 From: Christopher Baldassano Date: Wed, 13 Dec 2017 10:54:00 -0500 Subject: [PATCH 5/5] Moving p_from_null to utils and adding tests --- brainiak/isfc.py | 45 +---------------------------- brainiak/utils/utils.py | 47 +++++++++++++++++++++++++++++++ tests/isfc/test_isfc.py | 2 +- tests/utils/test_utils.py | 59 +++++++++++++++++++++++++++++++++++++++ 4 files changed, 108 insertions(+), 45 deletions(-) diff --git a/brainiak/isfc.py b/brainiak/isfc.py index 820bca7b8..32bb96f7d 100644 --- a/brainiak/isfc.py +++ b/brainiak/isfc.py @@ -32,7 +32,7 @@ from brainiak.fcma.util import compute_correlation import numpy as np from scipy import stats -from .utils.utils import phase_randomize, ecdf +from .utils.utils import phase_randomize, p_from_null def isc(D, collapse_subj=True, return_p=False, num_perm=1000, @@ -187,46 +187,3 @@ def isfc(D, collapse_subj=True, return_p=False, return ISFC[..., 0], p else: return ISFC[..., 0] - - -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). - - 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/brainiak/utils/utils.py b/brainiak/utils/utils.py index ee5032c52..4b60497e7 100644 --- a/brainiak/utils/utils.py +++ b/brainiak/utils/utils.py @@ -746,3 +746,50 @@ 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/tests/isfc/test_isfc.py b/tests/isfc/test_isfc.py index 7eb4b5391..e4d516be0 100644 --- a/tests/isfc/test_isfc.py +++ b/tests/isfc/test_isfc.py @@ -20,6 +20,7 @@ def test_ISC(): (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" @@ -41,7 +42,6 @@ def test_ISFC(): assert D.shape == (4, 5, 2), "Loaded data has incorrect shape" - np.random.seed(0) (ISFC, p) = brainiak.isfc.isfc(D, return_p=True, num_perm=100, two_sided=True, random_state=0) 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"