diff --git a/brainiak/isc.py b/brainiak/isc.py index 27d14a223..15201c241 100644 --- a/brainiak/isc.py +++ b/brainiak/isc.py @@ -486,7 +486,7 @@ def compute_summary_statistic(iscs, summary_statistic='mean', axis=None): The implementation is based on the work in [SilverDunlap1987]_. - .. [SilverDunlap1987] "Averaging corrlelation coefficients: should + .. [SilverDunlap1987] "Averaging correlation coefficients: should Fisher's z transformation be used?", N. C. Silver, W. P. Dunlap, 1987, Journal of Applied Psychology, 72, 146-148. https://doi.org/10.1037/0021-9010.72.1.146 @@ -643,7 +643,8 @@ def _threshold_nans(data, tolerate_nans): def bootstrap_isc(iscs, pairwise=False, summary_statistic='median', - n_bootstraps=1000, ci_percentile=95, random_state=None): + n_bootstraps=1000, ci_percentile=95, side='right', + random_state=None): """One-sample group-level bootstrap hypothesis test for ISCs @@ -653,16 +654,19 @@ def bootstrap_isc(iscs, pairwise=False, summary_statistic='median', should be either N ISC values for N subjects in the leave-one-out appraoch (pairwise=False), N(N-1)/2 ISC values for N subjects in the pairwise approach (pairwise=True). In the pairwise approach, ISC values should - correspond to the vectorized upper triangle of a square corrlation matrix + correspond to the vectorized upper triangle of a square correlation matrix (see scipy.stats.distance.squareform). Shifts bootstrap distribution by - actual summary statistic (effectively to zero) for two-tailed null - hypothesis test (Hall & Wilson, 1991). Uses subject-wise (not pair-wise) - resampling in the pairwise approach. Returns the observed ISC, the - confidence interval, and a p-value for the bootstrap hypothesis test, as - well as the bootstrap distribution of summary statistics. According to - Chen et al., 2016, this is the preferred nonparametric approach for - controlling false positive rates (FPR) for one-sample tests in the pairwise - approach. + actual summary statistic (effectively to zero) for null hypothesis test + (Hall & Wilson, 1991). Uses subject-wise (not pair-wise) resampling in the + pairwise approach. Returns the observed ISC, the confidence interval, and + a p-value for the bootstrap hypothesis test, as well as the bootstrap + distribution of summary statistics. The p-value corresponds to either a + 'two-sided', 'left'-, or 'right'-sided (default) test, as specified by + side. According to Chen et al., 2016, this is the preferred nonparametric + approach for controlling false positive rates (FPR) for one-sample tests + in the pairwise approach. The efficacy of this approach for controlling + FPRs in the leave-one-out approach has not yet been systematically + evaluated. The implementation is based on the work in [Chen2016]_ and [HallWilson1991]_. @@ -688,6 +692,9 @@ def bootstrap_isc(iscs, pairwise=False, summary_statistic='median', ci_percentile : int, default: 95 Percentile for computing confidence intervals + side : str + Perform one-sided ('left' or 'right') or 'two-sided' test + random_state = int or None, default: None Initial random seed @@ -790,7 +797,7 @@ def bootstrap_isc(iscs, pairwise=False, summary_statistic='median', # Get p-value for actual median from shifted distribution p = p_from_null(observed, shifted, - side='two-sided', exact=False, + side=side, exact=False, axis=0) return observed, ci, p, distribution @@ -1049,34 +1056,39 @@ def _permute_two_sample_iscs(iscs, group_parameters, i, pairwise=False, def permutation_isc(iscs, group_assignment=None, pairwise=False, # noqa: C901 summary_statistic='median', n_permutations=1000, - random_state=None): + side='right', random_state=None): """Group-level permutation test for ISCs For ISCs from one or more voxels or ROIs, permute group assignments to construct a permutation distribution. Input is a list or ndarray of - ISCs for a single voxel/ROI, or an ISCs-by-voxels ndarray. If two groups, - ISC values should stacked along first dimension (vertically), and a - group_assignment list (or 1d array) of same length as the number of - subjects should be provided to indicate groups. If no group_assignment - is provided, one-sample test is performed using a sign-flipping procedure. - Performs exact test if number of possible permutations (2**N for one-sample - sign-flipping, N! for two-sample shuffling) is less than or equal to number - of requested permutation; otherwise, performs approximate permutation test - using Monte Carlo resampling. ISC values should either be N ISC values for - N subjects in the leave-one-out approach (pairwise=False) or N(N-1)/2 ISC - values for N subjects in the pairwise approach (pairwise=True). In the - pairwise approach, ISC values should correspond to the vectorized upper - triangle of a square corrlation matrix (scipy.stats.distance.squareform). - Note that in the pairwise approach, group_assignment order should match the - row/column order of the subject-by-subject square ISC matrix even though - the input ISCs should be supplied as the vectorized upper triangle of the - square ISC matrix. Returns the observed ISC and permutation-based p-value - (two-tailed test), as well as the permutation distribution of summary - statistic. According to Chen et al., 2016, this is the preferred - nonparametric approach for controlling false positive rates (FPR) for - two-sample tests. This approach may yield inflated FPRs for one-sample - tests. + ISCs for a single voxel/ROI, or an ISCs-by-voxels ndarray. In the + leave-one-out approach, ISC values for two groups should be stacked + along first dimension (vertically) and a group_assignment list (or 1d + array) of same length as the number of subjects should be provided to + indicate groups. In the pairwise approach, pairwise ISCs should be + computed the across both groups at once; i.e. the pairwise ISC matrix + should be shaped N x N where N is the total number of subjects across + both groups, and should contain between-group ISC pairs. Pairwise ISC + input should correspond to the vectorized upper triangle of the square + pairwise ISC correlation matrix containing both groups. In the pairwise + approach, group_assignment order should match the row/column order of the + subject-by-subject square ISC matrix even though the input ISCs should be + supplied as the vectorized upper triangle of the square ISC matrix. If no + group_assignment is provided, one-sample test is performed using a sign- + flipping procedure. Performs exact test if number of possible permutations + (2**N for one-sample sign-flipping, N! for two-sample shuffling) is less + than or equal to number of requested permutation; otherwise, performs + approximate permutation test using Monte Carlo resampling. ISC values + should either be N ISC values for N subjects in the leave-one-out approach + (pairwise=False) or N(N-1)/2 ISC values for N subjects in the pairwise + approach (pairwise=True). Returns the observed ISC and permutation-based + p-value as well as the permutation distribution of summary statistic. + The p-value corresponds to either a 'two-sided', 'left'-, or 'right'-sided + (default) test, as specified by side. According to Chen et al., 2016, + this is the preferred nonparametric approach for controlling false + positive rates (FPR) for two-sample tests. This approach may yield + inflated FPRs for one-sample tests. The implementation is based on the work in [Chen2016]_. @@ -1097,6 +1109,9 @@ def permutation_isc(iscs, group_assignment=None, pairwise=False, # noqa: C901 n_permutations : int, default: 1000 Number of permutation iteration (randomizing group assignment) + side : str + Perform one-sided ('left' or 'right') or 'two-sided' test + random_state = int, None, or np.random.RandomState, default: None Initial random seed @@ -1224,18 +1239,19 @@ def permutation_isc(iscs, group_assignment=None, pairwise=False, # noqa: C901 # Get p-value for actual median from shifted distribution if exact_permutations: p = p_from_null(observed, distribution, - side='two-sided', exact=True, + side=side, exact=True, axis=0) else: p = p_from_null(observed, distribution, - side='two-sided', exact=False, + side=side, exact=False, axis=0) return observed, p, distribution def timeshift_isc(data, pairwise=False, summary_statistic='median', - n_shifts=1000, tolerate_nans=True, random_state=None): + n_shifts=1000, side='right', tolerate_nans=True, + random_state=None): """Circular time-shift randomization for one-sample ISC test @@ -1261,8 +1277,10 @@ def timeshift_isc(data, pairwise=False, summary_statistic='median', will not affect the pairwise approach; however, if a threshold float is provided, voxels that do not reach this threshold will be excluded. Note that accommodating NaNs may be notably slower than setting tolerate_nans to - False. Returns the observed ISC and p-values (two-tailed test), as well as - the null distribution of ISCs computed on randomly time-shifted data. + False. Returns the observed ISC and p-values, as well as the null + distribution of ISCs computed on randomly time-shifted data. The p-value + corresponds to either a 'two-sided', 'left'-, or 'right'-sided (default) + test, as specified by side. The implementation is based on the work in [Kauppi2010]_ and [Kauppi2014]_. @@ -1287,6 +1305,9 @@ def timeshift_isc(data, pairwise=False, summary_statistic='median', n_shifts : int, default: 1000 Number of randomly shifted samples + side : str + Perform one-sided ('left' or 'right') or 'two-sided' test + tolerate_nans : bool or float, default: True Accommodate NaNs (when averaging in leave-one-out approach) @@ -1377,14 +1398,15 @@ def timeshift_isc(data, pairwise=False, summary_statistic='median', # Get p-value for actual median from shifted distribution p = p_from_null(observed, distribution, - side='two-sided', exact=False, + side=side, exact=False, axis=0) return observed, p, distribution def phaseshift_isc(data, pairwise=False, summary_statistic='median', - n_shifts=1000, tolerate_nans=True, random_state=None): + n_shifts=1000, side='right', tolerate_nans=True, + random_state=None): """Phase randomization for one-sample ISC test @@ -1409,9 +1431,11 @@ def phaseshift_isc(data, pairwise=False, summary_statistic='median', N-1 subjects will be assigned NaN. Setting tolerate_nans to True or False will not affect the pairwise approach; however, if a threshold float is provided, voxels that do not reach this threshold will be excluded. Note - that accommodating NaNs may be notably slower than setting tolerate_nans to - False. Returns the observed ISC and p-values (two-tailed test), as well as - the null distribution of ISCs computed on phase-randomized data. + that accommodating NaNs may be notably slower than setting tolerate_nans + to False. Returns the observed ISC and p-values, as well as the null + distribution of ISCs computed on phase-randomized data. The p-value + corresponds to either a 'two-sided', 'left'-, or 'right'-sided (default) + test, as specified by side. The implementation is based on the work in [Lerner2011]_ and [Simony2016]_. @@ -1435,6 +1459,9 @@ def phaseshift_isc(data, pairwise=False, summary_statistic='median', n_shifts : int, default: 1000 Number of randomly shifted samples + side : str + Perform one-sided ('left' or 'right') or 'two-sided' test + tolerate_nans : bool or float, default: True Accommodate NaNs (when averaging in leave-one-out approach) @@ -1513,7 +1540,7 @@ def phaseshift_isc(data, pairwise=False, summary_statistic='median', # Get p-value for actual median from shifted distribution p = p_from_null(observed, distribution, - side='two-sided', exact=False, + side=side, exact=False, axis=0) return observed, p, distribution diff --git a/brainiak/utils/utils.py b/brainiak/utils/utils.py index 95f8ae736..08402d779 100644 --- a/brainiak/utils/utils.py +++ b/brainiak/utils/utils.py @@ -828,10 +828,10 @@ def p_from_null(observed, distribution, distribution : ndarray Null distribution of test statistic - side : str, default:'two-sided' + side : str, default: 'two-sided' Perform one-sided ('left' or 'right') or 'two-sided' test - axis: None or int, default:None + axis: None or int, default: None Axis indicating resampling iterations in input distribution Returns diff --git a/tests/isc/test_isc.py b/tests/isc/test_isc.py index fc4ea17c1..0f0256b81 100644 --- a/tests/isc/test_isc.py +++ b/tests/isc/test_isc.py @@ -45,10 +45,10 @@ def correlated_timeseries(n_subjects, n_TRs, noise=0, uncorrelated[:, 0, :].T), checks=False)) unc_mean = np.mean(squareform(np.corrcoef( uncorrelated[:, 0, :].T), checks=False)) - if unc_max < .3 and np.abs(unc_mean) < .001: + if unc_max < .25 and np.abs(unc_mean) < .001: correlated = False data = np.repeat(np.column_stack((signal, signal))[..., np.newaxis], - 20, axis=2) + n_subjects, axis=2) data = np.concatenate((data, uncorrelated), axis=1) data = data + np.random.randn(n_TRs, 3, n_subjects) * noise return data @@ -301,6 +301,7 @@ def test_bootstrap_isc(): assert np.all(iscs[:, -1] < .5) assert p[0] < .05 and p[1] < .05 assert p[2] > .01 + print(p) iscs = isc(data, pairwise=True) observed, ci, p, distribution = bootstrap_isc(iscs, pairwise=True)