diff --git a/NEWS.rst b/NEWS.rst index 3af1dcb58..1cc6725e1 100644 --- a/NEWS.rst +++ b/NEWS.rst @@ -17,6 +17,8 @@ Features - fmrisim: Added fitting procedure for SFNR and SNR parameters. Updated AR to be ARMA, involving both the generation and estimation. (`#372 `_) +- isc: Revamped ISC/ISFC functionality with more statistical tests. (`#384 + `_) - reprsimil: Added an option in BRSA to set the prior of SNR to be "equal" across all voxels. (`#387 `_) diff --git a/brainiak/isc.py b/brainiak/isc.py index 1678d92ea..b949d0dc5 100644 --- a/brainiak/isc.py +++ b/brainiak/isc.py @@ -51,10 +51,10 @@ import logging from scipy.spatial.distance import squareform from scipy.stats import pearsonr -from scipy.fftpack import fft, ifft import itertools as it from brainiak.fcma.util import compute_correlation -from brainiak.utils.utils import compute_p_from_null_distribution +from brainiak.utils.utils import (phase_randomize, p_from_null, + _check_timeseries_input) logger = logging.getLogger(__name__) @@ -65,6 +65,7 @@ "isc", "permutation_isc", "phaseshift_isc", + "squareform_isfc", "timeshift_isc", ] @@ -72,7 +73,7 @@ MAX_RANDOM_SEED = 2**32 - 1 -def isc(data, pairwise=False, summary_statistic=None): +def isc(data, pairwise=False, summary_statistic=None, tolerate_nans=True): """Intersubject correlation For each voxel or ROI, compute the Pearson correlation between each @@ -92,8 +93,22 @@ def isc(data, pairwise=False, summary_statistic=None): the last dimension is assumed to correspond to subjects. If only two subjects are supplied, simply compute Pearson correlation (precludes averaging in leave-one-out approach, and does not apply summary statistic). - Output is an ndarray where the first dimension is the number of subjects - or pairs and the second dimension is the number of voxels (or ROIs). + When using leave-one-out approach, NaNs are ignored when computing mean + time series of N-1 subjects (default: tolerate_nans=True). Alternatively, + you may supply a float between 0 and 1 indicating a threshold proportion + of N subjects with non-NaN values required when computing the average time + series for a given voxel. For example, if tolerate_nans=.8, ISCs will be + computed for any voxel where >= 80% of subjects have non-NaN values, + while voxels with < 80% non-NaN values will be assigned NaNs. If set to + False, NaNs are not tolerated and voxels with one or more NaNs among the + 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. Output is an ndarray where the first dimension is the number of + subjects or pairs and the second dimension is the number of voxels (or + ROIs). If only two subjects are supplied or a summary statistic is invoked, + the output is a ndarray n_voxels long. The implementation is based on the work in [Hasson2004]_. @@ -102,12 +117,15 @@ def isc(data, pairwise=False, summary_statistic=None): data : list or ndarray (n_TRs x n_voxels x n_subjects) fMRI data for which to compute ISC - pairwise : bool, default:False + pairwise : bool, default: False Whether to use pairwise (True) or leave-one-out (False) approach - summary_statistic : None or str, default:None + summary_statistic : None or str, default: None Return all ISCs or collapse using 'mean' or 'median' + tolerate_nans : bool or float, default: True + Accommodate NaNs (when averaging in leave-one-out approach) + Returns ------- iscs : subjects or pairs by voxels ndarray @@ -123,9 +141,16 @@ def isc(data, pairwise=False, summary_statistic=None): logger.info("Only two subjects! Simply computing Pearson correlation.") summary_statistic = None + # Check tolerate_nans input and use either mean/nanmean and exclude voxels + if tolerate_nans: + mean = np.nanmean + else: + mean = np.mean + data, mask = _threshold_nans(data, tolerate_nans) + # Loop over each voxel or ROI voxel_iscs = [] - for v in np.arange(n_voxels): + for v in np.arange(data.shape[1]): voxel_data = data[:, v, :].T if n_subjects == 2: iscs = pearsonr(voxel_data[0, :], voxel_data[1, :])[0] @@ -133,12 +158,16 @@ def isc(data, pairwise=False, summary_statistic=None): iscs = squareform(np.corrcoef(voxel_data), checks=False) elif not pairwise: iscs = np.array([pearsonr(subject, - np.mean(np.delete(voxel_data, - s, axis=0), - axis=0))[0] + mean(np.delete(voxel_data, + s, axis=0), + axis=0))[0] for s, subject in enumerate(voxel_data)]) voxel_iscs.append(iscs) - iscs = np.column_stack(voxel_iscs) + iscs_stack = np.column_stack(voxel_iscs) + + # Get ISCs back into correct shape after masking out NaNs + iscs = np.full((iscs_stack.shape[0], n_voxels), np.nan) + iscs[:, np.where(mask)[0]] = iscs_stack # Summarize results (if requested) if summary_statistic: @@ -146,22 +175,31 @@ def isc(data, pairwise=False, summary_statistic=None): summary_statistic=summary_statistic, axis=0)[np.newaxis, :] + # Throw away first dimension if singleton + if iscs.shape[0] == 1: + iscs = iscs[0] + return iscs -def isfc(data, pairwise=False, summary_statistic=None): +def isfc(data, targets=None, pairwise=False, summary_statistic=None, + vectorize_isfcs=True, tolerate_nans=True): """Intersubject functional correlation (ISFC) - For each voxel or ROI, compute the Pearson correlation between each - subject's response time series and other subjects' response time series - for all voxels or ROIs. If pairwise is False (default), use the - leave-one-out approach, where correlation is computed between each - subject and the average of the other subjects. If pairwise is True, - compute correlations between all pairs of subjects. If summary_statistic - is None, return N ISFC values for N subjects (leave-one-out) or N(N-1)/2 - ISFC values for each pair of N subjects, corresponding to the upper - triangle of the correlation matrix (see scipy.spatial.distance.squareform). + For each input voxel or ROI, compute the Pearson correlation between each + subject's response time series and all input voxels or ROIs in other + subjects. If a targets array is provided, instead compute ISFCs between + each input voxel time series and each voxel time series in targets across + subjects (resulting in asymmetric ISFC values). The targets array must have + the same number TRs and subjects as the input data. If pairwise is False + (default), use the leave-one-out approach, where correlation is computed + between each subject and the average of the other subjects. If pairwise is + True, compute correlations between all pairs of subjects. If a targets + array is provided, only the leave-one-out approach is supported. If + summary_statistic is None, return N ISFC values for N subjects (leave-one- + out) or N(N-1)/2 ISFC values for each pair of N subjects, corresponding to + the triangle of the correlation matrix (scipy.spatial.distance.squareform). Alternatively, use either 'mean' or 'median' to compute summary statistic of ISFCs (Fisher Z is applied if using mean). Input should be n_TRs by n_voxels by n_subjects array (e.g., brainiak.image.MaskedMultiSubjectData) @@ -170,9 +208,28 @@ def isfc(data, pairwise=False, summary_statistic=None): the last dimension is assumed to correspond to subjects. If only two subjects are supplied, simply compute ISFC between these two subjects (precludes averaging in leave-one-out approach, and does not apply summary - statistic). Output is n_voxels by n_voxels array if summary_statistic is - supplied; otherwise output is n_voxels by n_voxels by n_subjects (or - n_pairs) array. + statistic). Returns vectorized upper triangle of ISFC matrices for each + subject or pair when vectorized_isfcs=True, or full (redundant) 2D ISFC + matrices when vectorized_isfcs=False. When using leave-one-out approach, + NaNs are ignored when computing mean time series of N-1 subjects (default: + tolerate_nans=True). Alternatively, you may supply a float between 0 and + 1 indicating a threshold proportion of N subjects with non-NaN values + required when computing the average time series for a given voxel. For + example, if tolerate_nans=.8, ISCs will be computed for any voxel where + >= 80% of subjects have non-NaN values, while voxels with < 80% non-NaN + values will be assigned NaNs. If set to False, NaNs are not tolerated + and voxels with one or more NaNs among the 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. Output is either + a tuple comprising condensed off-diagonal ISFC values and the diagonal + ISC values if vectorize_isfcs=True, or a single ndarray with shape + n_subjects (or n_pairs) by n_voxels by n_voxels 3D array if + vectorize_isfcs=False (see brainiak.isc.squareform_isfc). If targets array + is provided (yielding asymmetric ISFCs), output ISFCs are not vectorized, + resulting in an n_subjects by n_voxels by n_targets ISFC array. If + summary_statistic is supplied, output is collapsed along first dimension. The implementation is based on the work in [Simony2016]_. @@ -181,133 +238,111 @@ def isfc(data, pairwise=False, summary_statistic=None): data : list or ndarray (n_TRs x n_voxels x n_subjects) fMRI data for which to compute ISFC + targets : list or ndarray (n_TRs x n_voxels x n_subjects), optional + fMRI data to use as targets for ISFC + pairwise : bool, default: False Whether to use pairwise (True) or leave-one-out (False) approach - summary_statistic : None or str, default:None + summary_statistic : None or str, default: None Return all ISFCs or collapse using 'mean' or 'median' + vectorize_isfcs : bool, default: True + Return tuple of condensed ISFCs and ISCs (True) or square (redundant) + ISFCs (False) + + tolerate_nans : bool or float, default: True + Accommodate NaNs (when averaging in leave-one-out approach) + Returns ------- - isfcs : subjects or pairs by voxels ndarray - ISFC for each subject or pair (or summary statistic) per voxel + isfcs : ndarray or tuple of ndarrays + ISFCs for each subject or pair (or summary statistic) per voxel pair """ # Check response time series input format data, n_TRs, n_voxels, n_subjects = _check_timeseries_input(data) - # Handle just two subjects properly - if n_subjects == 2: + # Check for optional targets input array + targets, t_n_TRs, t_n_voxels, t_n_subejcts, symmetric = ( + _check_targets_input(targets, data)) + if not symmetric: + pairwise = False + + # Check tolerate_nans input and use either mean/nanmean and exclude voxels + if tolerate_nans: + mean = np.nanmean + else: + mean = np.mean + data, mask = _threshold_nans(data, tolerate_nans) + targets, targets_mask = _threshold_nans(targets, tolerate_nans) + + # Handle just two subjects properly (for symmetric approach) + if symmetric and n_subjects == 2: isfcs = compute_correlation(np.ascontiguousarray(data[..., 0].T), - np.ascontiguousarray(data[..., 1].T)) + np.ascontiguousarray(data[..., 1].T), + return_nans=True) isfcs = (isfcs + isfcs.T) / 2 - assert isfcs.shape == (n_voxels, n_voxels) + isfcs = isfcs[..., np.newaxis] summary_statistic = None logger.info("Only two subjects! Computing ISFC between them.") - # Compute all pairwise ISFCs + # Compute all pairwise ISFCs (only for symmetric approach) elif pairwise: isfcs = [] for pair in it.combinations(np.arange(n_subjects), 2): isfc_pair = compute_correlation(np.ascontiguousarray( data[..., pair[0]].T), np.ascontiguousarray( - data[..., pair[1]].T)) - isfc_pair = (isfc_pair + isfc_pair.T) / 2 + targets[..., pair[1]].T), + return_nans=True) + if symmetric: + isfc_pair = (isfc_pair + isfc_pair.T) / 2 isfcs.append(isfc_pair) isfcs = np.dstack(isfcs) - assert isfcs.shape == (n_voxels, n_voxels, - n_subjects * (n_subjects - 1) / 2) # Compute ISFCs using leave-one-out approach - elif not pairwise: + else: # Roll subject axis for loop data = np.rollaxis(data, 2, 0) + targets = np.rollaxis(targets, 2, 0) # Compute leave-one-out ISFCs isfcs = [compute_correlation(np.ascontiguousarray(subject.T), - np.ascontiguousarray(np.mean( - np.delete(data, s, axis=0), - axis=0).T)) + np.ascontiguousarray(mean( + np.delete(targets, s, axis=0), + axis=0).T), + return_nans=True) for s, subject in enumerate(data)] # Transpose and average ISFC matrices for both directions - isfcs = np.dstack([(isfc_matrix + isfc_matrix.T) / 2 - for isfc_matrix in isfcs]) - assert isfcs.shape == (n_voxels, n_voxels, n_subjects) + isfcs = np.dstack([(isfc_matrix + isfc_matrix.T) / 2 if + symmetric else isfc_matrix for + isfc_matrix in isfcs]) + + # Get ISCs back into correct shape after masking out NaNs + isfcs_all = np.full((n_voxels, t_n_voxels, isfcs.shape[2]), np.nan) + isfcs_all[np.ix_(np.where(mask)[0], np.where(targets_mask)[0])] = isfcs + isfcs = np.moveaxis(isfcs_all, 2, 0) # Summarize results (if requested) if summary_statistic: isfcs = compute_summary_statistic(isfcs, summary_statistic=summary_statistic, - axis=2) - - return isfcs - - -def _check_timeseries_input(data): - - """Checks response time series input data for ISC analysis - - Input data should be a n_TRs by n_voxels by n_subjects ndarray - (e.g., brainiak.image.MaskedMultiSubjectData) or a list where each - item is a n_TRs by n_voxels ndarray for a given subject. Multiple - input ndarrays must be the same shape. If a 2D array is supplied, - the last dimension is assumed to correspond to subjects. This - function is only intended to be used internally by other - functions in this module (e.g., isc, isfc). - - Parameters - ---------- - data : ndarray or list - Time series data - - Returns - ------- - iscs : ndarray - Array of ISC values + axis=0) - n_TRs : int - Number of time points (TRs) + # Throw away first dimension if singleton + if isfcs.shape[0] == 1: + isfcs = isfcs[0] - n_voxels : int - Number of voxels (or ROIs) - - n_subjects : int - Number of subjects - - """ - - # Convert list input to 3d and check shapes - if type(data) == list: - data_shape = data[0].shape - for i, d in enumerate(data): - if d.shape != data_shape: - raise ValueError("All ndarrays in input list " - "must be the same shape!") - if d.ndim == 1: - data[i] = d[:, np.newaxis] - data = np.dstack(data) - - # Convert input ndarray to 3d and check shape - elif isinstance(data, np.ndarray): - if data.ndim == 2: - data = data[:, np.newaxis, :] - elif data.ndim == 3: - pass - else: - raise ValueError("Input ndarray should have 2 " - "or 3 dimensions (got {0})!".format(data.ndim)) - - # Infer subjects, TRs, voxels and log for user to check - n_TRs, n_voxels, n_subjects = data.shape - logger.info("Assuming {0} subjects with {1} time points " - "and {2} voxel(s) or ROI(s) for ISC analysis.".format( - n_subjects, n_TRs, n_voxels)) - - return data, n_TRs, n_voxels, n_subjects + # Optionally squareform to vectorize ISFC matrices (only if symmetric) + if vectorize_isfcs and symmetric: + isfcs, iscs = squareform_isfc(isfcs) + return isfcs, iscs + else: + return isfcs def _check_isc_input(iscs, pairwise=False): @@ -365,6 +400,60 @@ def _check_isc_input(iscs, pairwise=False): return iscs, n_subjects, n_voxels +def _check_targets_input(targets, data): + + """Checks ISFC targets input array + + For ISFC analysis, targets input array should either be a list + of n_TRs by n_targets arrays (where each array corresponds to + a subject), or an n_TRs by n_targets by n_subjects ndarray. This + function also checks the shape of the targets array against the + input data array. + + Parameters + ---------- + data : list or ndarray (n_TRs x n_voxels x n_subjects) + fMRI data for which to compute ISFC + + targets : list or ndarray (n_TRs x n_voxels x n_subjects) + fMRI data to use as targets for ISFC + + Returns + ------- + targets : ndarray (n_TRs x n_voxels x n_subjects) + ISFC targets with standadized structure + + n_TRs : int + Number of time points (TRs) for targets array + + n_voxels : int + Number of voxels (or ROIs) for targets array + + n_subjects : int + Number of subjects for targets array + + symmetric : bool + Indicator for symmetric vs. asymmetric + """ + + if isinstance(targets, np.ndarray) or isinstance(targets, list): + targets, n_TRs, n_voxels, n_subjects = ( + _check_timeseries_input(targets)) + if data.shape[0] != n_TRs: + raise ValueError("Targets array must have same number of " + "TRs as input data") + if data.shape[2] != n_subjects: + raise ValueError("Targets array must have same number of " + "subjects as input data") + symmetric = False + else: + targets = data + n_TRs, n_voxels, n_subjects = data.shape + symmetric = True + + return targets, n_TRs, n_voxels, n_subjects, symmetric + + def compute_summary_statistic(iscs, summary_statistic='mean', axis=None): """Computes summary statistics for ISCs @@ -385,7 +474,7 @@ def compute_summary_statistic(iscs, summary_statistic='mean', axis=None): iscs : list or ndarray ISC values - summary_statistic : str, default:'mean' + summary_statistic : str, default: 'mean' Summary statistic, 'mean' or 'median' axis : None or int or tuple of ints, optional @@ -411,6 +500,126 @@ def compute_summary_statistic(iscs, summary_statistic='mean', axis=None): return statistic +def squareform_isfc(isfcs, iscs=None): + + """Converts square ISFCs to condensed ISFCs (and ISCs), and vice-versa + + If input is a 2- or 3-dimensional array of square ISFC matrices, converts + this to the condensed off-diagonal ISFC values (i.e., the vectorized + triangle) and the diagonal ISC values. In this case, input must be a + single array of shape either n_voxels x n_voxels or n_subjects (or + n_pairs) x n_voxels x n_voxels. The condensed ISFC values are vectorized + according to scipy.spatial.distance.squareform, yielding n_voxels * + (n_voxels - 1) / 2 values comprising every voxel pair. Alternatively, if + input is an array of condensed off-diagonal ISFC values and an array of + diagonal ISC values, the square (redundant) ISFC values are returned. + This function mimics scipy.spatial.distance.squareform, but is intended + to retain the diagonal ISC values. + + Parameters + ---------- + isfcs : ndarray + Either condensed or redundant ISFC values + + iscs: ndarray, optional + Diagonal ISC values, required when input is condensed + + Returns + ------- + isfcs : ndarray or tuple of ndarrays + If condensed ISFCs are passed, a single redundant ISFC array is + returned; if redundant ISFCs are passed, both a condensed off- + diagonal ISFC array and the diagonal ISC values are returned + """ + + # Check if incoming ISFCs are square (redundant) + if not type(iscs) == np.ndarray and isfcs.shape[-2] == isfcs.shape[-1]: + if isfcs.ndim == 2: + isfcs = isfcs[np.newaxis, ...] + if isfcs.ndim == 3: + iscs = np.diagonal(isfcs, axis1=1, axis2=2) + isfcs = np.vstack([squareform(isfc, checks=False)[np.newaxis, :] + for isfc in isfcs]) + else: + raise ValueError("Square (redundant) ISFCs must be square " + "with multiple subjects or pairs of subjects " + "indexed by the first dimension") + if isfcs.shape[0] == iscs.shape[0] == 1: + isfcs, iscs = isfcs[0], iscs[0] + return isfcs, iscs + + # Otherwise, convert from condensed to redundant + else: + if isfcs.ndim == iscs.ndim == 1: + isfcs, iscs = isfcs[np.newaxis, :], iscs[np.newaxis, :] + isfcs_stack = [] + for isfc, isc in zip(isfcs, iscs): + isfc_sq = squareform(isfc, checks=False) + np.fill_diagonal(isfc_sq, isc) + isfcs_stack.append(isfc_sq[np.newaxis, ...]) + isfcs = np.vstack(isfcs_stack) + if isfcs.shape[0] == 1: + isfcs = isfcs[0] + return isfcs + + +def _threshold_nans(data, tolerate_nans): + + """Thresholds data based on proportion of subjects with NaNs + + Takes in data and a threshold value (float between 0.0 and 1.0) determining + the permissible proportion of subjects with non-NaN values. For example, if + threshold=.8, any voxel where >= 80% of subjects have non-NaN values will + be left unchanged, while any voxel with < 80% non-NaN values will be + assigned all NaN values and included in the nan_mask output. Note that the + output data has not been masked and will be same shape as the input data, + but may have a different number of NaNs based on the threshold. + + Parameters + ---------- + data : ndarray (n_TRs x n_voxels x n_subjects) + fMRI time series data + + tolerate_nans : bool or float (0.0 <= threshold <= 1.0) + Proportion of subjects with non-NaN values required to keep voxel + + Returns + ------- + data : ndarray (n_TRs x n_voxels x n_subjects) + fMRI time series data with adjusted NaNs + + nan_mask : ndarray (n_voxels,) + Boolean mask array of voxels with too many NaNs based on threshold + + """ + + nans = np.all(np.any(np.isnan(data), axis=0), axis=1) + + # Check tolerate_nans input and use either mean/nanmean and exclude voxels + if tolerate_nans is True: + logger.info("ISC computation will tolerate all NaNs when averaging") + + elif type(tolerate_nans) is float: + if not 0.0 <= tolerate_nans <= 1.0: + raise ValueError("If threshold to tolerate NaNs is a float, " + "it must be between 0.0 and 1.0; got {0}".format( + tolerate_nans)) + nans += ~(np.sum(~np.any(np.isnan(data), axis=0), axis=1) >= + data.shape[-1] * tolerate_nans) + logger.info("ISC computation will tolerate voxels with at least " + "{0} non-NaN values: {1} voxels do not meet " + "threshold".format(tolerate_nans, + np.sum(nans))) + + else: + logger.info("ISC computation will not tolerate NaNs when averaging") + + mask = ~nans + data = data[:, mask, :] + + return data, mask + + def bootstrap_isc(iscs, pairwise=False, summary_statistic='median', n_bootstraps=1000, ci_percentile=95, random_state=None): @@ -445,19 +654,19 @@ def bootstrap_isc(iscs, pairwise=False, summary_statistic='median', iscs : list or ndarray, ISCs by voxels array ISC values for one or more voxels - pairwise : bool, default:False + pairwise : bool, default: False Indicator of pairwise or leave-one-out, should match ISCs structure - summary_statistic : str, default:'median' + summary_statistic : str, default: 'median' Summary statistic, either 'median' (default) or 'mean' - n_bootstraps : int, default:1000 + n_bootstraps : int, default: 1000 Number of bootstrap samples (subject-level with replacement) - ci_percentile : int, default:95 + ci_percentile : int, default: 95 Percentile for computing confidence intervals - random_state = int or None, default:None + random_state = int or None, default: None Initial random seed Returns @@ -486,7 +695,7 @@ def bootstrap_isc(iscs, pairwise=False, summary_statistic='median', # Compute summary statistic for observed ISCs observed = compute_summary_statistic(iscs, summary_statistic=summary_statistic, - axis=0)[np.newaxis, :] + axis=0) # Set up an empty list to build our bootstrap distribution distribution = [] @@ -548,7 +757,6 @@ def bootstrap_isc(iscs, pairwise=False, summary_statistic='median', # Convert distribution to numpy array distribution = np.array(distribution) - assert distribution.shape == (n_bootstraps, n_voxels) # Compute CIs of median from bootstrap distribution (default: 95%) ci = (np.percentile(distribution, (100 - ci_percentile)/2, axis=0), @@ -559,12 +767,9 @@ def bootstrap_isc(iscs, pairwise=False, summary_statistic='median', shifted = distribution - observed # Get p-value for actual median from shifted distribution - p = compute_p_from_null_distribution(observed, shifted, - side='two-sided', exact=False, - axis=0) - - # Reshape p-values to fit with data shape - p = p[np.newaxis, :] + p = p_from_null(observed, shifted, + side='two-sided', exact=False, + axis=0) return observed, ci, p, distribution @@ -685,16 +890,16 @@ def _permute_one_sample_iscs(iscs, group_parameters, i, pairwise=False, i : int Permutation iteration - pairwise : bool, default:False + pairwise : bool, default: False Indicator of pairwise or leave-one-out, should match ISCs variable - summary_statistic : str, default:'median' + summary_statistic : str, default: 'median' Summary statistic, either 'median' (default) or 'mean' exact_permutations : list List of permutations - prng = None or np.random.RandomState, default:None + prng = None or np.random.RandomState, default: None Initial random seed Returns @@ -753,16 +958,16 @@ def _permute_two_sample_iscs(iscs, group_parameters, i, pairwise=False, i : int Permutation iteration - pairwise : bool, default:False + pairwise : bool, default: False Indicator of pairwise or leave-one-out, should match ISCs variable - summary_statistic : str, default:'median' + summary_statistic : str, default: 'median' Summary statistic, either 'median' (default) or 'mean' exact_permutations : list List of permutations - prng = None or np.random.RandomState, default:None + prng = None or np.random.RandomState, default: None Initial random seed Indicator of pairwise or leave-one-out, should match ISCs variable @@ -861,16 +1066,16 @@ def permutation_isc(iscs, group_assignment=None, pairwise=False, # noqa: C901 group_assignment : list or ndarray, group labels Group labels matching order of ISC input - pairwise : bool, default:False + pairwise : bool, default: False Indicator of pairwise or leave-one-out, should match ISCs variable - summary_statistic : str, default:'median' + summary_statistic : str, default: 'median' Summary statistic, either 'median' (default) or 'mean' - n_permutations : int, default:1000 + n_permutations : int, default: 1000 Number of permutation iteration (randomizing group assignment) - random_state = int, None, or np.random.RandomState, default:None + random_state = int, None, or np.random.RandomState, default: None Initial random seed Returns @@ -949,7 +1154,7 @@ def permutation_isc(iscs, group_assignment=None, pairwise=False, # noqa: C901 group_parameters['group_labels'][1], :], summary_statistic=summary_statistic, axis=0)) - observed = np.array(observed)[np.newaxis, :] + observed = np.array(observed) # Set up an empty list to build our permutation distribution distribution = [] @@ -993,26 +1198,22 @@ def permutation_isc(iscs, group_assignment=None, pairwise=False, # noqa: C901 # Convert distribution to numpy array distribution = np.array(distribution) - assert distribution.shape == (n_permutations, n_voxels) # Get p-value for actual median from shifted distribution if exact_permutations: - p = compute_p_from_null_distribution(observed, distribution, - side='two-sided', exact=True, - axis=0) + p = p_from_null(observed, distribution, + side='two-sided', exact=True, + axis=0) else: - p = compute_p_from_null_distribution(observed, distribution, - side='two-sided', exact=False, - axis=0) - - # Reshape p-values to fit with data shape - p = p[np.newaxis, :] + p = p_from_null(observed, distribution, + side='two-sided', exact=False, + axis=0) return observed, p, distribution def timeshift_isc(data, pairwise=False, summary_statistic='median', - n_shifts=1000, random_state=None): + n_shifts=1000, tolerate_nans=True, random_state=None): """Circular time-shift randomization for one-sample ISC test @@ -1026,7 +1227,19 @@ def timeshift_isc(data, pairwise=False, summary_statistic='median', each item is a time-points by voxels ndarray for a given subject. Multiple input ndarrays must be the same shape. If a single ndarray is supplied, the last dimension is assumed to correspond to subjects. - Returns the observed ISC and p-values (two-tailed test), as well as + When using leave-one-out approach, NaNs are ignored when computing mean + time series of N-1 subjects (default: tolerate_nans=True). Alternatively, + you may supply a float between 0 and 1 indicating a threshold proportion + of N subjects with non-NaN values required when computing the average time + series for a given voxel. For example, if tolerate_nans=.8, ISCs will be + computed for any voxel where >= 80% of subjects have non-NaN values, + while voxels with < 80% non-NaN values will be assigned NaNs. If set to + False, NaNs are not tolerated and voxels with one or more NaNs among the + 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 randomly time-shifted data. The implementation is based on the work in [Kauppi2010]_ and @@ -1046,13 +1259,16 @@ def timeshift_isc(data, pairwise=False, summary_statistic='median', pairwise : bool, default: False Whether to use pairwise (True) or leave-one-out (False) approach - summary_statistic : str, default:'median' + summary_statistic : str, default: 'median' Summary statistic, either 'median' (default) or 'mean' - n_shifts : int, default:1000 + n_shifts : int, default: 1000 Number of randomly shifted samples - random_state = int, None, or np.random.RandomState, default:None + tolerate_nans : bool or float, default: True + Accommodate NaNs (when averaging in leave-one-out approach) + + random_state = int, None, or np.random.RandomState, default: None Initial random seed Returns @@ -1072,7 +1288,8 @@ def timeshift_isc(data, pairwise=False, summary_statistic='median', # Get actual observed ISC observed = isc(data, pairwise=pairwise, - summary_statistic=summary_statistic) + summary_statistic=summary_statistic, + tolerate_nans=tolerate_nans) # Roll axis to get subjects in first dimension for loop if pairwise: @@ -1105,7 +1322,8 @@ def timeshift_isc(data, pairwise=False, summary_statistic='median', # Compute null ISC on shifted data for pairwise approach shifted_isc = isc(shifted_data, pairwise=pairwise, - summary_statistic=summary_statistic) + summary_statistic=summary_statistic, + tolerate_nans=tolerate_nans) # In leave-one-out, apply shift only to each left-out participant elif not pairwise: @@ -1117,7 +1335,8 @@ def timeshift_isc(data, pairwise=False, summary_statistic='median', nonshifted_mean = np.mean(np.delete(data, s, 2), axis=2) loo_isc = isc(np.dstack((shifted_subject, nonshifted_mean)), pairwise=False, - summary_statistic=None) + summary_statistic=None, + tolerate_nans=tolerate_nans) shifted_isc.append(loo_isc) # Get summary statistics across left-out subjects @@ -1133,21 +1352,17 @@ def timeshift_isc(data, pairwise=False, summary_statistic='median', # Convert distribution to numpy array distribution = np.vstack(distribution) - assert distribution.shape == (n_shifts, n_voxels) # Get p-value for actual median from shifted distribution - p = compute_p_from_null_distribution(observed, distribution, - side='two-sided', exact=False, - axis=0) - - # Reshape p-values to fit with data shape - p = p[np.newaxis, :] + p = p_from_null(observed, distribution, + side='two-sided', exact=False, + axis=0) return observed, p, distribution def phaseshift_isc(data, pairwise=False, summary_statistic='median', - n_shifts=1000, random_state=None): + n_shifts=1000, tolerate_nans=True, random_state=None): """Phase randomization for one-sample ISC test @@ -1161,7 +1376,19 @@ def phaseshift_isc(data, pairwise=False, summary_statistic='median', each item is a time-points by voxels ndarray for a given subject. Multiple input ndarrays must be the same shape. If a single ndarray is supplied, the last dimension is assumed to correspond to subjects. - Returns the observed ISC and p-values (two-tailed test), as well as + When using leave-one-out approach, NaNs are ignored when computing mean + time series of N-1 subjects (default: tolerate_nans=True). Alternatively, + you may supply a float between 0 and 1 indicating a threshold proportion + of N subjects with non-NaN values required when computing the average time + series for a given voxel. For example, if tolerate_nans=.8, ISCs will be + computed for any voxel where >= 80% of subjects have non-NaN values, + while voxels with < 80% non-NaN values will be assigned NaNs. If set to + False, NaNs are not tolerated and voxels with one or more NaNs among the + 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. The implementation is based on the work in [Lerner2011]_ and @@ -1180,13 +1407,16 @@ def phaseshift_isc(data, pairwise=False, summary_statistic='median', pairwise : bool, default: False Whether to use pairwise (True) or leave-one-out (False) approach - summary_statistic : str, default:'median' + summary_statistic : str, default: 'median' Summary statistic, either 'median' (default) or 'mean' - n_shifts : int, default:1000 + n_shifts : int, default: 1000 Number of randomly shifted samples - random_state = int, None, or np.random.RandomState, default:None + tolerate_nans : bool or float, default: True + Accommodate NaNs (when averaging in leave-one-out approach) + + random_state = int, None, or np.random.RandomState, default: None Initial random seed Returns @@ -1206,7 +1436,8 @@ def phaseshift_isc(data, pairwise=False, summary_statistic='median', # Get actual observed ISC observed = isc(data, pairwise=pairwise, - summary_statistic=summary_statistic) + summary_statistic=summary_statistic, + tolerate_nans=tolerate_nans) # Iterate through randomized shifts to create null distribution distribution = [] @@ -1218,58 +1449,32 @@ def phaseshift_isc(data, pairwise=False, summary_statistic='median', else: prng = np.random.RandomState(random_state) - # Get randomized phase shifts - if n_TRs % 2 == 0: - # Why are we indexing from 1 not zero here? n_TRs / -1 long? - pos_freq = np.arange(1, data.shape[0] // 2) - neg_freq = np.arange(data.shape[0] - 1, data.shape[0] // 2, -1) - else: - pos_freq = np.arange(1, (data.shape[0] - 1) // 2 + 1) - neg_freq = np.arange(data.shape[0] - 1, - (data.shape[0] - 1) // 2, -1) - - phase_shifts = prng.rand(len(pos_freq), 1, n_subjects) * 2 * np.math.pi + # Get shifted version of data + shifted_data = phase_randomize(data, random_state=prng) # In pairwise approach, apply all shifts then compute pairwise ISCs if pairwise: - # Fast Fourier transform along time dimension of data - fft_data = fft(data, axis=0) - - # Shift pos and neg frequencies symmetrically, to keep signal real - fft_data[pos_freq, :, :] *= np.exp(1j * phase_shifts) - fft_data[neg_freq, :, :] *= np.exp(-1j * phase_shifts) - - # Inverse FFT to put data back in time domain for ISC - shifted_data = np.real(ifft(fft_data, axis=0)) - # Compute null ISC on shifted data for pairwise approach shifted_isc = isc(shifted_data, pairwise=True, - summary_statistic=summary_statistic) + summary_statistic=summary_statistic, + tolerate_nans=tolerate_nans) # In leave-one-out, apply shift only to each left-out participant elif not pairwise: - # Roll subject axis in phaseshifts for loop - phase_shifts = np.rollaxis(phase_shifts, 2, 0) + # Roll subject axis of phase-randomized data + shifted_data = np.rollaxis(shifted_data, 2, 0) shifted_isc = [] - for s, shift in enumerate(phase_shifts): - - # Apply FFT to left-out subject - fft_subject = fft(data[:, :, s], axis=0) - - # Shift pos and neg frequencies symmetrically, keep signal real - fft_subject[pos_freq, :] *= np.exp(1j * shift) - fft_subject[neg_freq, :] *= np.exp(-1j * shift) - - # Inverse FFT to put data back in time domain for ISC - shifted_subject = np.real(ifft(fft_subject, axis=0)) + for s, shifted_subject in enumerate(shifted_data): # ISC of shifted left-out subject vs mean of N-1 subjects - nonshifted_mean = np.mean(np.delete(data, s, 2), axis=2) + nonshifted_mean = np.mean(np.delete(data, s, axis=2), + axis=2) loo_isc = isc(np.dstack((shifted_subject, nonshifted_mean)), - pairwise=False, summary_statistic=None) + pairwise=False, summary_statistic=None, + tolerate_nans=tolerate_nans) shifted_isc.append(loo_isc) # Get summary statistics across left-out subjects @@ -1283,14 +1488,10 @@ def phaseshift_isc(data, pairwise=False, summary_statistic='median', # Convert distribution to numpy array distribution = np.vstack(distribution) - assert distribution.shape == (n_shifts, n_voxels) # Get p-value for actual median from shifted distribution - p = compute_p_from_null_distribution(observed, distribution, - side='two-sided', exact=False, - axis=0) - - # Reshape p-values to fit with data shape - p = p[np.newaxis, :] + p = p_from_null(observed, distribution, + side='two-sided', exact=False, + axis=0) return observed, p, distribution diff --git a/brainiak/utils/utils.py b/brainiak/utils/utils.py index c44ac5488..bdb119cc6 100644 --- a/brainiak/utils/utils.py +++ b/brainiak/utils/utils.py @@ -17,9 +17,7 @@ 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 import logging logger = logging.getLogger(__name__) @@ -31,10 +29,8 @@ __all__ = [ "center_mass_exp", - "compute_p_from_null_distribution", "concatenate_not_none", "cov2corr", - "ecdf", "from_tri_2_sym", "from_sym_2_tri", "gen_design", @@ -697,25 +693,32 @@ def usable_cpu_count(): return result -def phase_randomize(D, random_state=0): - """Randomly shift signal phases +def phase_randomize(data, voxelwise=False, random_state=None): + """Randomize phase of time series across subjects - 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. + For each subject, apply Fourier transform to voxel time series + and then randomly shift the phase of each frequency before inverting + back into the time domain. This yields time series with the same power + spectrum (and thus the same autocorrelation) as the original time series + but will remove any meaningful temporal relationships among time series + across subjects. By default (voxelwise=False), the same phase shift is + applied across all voxels; however if voxelwise=True, different random + phase shifts are applied to each voxel. The typical input is a time by + voxels by subjects ndarray. The first dimension is assumed to be the + time dimension and will be phase randomized. If a 2-dimensional ndarray + is provided, the last dimension is assumed to be subjects, and different + phase randomizations will be applied to each subject. - 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. + The implementation is based on the work in [Lerner2011]_ and + [Simony2016]_. Parameters ---------- - D : voxel by time by subject ndarray - fMRI data to be phase randomized + data : ndarray (n_TRs x n_voxels x n_subjects) + Data to be phase randomized (per subject) + + voxelwise : bool, default: False + Apply same (False) or different (True) randomizations across voxels random_state : RandomState or an int seed (0 by default) A random number generator instance to define the state of the @@ -723,123 +726,59 @@ def phase_randomize(D, random_state=0): Returns ---------- - ndarray of same shape as D - phase randomized timecourses + shifted_data : ndarray (n_TRs x n_voxels x n_subjects) + Phase-randomized time series """ - 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)) + # Check if input is 2-dimensional + data_ndim = data.ndim + # Get basic shape of data + data, n_TRs, n_voxels, n_subjects = _check_timeseries_input(data) -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 + # Random seed to be deterministically re-randomized at each iteration + if isinstance(random_state, np.random.RandomState): + prng = random_state + else: + prng = np.random.RandomState(random_state) - 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) + # Get randomized phase shifts + if n_TRs % 2 == 0: + # Why are we indexing from 1 not zero here? n_TRs / -1 long? + pos_freq = np.arange(1, data.shape[0] // 2) + neg_freq = np.arange(data.shape[0] - 1, data.shape[0] // 2, -1) + else: + pos_freq = np.arange(1, (data.shape[0] - 1) // 2 + 1) + neg_freq = np.arange(data.shape[0] - 1, + (data.shape[0] - 1) // 2, -1) - def ecdf_fun(q): - return yp[np.searchsorted(xp, q, side="right")] + if not voxelwise: + phase_shifts = (prng.rand(len(pos_freq), 1, n_subjects) + * 2 * np.math.pi) + else: + phase_shifts = (prng.rand(len(pos_freq), n_voxels, n_subjects) + * 2 * np.math.pi) - return ecdf_fun + # Fast Fourier transform along time dimension of data + fft_data = fft(data, axis=0) + # Shift pos and neg frequencies symmetrically, to keep signal real + fft_data[pos_freq, :, :] *= np.exp(1j * phase_shifts) + fft_data[neg_freq, :, :] *= np.exp(-1j * phase_shifts) -def p_from_null(X, two_sided=False, - max_null_input=None, min_null_input=None): - """Compute p value of true result from null distribution + # Inverse FFT to put data back in time domain + shifted_data = np.real(ifft(fft_data, axis=0)) - 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). + # Go back to 2-dimensions if input was 2-dimensional + if data_ndim == 2: + shifted_data = shifted_data[:, 0, :] - 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. + return shifted_data - 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:] - If max_null_input and min_null_input are provided, - X should contain only the real result - - 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) - - max_null_input : ndarray with num_perm (see `brainiak.isfc`) entries - By default this array is derived from the X input array, - which can be very large and takes up huge memory space. - To save memory, the function which calls p_from_null - should provide this array as input. - - min_null_input : ndarray with num_perm (see `brainiak.isfc`) entries - See max_null_input - - Returns - ------- - p : ndarray the same shape as X, without the last dimension - p values for each true X value under the null distribution - """ - if (min_null_input is None) or (max_null_input is None): - real_data = X[..., 0] - 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) - else: - real_data = X - # maximum & minimum in each null dataset should be provided as input - max_null = max_null_input - min_null = min_null_input - # 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(real_data), - min_null_ecdf(real_data)) - p = np.minimum(p, 1) - else: - p = 1 - max_null_ecdf(real_data) - - return p - - -def compute_p_from_null_distribution(observed, distribution, - side='two-sided', exact=False, - axis=None): +def p_from_null(observed, distribution, + side='two-sided', exact=False, + axis=None): """Compute p-value from null distribution Returns the p-value for an observed test statistic given a null @@ -888,13 +827,13 @@ def compute_p_from_null_distribution(observed, distribution, logger.info("Assuming {0} resampling iterations".format(n_samples)) if side == 'two-sided': - # numerator for two-sided test + # Numerator for two-sided test numerator = np.sum(np.abs(distribution) >= np.abs(observed), axis=axis) elif side == 'left': - # numerator for one-sided test in left tail + # Numerator for one-sided test in left tail numerator = np.sum(distribution <= observed, axis=axis) elif side == 'right': - # numerator for one-sided test in right tail + # Numerator for one-sided test in right tail numerator = np.sum(distribution >= observed, axis=axis) # If exact test all possible permutations and do not adjust @@ -907,3 +846,66 @@ def compute_p_from_null_distribution(observed, distribution, p = (numerator + 1) / (n_samples + 1) return p + + +def _check_timeseries_input(data): + + """Checks response time series input data (e.g., for ISC analysis) + + Input data should be a n_TRs by n_voxels by n_subjects ndarray + (e.g., brainiak.image.MaskedMultiSubjectData) or a list where each + item is a n_TRs by n_voxels ndarray for a given subject. Multiple + input ndarrays must be the same shape. If a 2D array is supplied, + the last dimension is assumed to correspond to subjects. This + function is generally intended to be used internally by other + functions module (e.g., isc, isfc in brainiak.isc). + + Parameters + ---------- + data : ndarray or list + Time series data + + Returns + ------- + data : ndarray + Input time series data with standardized structure + + n_TRs : int + Number of time points (TRs) + + n_voxels : int + Number of voxels (or ROIs) + + n_subjects : int + Number of subjects + + """ + + # Convert list input to 3d and check shapes + if type(data) == list: + data_shape = data[0].shape + for i, d in enumerate(data): + if d.shape != data_shape: + raise ValueError("All ndarrays in input list " + "must be the same shape!") + if d.ndim == 1: + data[i] = d[:, np.newaxis] + data = np.dstack(data) + + # Convert input ndarray to 3d and check shape + elif isinstance(data, np.ndarray): + if data.ndim == 2: + data = data[:, np.newaxis, :] + elif data.ndim == 3: + pass + else: + raise ValueError("Input ndarray should have 2 " + "or 3 dimensions (got {0})!".format(data.ndim)) + + # Infer subjects, TRs, voxels and log for user to check + n_TRs, n_voxels, n_subjects = data.shape + logger.info("Assuming {0} subjects with {1} time points " + "and {2} voxel(s) or ROI(s) for ISC analysis.".format( + n_subjects, n_TRs, n_voxels)) + + return data, n_TRs, n_voxels, n_subjects diff --git a/docs/newsfragments/396.feature b/docs/newsfragments/396.feature new file mode 100644 index 000000000..e6731e504 --- /dev/null +++ b/docs/newsfragments/396.feature @@ -0,0 +1 @@ +utils.phase_randomize outputs phase-scrambled input data, not tied to ISC/ISFC diff --git a/docs/newsfragments/397.feature b/docs/newsfragments/397.feature new file mode 100644 index 000000000..fabbc3014 --- /dev/null +++ b/docs/newsfragments/397.feature @@ -0,0 +1 @@ +utils.p_from_null no longer estimates distribution, simply returns p-value diff --git a/docs/newsfragments/398.feature b/docs/newsfragments/398.feature new file mode 100644 index 000000000..9935b6d99 --- /dev/null +++ b/docs/newsfragments/398.feature @@ -0,0 +1 @@ +ISC/ISFC analyses will now tolerate NaNs or tolerate a proportion of NaNs diff --git a/docs/newsfragments/399.feature b/docs/newsfragments/399.feature new file mode 100644 index 000000000..d51d1b372 --- /dev/null +++ b/docs/newsfragments/399.feature @@ -0,0 +1 @@ +ISFC will now output either vectorized triangle and diagonal of square matrices diff --git a/docs/newsfragments/405.feature b/docs/newsfragments/405.feature new file mode 100644 index 000000000..5bc21b1cd --- /dev/null +++ b/docs/newsfragments/405.feature @@ -0,0 +1 @@ +Updates NaN-toleration, ISFC output shape, and targets array for ISFC diff --git a/docs/newsfragments/409.feature b/docs/newsfragments/409.feature new file mode 100644 index 000000000..1382a649f --- /dev/null +++ b/docs/newsfragments/409.feature @@ -0,0 +1 @@ +Asymmetric ISFC analysis can now be performed based on a targets array diff --git a/tests/isc/mask.nii.gz b/tests/isc/mask.nii.gz deleted file mode 100644 index 83b31ffbd..000000000 Binary files a/tests/isc/mask.nii.gz and /dev/null differ diff --git a/tests/isc/subj1.nii.gz b/tests/isc/subj1.nii.gz deleted file mode 100644 index 9cf26af98..000000000 Binary files a/tests/isc/subj1.nii.gz and /dev/null differ diff --git a/tests/isc/subj2.nii.gz b/tests/isc/subj2.nii.gz deleted file mode 100644 index ec6bb9a44..000000000 Binary files a/tests/isc/subj2.nii.gz and /dev/null differ diff --git a/tests/isc/test_isc.py b/tests/isc/test_isc.py index 6bd681cc3..fc4ea17c1 100644 --- a/tests/isc/test_isc.py +++ b/tests/isc/test_isc.py @@ -1,7 +1,9 @@ import numpy as np import logging +import pytest from brainiak.isc import (isc, isfc, bootstrap_isc, permutation_isc, - timeshift_isc, phaseshift_isc) + squareform_isfc, timeshift_isc, + phaseshift_isc) from scipy.spatial.distance import squareform logger = logging.getLogger(__name__) @@ -114,20 +116,22 @@ def test_isc_options(): iscs_loo = isc(data, pairwise=False, summary_statistic=None) assert iscs_loo.shape == (n_subjects, n_voxels) + # Just two subjects + iscs_loo = isc(data[..., :2], pairwise=False, summary_statistic=None) + assert iscs_loo.shape == (n_voxels,) + iscs_pw = isc(data, pairwise=True, summary_statistic=None) assert iscs_pw.shape == (n_subjects*(n_subjects-1)/2, n_voxels) # Check summary statistics isc_mean = isc(data, pairwise=False, summary_statistic='mean') - assert isc_mean.shape == (1, n_voxels) + assert isc_mean.shape == (n_voxels,) isc_median = isc(data, pairwise=False, summary_statistic='median') - assert isc_median.shape == (1, n_voxels) + assert isc_median.shape == (n_voxels,) - try: + with pytest.raises(ValueError): isc(data, pairwise=False, summary_statistic='min') - except ValueError: - logger.info("Correctly caught unexpected summary statistic") logger.info("Finished testing ISC options") @@ -150,6 +154,97 @@ def test_isc_output(): logger.info("Finished testing ISC outputs") +# Check for proper handling of NaNs in ISC +def test_isc_nans(): + + # Set parameters for toy time series data + n_subjects = 20 + n_TRs = 60 + n_voxels = 30 + random_state = 42 + + logger.info("Testing ISC options") + + data = simulated_timeseries(n_subjects, n_TRs, + n_voxels=n_voxels, data_type='array', + random_state=random_state) + + # Inject NaNs into data + data[0, 0, 0] = np.nan + + # Don't tolerate NaNs, should lose zeroeth voxel + iscs_loo = isc(data, pairwise=False, tolerate_nans=False) + assert np.sum(np.isnan(iscs_loo)) == n_subjects + + # Tolerate all NaNs, only subject with NaNs yields NaN + iscs_loo = isc(data, pairwise=False, tolerate_nans=True) + assert np.sum(np.isnan(iscs_loo)) == 1 + + # Pairwise approach shouldn't care + iscs_pw_T = isc(data, pairwise=True, tolerate_nans=True) + iscs_pw_F = isc(data, pairwise=True, tolerate_nans=False) + assert np.allclose(iscs_pw_T, iscs_pw_F, equal_nan=True) + + assert (np.sum(np.isnan(iscs_pw_T)) == + np.sum(np.isnan(iscs_pw_F)) == + n_subjects - 1) + + # Set proportion of nans to reject (70% and 90% non-NaN) + data[0, 0, :] = np.nan + data[0, 1, :n_subjects - int(n_subjects * .7)] = np.nan + data[0, 2, :n_subjects - int(n_subjects * .9)] = np.nan + + iscs_loo_T = isc(data, pairwise=False, tolerate_nans=True) + iscs_loo_F = isc(data, pairwise=False, tolerate_nans=False) + iscs_loo_95 = isc(data, pairwise=False, tolerate_nans=.95) + iscs_loo_90 = isc(data, pairwise=False, tolerate_nans=.90) + iscs_loo_80 = isc(data, pairwise=False, tolerate_nans=.8) + iscs_loo_70 = isc(data, pairwise=False, tolerate_nans=.7) + iscs_loo_60 = isc(data, pairwise=False, tolerate_nans=.6) + + assert (np.sum(np.isnan(iscs_loo_F)) == + np.sum(np.isnan(iscs_loo_95)) == 60) + assert (np.sum(np.isnan(iscs_loo_80)) == + np.sum(np.isnan(iscs_loo_90)) == 42) + assert (np.sum(np.isnan(iscs_loo_T)) == + np.sum(np.isnan(iscs_loo_60)) == + np.sum(np.isnan(iscs_loo_70)) == 28) + assert np.array_equal(np.sum(np.isnan(iscs_loo_F), axis=0), + np.sum(np.isnan(iscs_loo_95), axis=0)) + assert np.array_equal(np.sum(np.isnan(iscs_loo_80), axis=0), + np.sum(np.isnan(iscs_loo_90), axis=0)) + assert np.all((np.array_equal( + np.sum(np.isnan(iscs_loo_T), axis=0), + np.sum(np.isnan(iscs_loo_60), axis=0)), + np.array_equal( + np.sum(np.isnan(iscs_loo_T), axis=0), + np.sum(np.isnan(iscs_loo_70), axis=0)), + np.array_equal( + np.sum(np.isnan(iscs_loo_60), axis=0), + np.sum(np.isnan(iscs_loo_70), axis=0)))) + + data = simulated_timeseries(n_subjects, n_TRs, + n_voxels=n_voxels, data_type='array', + random_state=random_state) + + # Make sure voxel with NaNs across all subjects is always removed + data[0, 0, :] = np.nan + iscs_loo_T = isc(data, pairwise=False, tolerate_nans=True) + iscs_loo_F = isc(data, pairwise=False, tolerate_nans=False) + assert np.allclose(iscs_loo_T, iscs_loo_F, equal_nan=True) + assert (np.sum(np.isnan(iscs_loo_T)) == + np.sum(np.isnan(iscs_loo_F)) == + n_subjects) + + iscs_pw_T = isc(data, pairwise=True, tolerate_nans=True) + iscs_pw_F = isc(data, pairwise=True, tolerate_nans=False) + assert np.allclose(iscs_pw_T, iscs_pw_F, equal_nan=True) + + assert (np.sum(np.isnan(iscs_pw_T)) == + np.sum(np.isnan(iscs_pw_F)) == + n_subjects * (n_subjects - 1) / 2) + + # Test one-sample bootstrap test def test_bootstrap_isc(): @@ -204,15 +299,15 @@ def test_bootstrap_isc(): observed, ci, p, distribution = bootstrap_isc(iscs, pairwise=False) assert np.all(iscs[:, :2] > .5) assert np.all(iscs[:, -1] < .5) - assert p[0, 0] < .05 and p[0, 1] < .05 - assert p[0, 2] > .01 + assert p[0] < .05 and p[1] < .05 + assert p[2] > .01 iscs = isc(data, pairwise=True) observed, ci, p, distribution = bootstrap_isc(iscs, pairwise=True) assert np.all(iscs[:, :2] > .5) assert np.all(iscs[:, -1] < .5) - assert p[0, 0] < .05 and p[0, 1] < .05 - assert p[0, 2] > .01 + assert p[0] < .05 and p[1] < .05 + assert p[2] > .01 # Check that ISC computation and bootstrap observed are same iscs = isc(data, pairwise=False) @@ -347,29 +442,31 @@ def test_permutation_isc(): observed, p, distribution = permutation_isc(iscs, pairwise=False) assert np.all(iscs[:, :2] > .5) assert np.all(iscs[:, -1] < .5) - assert p[0, 0] < .05 and p[0, 1] < .05 - assert p[0, 2] > .01 + assert p[0] < .05 and p[1] < .05 + assert p[2] > .01 iscs = isc(data, pairwise=True) observed, p, distribution = permutation_isc(iscs, pairwise=True) assert np.all(iscs[:, :2] > .5) assert np.all(iscs[:, -1] < .5) - assert p[0, 0] < .05 and p[0, 1] < .05 - assert p[0, 2] > .01 + assert p[0] < .05 and p[1] < .05 + assert p[2] > .01 # Check that ISC computation and permutation observed are same iscs = isc(data, pairwise=False) observed, p, distribution = permutation_isc(iscs, pairwise=False, summary_statistic='median') - assert np.array_equal(observed, isc(data, pairwise=False, - summary_statistic='median')) + assert np.allclose(observed, isc(data, pairwise=False, + summary_statistic='median'), + rtol=1e-03) # Check that ISC computation and permuation observed are same iscs = isc(data, pairwise=True) observed, p, distribution = permutation_isc(iscs, pairwise=True, summary_statistic='mean') - assert np.array_equal(observed, isc(data, pairwise=True, - summary_statistic='mean')) + assert np.allclose(observed, isc(data, pairwise=True, + summary_statistic='mean'), + rtol=1e-03) logger.info("Finished testing permutaton test") @@ -410,29 +507,31 @@ def test_timeshift_isc(): observed, p, distribution = timeshift_isc(data, pairwise=False) assert np.all(iscs[:, :2] > .5) assert np.all(iscs[:, -1] < .5) - assert p[0, 0] < .05 and p[0, 1] < .05 - assert p[0, 2] > .01 + assert p[0] < .05 and p[1] < .05 + assert p[2] > .01 iscs = isc(data, pairwise=True) observed, p, distribution = timeshift_isc(data, pairwise=True) assert np.all(iscs[:, :2] > .5) assert np.all(iscs[:, -1] < .5) - assert p[0, 0] < .05 and p[0, 1] < .05 - assert p[0, 2] > .01 + assert p[0] < .05 and p[1] < .05 + assert p[2] > .01 # Check that ISC computation and permutation observed are same iscs = isc(data, pairwise=False) observed, p, distribution = timeshift_isc(data, pairwise=False, summary_statistic='median') - assert np.array_equal(observed, isc(data, pairwise=False, - summary_statistic='median')) + assert np.allclose(observed, isc(data, pairwise=False, + summary_statistic='median'), + rtol=1e-03) # Check that ISC computation and permuation observed are same iscs = isc(data, pairwise=True) observed, p, distribution = timeshift_isc(data, pairwise=True, summary_statistic='mean') - assert np.array_equal(observed, isc(data, pairwise=True, - summary_statistic='mean')) + assert np.allclose(observed, isc(data, pairwise=True, + summary_statistic='mean'), + rtol=1e-03) logger.info("Finished testing circular time-shift") @@ -467,29 +566,31 @@ def test_phaseshift_isc(): observed, p, distribution = phaseshift_isc(data, pairwise=False) assert np.all(iscs[:, :2] > .5) assert np.all(iscs[:, -1] < .5) - assert p[0, 0] < .05 and p[0, 1] < .05 - assert p[0, 2] > .01 + assert p[0] < .05 and p[1] < .05 + assert p[2] > .01 iscs = isc(data, pairwise=True) observed, p, distribution = phaseshift_isc(data, pairwise=True) assert np.all(iscs[:, :2] > .5) assert np.all(iscs[:, -1] < .5) - assert p[0, 0] < .05 and p[0, 1] < .05 - assert p[0, 2] > .01 + assert p[0] < .05 and p[1] < .05 + assert p[2] > .01 # Check that ISC computation and permutation observed are same iscs = isc(data, pairwise=False) observed, p, distribution = phaseshift_isc(data, pairwise=False, summary_statistic='median') - assert np.array_equal(observed, isc(data, pairwise=False, - summary_statistic='median')) + assert np.allclose(observed, isc(data, pairwise=False, + summary_statistic='median'), + rtol=1e-03) # Check that ISC computation and permuation observed are same iscs = isc(data, pairwise=True) observed, p, distribution = phaseshift_isc(data, pairwise=True, summary_statistic='mean') - assert np.array_equal(observed, isc(data, pairwise=True, - summary_statistic='mean')) + assert np.allclose(observed, isc(data, pairwise=True, + summary_statistic='mean'), + rtol=1e-03) logger.info("Finished testing phase randomization") @@ -506,51 +607,350 @@ def test_isfc_options(): data = simulated_timeseries(n_subjects, n_TRs, n_voxels=n_voxels, data_type='array') - isfcs = isfc(data, pairwise=False, summary_statistic=None) + isfcs, iscs = isfc(data, pairwise=False, summary_statistic=None) + assert isfcs.shape == (n_subjects, n_voxels * (n_voxels - 1) / 2) + assert iscs.shape == (n_subjects, n_voxels) + + # Without vectorized upper triangle + isfcs = isfc(data, pairwise=False, summary_statistic=None, + vectorize_isfcs=False) + assert isfcs.shape == (n_subjects, n_voxels, n_voxels) # Just two subjects - isfcs = isfc(data[..., :2], pairwise=False, summary_statistic=None) + isfcs, iscs = isfc(data[..., :2], pairwise=False, summary_statistic=None) + assert isfcs.shape == (n_voxels * (n_voxels - 1) / 2,) + assert iscs.shape == (n_voxels,) + + isfcs = isfc(data[..., :2], pairwise=False, summary_statistic=None, + vectorize_isfcs=False) + assert isfcs.shape == (n_voxels, n_voxels) # ISFC with pairwise approach - isfcs = isfc(data, pairwise=True, summary_statistic=None) + isfcs, iscs = isfc(data, pairwise=True, summary_statistic=None) + assert isfcs.shape == (n_subjects * (n_subjects - 1) / 2, + n_voxels * (n_voxels - 1) / 2) + assert iscs.shape == (n_subjects * (n_subjects - 1) / 2, + n_voxels) + + isfcs = isfc(data, pairwise=True, summary_statistic=None, + vectorize_isfcs=False) + assert isfcs.shape == (n_subjects * (n_subjects - 1) / 2, + n_voxels, n_voxels) # ISFC with summary statistics - isfcs = isfc(data, pairwise=True, summary_statistic='mean') - isfcs = isfc(data, pairwise=True, summary_statistic='median') + isfcs, iscs = isfc(data, pairwise=True, summary_statistic='mean') + isfcs, iscs = isfc(data, pairwise=True, summary_statistic='median') # Check output p-values data = correlated_timeseries(20, 60, noise=.5, random_state=42) - isfcs = isfc(data, pairwise=False) - assert np.all(isfcs[0, 1, :] > .5) and np.all(isfcs[1, 0, :] > .5) - assert np.all(isfcs[:2, 2, :] < .5) and np.all(isfcs[2, :2, :] < .5) + isfcs = isfc(data, pairwise=False, vectorize_isfcs=False) + assert np.all(isfcs[:, 0, 1] > .5) and np.all(isfcs[:, 1, 0] > .5) + assert np.all(isfcs[:, :2, 2] < .5) and np.all(isfcs[:, 2, :2] < .5) - isfcs = isfc(data, pairwise=True) - assert np.all(isfcs[0, 1, :] > .5) and np.all(isfcs[1, 0, :] > .5) - assert np.all(isfcs[:2, 2, :] < .5) and np.all(isfcs[2, :2, :] < .5) + isfcs = isfc(data, pairwise=True, vectorize_isfcs=False) + assert np.all(isfcs[:, 0, 1] > .5) and np.all(isfcs[:, 1, 0] > .5) + assert np.all(isfcs[:, :2, 2] < .5) and np.all(isfcs[:, 2, :2] < .5) # Check that ISC and ISFC diagonal are identical iscs = isc(data, pairwise=False) - isfcs = isfc(data, pairwise=False) + isfcs = isfc(data, pairwise=False, vectorize_isfcs=False) for s in np.arange(len(iscs)): - assert np.allclose(isfcs[..., s].diagonal(), iscs[s, :], rtol=1e-03) + assert np.allclose(isfcs[s, ...].diagonal(), iscs[s, :], rtol=1e-03) + isfcs, iscs_v = isfc(data, pairwise=False) + assert np.allclose(iscs, iscs_v, rtol=1e-03) - # Check that ISC and ISFC diagonal are identical + # Check that ISC and ISFC diagonal are identical (pairwise) iscs = isc(data, pairwise=True) - isfcs = isfc(data, pairwise=True) + isfcs = isfc(data, pairwise=True, vectorize_isfcs=False) for s in np.arange(len(iscs)): - assert np.allclose(isfcs[..., s].diagonal(), iscs[s, :], rtol=1e-03) + assert np.allclose(isfcs[s, ...].diagonal(), iscs[s, :], rtol=1e-03) + isfcs, iscs_v = isfc(data, pairwise=True) + assert np.allclose(iscs, iscs_v, rtol=1e-03) + + # Generate 'targets' data and use for ISFC + data = simulated_timeseries(n_subjects, n_TRs, + n_voxels=n_voxels, data_type='array') + n_targets = 15 + targets_data = simulated_timeseries(n_subjects, n_TRs, + n_voxels=n_targets, + data_type='array') + isfcs = isfc(data, targets=targets_data, pairwise=False, + vectorize_isfcs=False) + assert isfcs.shape == (n_subjects, n_voxels, n_targets) + + # Ensure 'square' output enforced + isfcs = isfc(data, targets=targets_data, pairwise=False, + vectorize_isfcs=True) + assert isfcs.shape == (n_subjects, n_voxels, n_targets) + + # Check list input for targets + targets_data = simulated_timeseries(n_subjects, n_TRs, + n_voxels=n_targets, + data_type='list') + isfcs = isfc(data, targets=targets_data, pairwise=False, + vectorize_isfcs=False) + assert isfcs.shape == (n_subjects, n_voxels, n_targets) + + # Check that mismatching subjects / TRs breaks targets + targets_data = simulated_timeseries(n_subjects, n_TRs, + n_voxels=n_targets, + data_type='array') + + with pytest.raises(ValueError): + isfcs = isfc(data, targets=targets_data[..., :-1], + pairwise=False, vectorize_isfcs=False) + assert isfcs.shape == (n_subjects, n_voxels, n_targets) + + with pytest.raises(ValueError): + isfcs = isfc(data, targets=targets_data[:-1, ...], + pairwise=False, vectorize_isfcs=False) + + # Check targets for only 2 subjects + isfcs = isfc(data[..., :2], targets=targets_data[..., :2], + pairwise=False, summary_statistic=None) + assert isfcs.shape == (2, n_voxels, n_targets) + + isfcs = isfc(data[..., :2], targets=targets_data[..., :2], + pairwise=True, summary_statistic=None) + assert isfcs.shape == (2, n_voxels, n_targets) + + # Check that supplying targets enforces leave-one-out + isfcs_pw = isfc(data, targets=targets_data, pairwise=True, + vectorize_isfcs=False, tolerate_nans=False) + assert isfcs_pw.shape == (n_subjects, n_voxels, n_targets) logger.info("Finished testing ISFC options") +# Check for proper handling of NaNs in ISFC +def test_isfc_nans(): + + # Set parameters for toy time series data + n_subjects = 20 + n_TRs = 60 + n_voxels = 30 + random_state = 42 + + logger.info("Testing ISC options") + + data = simulated_timeseries(n_subjects, n_TRs, + n_voxels=n_voxels, data_type='array', + random_state=random_state) + + # Inject NaNs into data + data[0, 0, 0] = np.nan + + # Don't tolerate NaNs, should lose zeroeth voxel + isfcs_loo = isfc(data, pairwise=False, vectorize_isfcs=False, + tolerate_nans=False) + assert np.sum(np.isnan(isfcs_loo)) == n_subjects * (n_voxels * 2 - 1) + + # With vectorized ISFCs + isfcs_loo, iscs_loo = isfc(data, pairwise=False, vectorize_isfcs=True, + tolerate_nans=False) + assert np.sum(np.isnan(isfcs_loo)) == n_subjects * (n_voxels - 1) + + # Tolerate all NaNs, only subject with NaNs yields NaN + isfcs_loo = isfc(data, pairwise=False, vectorize_isfcs=False, + tolerate_nans=True) + assert np.sum(np.isnan(isfcs_loo)) == n_voxels * 2 - 1 + + isfcs_loo, iscs_loo = isfc(data, pairwise=False, vectorize_isfcs=True, + tolerate_nans=True) + assert np.sum(np.isnan(isfcs_loo)) == n_voxels - 1 + + # Pairwise approach shouldn't care + isfcs_pw_T = isfc(data, pairwise=True, vectorize_isfcs=False, + tolerate_nans=True) + isfcs_pw_F = isfc(data, pairwise=True, vectorize_isfcs=False, + tolerate_nans=False) + assert np.allclose(isfcs_pw_T, isfcs_pw_F, equal_nan=True) + assert (np.sum(np.isnan(isfcs_pw_T)) == + np.sum(np.isnan(isfcs_pw_F)) == + (n_voxels * 2 - 1) * (n_subjects - 1)) + + isfcs_pw_T, iscs_pw_T = isfc(data, pairwise=True, vectorize_isfcs=True, + tolerate_nans=True) + isfcs_pw_F, iscs_pw_T = isfc(data, pairwise=True, vectorize_isfcs=True, + tolerate_nans=False) + assert np.allclose(isfcs_pw_T, isfcs_pw_F, equal_nan=True) + assert (np.sum(np.isnan(isfcs_pw_T)) == + np.sum(np.isnan(isfcs_pw_F)) == + (n_voxels - 1) * (n_subjects - 1)) + + # Set proportion of nans to reject (70% and 90% non-NaN) + data[0, 0, :] = np.nan + data[0, 1, :n_subjects - int(n_subjects * .7)] = np.nan + data[0, 2, :n_subjects - int(n_subjects * .9)] = np.nan + + isfcs_loo_T = isfc(data, pairwise=False, vectorize_isfcs=False, + tolerate_nans=True) + isfcs_loo_F = isfc(data, pairwise=False, vectorize_isfcs=False, + tolerate_nans=False) + isfcs_loo_95 = isfc(data, pairwise=False, vectorize_isfcs=False, + tolerate_nans=.95) + isfcs_loo_90 = isfc(data, pairwise=False, vectorize_isfcs=False, + tolerate_nans=.90) + isfcs_loo_80 = isfc(data, pairwise=False, vectorize_isfcs=False, + tolerate_nans=.8) + isfcs_loo_70 = isfc(data, pairwise=False, vectorize_isfcs=False, + tolerate_nans=.7) + isfcs_loo_60 = isfc(data, pairwise=False, vectorize_isfcs=False, + tolerate_nans=.6) + assert (np.sum(np.isnan(isfcs_loo_F)) == + np.sum(np.isnan(isfcs_loo_95)) == 3420) + assert (np.sum(np.isnan(isfcs_loo_80)) == + np.sum(np.isnan(isfcs_loo_90)) == 2430) + assert (np.sum(np.isnan(isfcs_loo_T)) == + np.sum(np.isnan(isfcs_loo_60)) == + np.sum(np.isnan(isfcs_loo_70)) == 1632) + assert np.array_equal(np.sum(np.isnan(isfcs_loo_F), axis=0), + np.sum(np.isnan(isfcs_loo_95), axis=0)) + assert np.array_equal(np.sum(np.isnan(isfcs_loo_80), axis=0), + np.sum(np.isnan(isfcs_loo_90), axis=0)) + assert np.all((np.array_equal( + np.sum(np.isnan(isfcs_loo_T), axis=0), + np.sum(np.isnan(isfcs_loo_60), axis=0)), + np.array_equal( + np.sum(np.isnan(isfcs_loo_T), axis=0), + np.sum(np.isnan(isfcs_loo_70), axis=0)), + np.array_equal( + np.sum(np.isnan(isfcs_loo_60), axis=0), + np.sum(np.isnan(isfcs_loo_70), axis=0)))) + + isfcs_loo_T, _ = isfc(data, pairwise=False, vectorize_isfcs=True, + tolerate_nans=True) + isfcs_loo_F, _ = isfc(data, pairwise=False, vectorize_isfcs=True, + tolerate_nans=False) + isfcs_loo_95, _ = isfc(data, pairwise=False, vectorize_isfcs=True, + tolerate_nans=.95) + isfcs_loo_90, _ = isfc(data, pairwise=False, vectorize_isfcs=True, + tolerate_nans=.90) + isfcs_loo_80, _ = isfc(data, pairwise=False, vectorize_isfcs=True, + tolerate_nans=.8) + isfcs_loo_70, _ = isfc(data, pairwise=False, vectorize_isfcs=True, + tolerate_nans=.7) + isfcs_loo_60, _ = isfc(data, pairwise=False, vectorize_isfcs=True, + tolerate_nans=.6) + assert (np.sum(np.isnan(isfcs_loo_F)) == + np.sum(np.isnan(isfcs_loo_95)) == 1680) + assert (np.sum(np.isnan(isfcs_loo_80)) == + np.sum(np.isnan(isfcs_loo_90)) == 1194) + assert (np.sum(np.isnan(isfcs_loo_T)) == + np.sum(np.isnan(isfcs_loo_60)) == + np.sum(np.isnan(isfcs_loo_70)) == 802) + assert np.array_equal(np.sum(np.isnan(isfcs_loo_F), axis=0), + np.sum(np.isnan(isfcs_loo_95), axis=0)) + assert np.array_equal(np.sum(np.isnan(isfcs_loo_80), axis=0), + np.sum(np.isnan(isfcs_loo_90), axis=0)) + assert np.all((np.array_equal( + np.sum(np.isnan(isfcs_loo_T), axis=0), + np.sum(np.isnan(isfcs_loo_60), axis=0)), + np.array_equal( + np.sum(np.isnan(isfcs_loo_T), axis=0), + np.sum(np.isnan(isfcs_loo_70), axis=0)), + np.array_equal( + np.sum(np.isnan(isfcs_loo_60), axis=0), + np.sum(np.isnan(isfcs_loo_70), axis=0)))) + + data = simulated_timeseries(n_subjects, n_TRs, + n_voxels=n_voxels, data_type='array', + random_state=random_state) + + # Make sure voxel with NaNs across all subjects is always removed + data[0, 0, :] = np.nan + isfcs_loo_T = isfc(data, pairwise=False, vectorize_isfcs=False, + tolerate_nans=True) + isfcs_loo_F = isfc(data, pairwise=False, vectorize_isfcs=False, + tolerate_nans=False) + assert np.allclose(isfcs_loo_T, isfcs_loo_F, equal_nan=True) + assert (np.sum(np.isnan(isfcs_loo_T)) == + np.sum(np.isnan(isfcs_loo_F)) == + 1180) + + isfcs_pw_T = isfc(data, pairwise=True, vectorize_isfcs=False, + tolerate_nans=True) + isfcs_pw_F = isfc(data, pairwise=True, vectorize_isfcs=False, + tolerate_nans=False) + assert np.allclose(isfcs_pw_T, isfcs_pw_F, equal_nan=True) + + assert (np.sum(np.isnan(isfcs_pw_T)) == + np.sum(np.isnan(isfcs_pw_T)) == + 11210) + + # Check for NaN-handling in targets + n_targets = 15 + data = simulated_timeseries(n_subjects, n_TRs, + n_voxels=n_voxels, data_type='array', + random_state=random_state) + targets_data = simulated_timeseries(n_subjects, n_TRs, + n_voxels=n_targets, + data_type='array') + + # Inject NaNs into targets_data + targets_data[0, 0, 0] = np.nan + + # Don't tolerate NaNs, should lose zeroeth voxel + isfcs_loo = isfc(data, targets=targets_data, pairwise=False, + vectorize_isfcs=False, tolerate_nans=False) + assert np.sum(np.isnan(isfcs_loo)) == (n_subjects - 1) * (n_targets * 2) + + # Single NaN in targets will get averaged out with tolerate + isfcs_loo = isfc(data, targets=targets_data, pairwise=False, + vectorize_isfcs=False, tolerate_nans=True) + assert np.sum(np.isnan(isfcs_loo)) == 0 + + +def test_squareform_isfc(): + + # Set parameters for toy time series data + n_subjects = 20 + n_TRs = 60 + n_voxels = 30 + random_state = 42 + + logger.info("Testing ISC options") + + data = simulated_timeseries(n_subjects, n_TRs, + n_voxels=n_voxels, data_type='array', + random_state=random_state) + + # Generate square redundant ISFCs + isfcs_r = isfc(data, vectorize_isfcs=False) + assert isfcs_r.shape == (n_subjects, n_voxels, n_voxels) + + # Squareform these into condensed ISFCs and ISCs + isfcs_c, iscs_c = squareform_isfc(isfcs_r) + assert isfcs_c.shape == (n_subjects, n_voxels * (n_voxels - 1) / 2) + assert iscs_c.shape == (n_subjects, n_voxels) + + # Go back the other way and check it's the same + isfcs_new = squareform_isfc(isfcs_c, iscs_c) + assert np.array_equal(isfcs_r, isfcs_new) + + # Check against ISC function + assert np.allclose(isc(data), iscs_c, rtol=1e-03) + + # Check for two subjects + isfcs_r = isfc(data[..., :2], vectorize_isfcs=False) + assert isfcs_r.shape == (n_voxels, n_voxels) + isfcs_c, iscs_c = squareform_isfc(isfcs_r) + assert isfcs_c.shape == (n_voxels * (n_voxels - 1) / 2,) + assert iscs_c.shape == (n_voxels,) + assert np.array_equal(isfcs_r, squareform_isfc(isfcs_c, iscs_c)) + + if __name__ == '__main__': test_isc_input() test_isc_options() test_isc_output() + test_isc_nans() test_bootstrap_isc() test_permutation_isc() test_timeshift_isc() test_phaseshift_isc() test_isfc_options() + test_isfc_nans() + test_squareform_isfc() logger.info("Finished all ISC tests") diff --git a/tests/utils/test_utils.py b/tests/utils/test_utils.py index 1348790db..f37c51a7e 100644 --- a/tests/utils/test_utils.py +++ b/tests/utils/test_utils.py @@ -181,69 +181,216 @@ def test_center_mass_exp(): 'close to its mid-point' +def test_p_from_null(): + import numpy as np + from brainiak.utils.utils import p_from_null + + # Create random null and observed value in tail + null = np.random.randn(10000) + observed = np.ceil(np.percentile(null, 97.5) * 1000) / 1000 + + # Check that we catch improper side + with pytest.raises(ValueError): + _ = p_from_null(observed, null, side='wrong') + + # Check two-tailed p-value for observed + p_ts = p_from_null(observed, null) + assert np.isclose(p_ts, 0.05, atol=1e-02) + + # Check two-tailed p-value for observed + p_right = p_from_null(observed, null, side='right') + assert np.isclose(p_right, 0.025, atol=1e-02) + assert np.isclose(p_right, p_ts / 2, atol=1e-02) + + # Check two-tailed p-value for observed + p_left = p_from_null(observed, null, side='left') + assert np.isclose(p_left, 0.975, atol=1e-02) + assert np.isclose(1 - p_left, p_right, atol=1e-02) + assert np.isclose(1 - p_left, p_ts / 2, atol=1e-02) + + # Check 2-dimensional input (i.e., samples by voxels) + null = np.random.randn(10000, 3) + observed = np.ceil(np.percentile(null, 97.5, axis=0) * 1000) / 1000 + + # Check two-tailed p-value for observed + p_ts = p_from_null(observed, null, axis=0) + assert np.allclose(p_ts, 0.05, atol=1e-02) + + # Check two-tailed p-value for observed + p_right = p_from_null(observed, null, side='right', axis=0) + assert np.allclose(p_right, 0.025, atol=1e-02) + assert np.allclose(p_right, p_ts / 2, atol=1e-02) + + # Check two-tailed p-value for observed + p_left = p_from_null(observed, null, side='left', axis=0) + assert np.allclose(p_left, 0.975, atol=1e-02) + assert np.allclose(1 - p_left, p_right, atol=1e-02) + assert np.allclose(1 - p_left, p_ts / 2, atol=1e-02) + + # Check for exact test + p_ts = p_from_null(observed, null, exact=True, axis=0) + assert np.allclose(p_ts, 0.05, atol=1e-02) + + # Check two-tailed p-value for exact + p_right = p_from_null(observed, null, side='right', + exact=True, axis=0) + assert np.allclose(p_right, 0.025, atol=1e-02) + assert np.allclose(p_right, p_ts / 2, atol=1e-02) + + # Check two-tailed p-value for exact + p_left = p_from_null(observed, null, side='left', + exact=True, axis=0) + assert np.allclose(p_left, 0.975, atol=1e-02) + assert np.allclose(1 - p_left, p_right, atol=1e-02) + assert np.allclose(1 - p_left, p_ts / 2, atol=1e-02) + + 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 + from brainiak.utils.utils import phase_randomize - # 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)) + data = np.repeat(np.repeat(np.random.randn(60)[:, np.newaxis, np.newaxis], + 30, axis=1), + 20, axis=2) + assert np.array_equal(data[..., 0], data[..., 1]) - 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] + # Phase-randomize data across subjects (same across voxels) + shifted_data = phase_randomize(data, voxelwise=False, random_state=1) + assert shifted_data.shape == data.shape + assert not np.array_equal(shifted_data[..., 0], shifted_data[..., 1]) + assert not np.array_equal(shifted_data[..., 0], data[..., 0]) - assert np.isclose(abs(freq), abs(freq_pr)).all(), \ - "Amplitude spectrum not preserved under phase randomization" + # Check that uneven n_TRs doesn't explode + _ = phase_randomize(data[:-1, ...]) - assert abs(p_corr) < 0.03, \ - "Phases still correlated after randomization" + # Check that random_state returns same shifts + shifted_data_ = phase_randomize(data, voxelwise=False, random_state=1) + assert np.array_equal(shifted_data, shifted_data_) + shifted_data_ = phase_randomize(data, voxelwise=False, random_state=2) + assert not np.array_equal(shifted_data, shifted_data_) -def test_ecdf(): - from brainiak.utils.utils import ecdf - import numpy as np + # Phase-randomize subjects and voxels + shifted_data = phase_randomize(data, voxelwise=True, random_state=1) + assert shifted_data.shape == data.shape + assert not np.array_equal(shifted_data[..., 0], shifted_data[..., 1]) + assert not np.array_equal(shifted_data[..., 0], data[..., 0]) + assert not np.array_equal(shifted_data[:, 0, 0], shifted_data[:, 1, 0]) - x = np.array([1, 4, 3]) - cdf_fun = ecdf(x) + # Try with 2-dimensional input + shifted_data = phase_randomize(data[..., 0], + voxelwise=True, + random_state=1) + assert shifted_data.ndim == 2 + assert not np.array_equal(shifted_data[:, 0], shifted_data[:, 1]) - 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" + # Create correlated noisy data + corr_data = np.repeat(np.random.randn(60)[:, np.newaxis, np.newaxis], + 2, axis=2) + np.random.randn(60, 1, 2) + # Get correlation and frequency domain for data + corr_r = pearsonr(corr_data[:, 0, 0], + corr_data[:, 0, 1])[0] + corr_freq = fft(corr_data, axis=0) -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] + # Phase-randomize time series and get correlation/frequency + shifted_data = phase_randomize(corr_data) + shifted_r = pearsonr(shifted_data[:, 0, 0], + shifted_data[:, 0, 1])[0] + shifted_freq = fft(shifted_data, axis=0) - Y = X[:, 0] - Y_max = np.max(X[:, 1:], axis=0) - Y_min = np.min(X[:, 1:], axis=0) + # Check that phase-randomization reduces correlation + assert np.abs(shifted_r) < np.abs(corr_r) - p_1side = p_from_null(X, two_sided=False) - assert np.isclose(p_1side, [0.25, 1]).all(), "One-sided p value incorrect" + # Check that amplitude spectrum is preserved + assert np.allclose(np.abs(shifted_freq), np.abs(corr_freq)) - p_2side = p_from_null(X, two_sided=True) - assert np.isclose(p_2side, [0.5, 0]).all(), "Two-sided p value incorrect" - p_2side_m = p_from_null(Y, two_sided=True, - max_null_input=Y_max, - min_null_input=Y_min) - assert np.isclose(p_2side, p_2side_m).all(), "p_null differs with max/min" +def test_check_timeseries_input(): + import numpy as np + from itertools import combinations + from brainiak.utils.utils import _check_timeseries_input + + # Set a fixed vector for comparison + vector = np.random.randn(60) + + # List of subjects with one voxel/ROI + list_1d = [vector for _ in np.arange(10)] + (data_list_1d, n_TRs, + n_voxels, n_subjects) = _check_timeseries_input(list_1d) + assert n_TRs == 60 + assert n_voxels == 1 + assert n_subjects == 10 + + # Array of subjects with one voxel/ROI + array_2d = np.hstack([vector[:, np.newaxis] + for _ in np.arange(10)]) + (data_array_2d, n_TRs, + n_voxels, n_subjects) = _check_timeseries_input(array_2d) + assert n_TRs == 60 + assert n_voxels == 1 + assert n_subjects == 10 + + # List of 2-dimensional arrays + list_2d = [vector[:, np.newaxis] for _ in np.arange(10)] + (data_list_2d, n_TRs, + n_voxels, n_subjects) = _check_timeseries_input(list_2d) + assert n_TRs == 60 + assert n_voxels == 1 + assert n_subjects == 10 + + # Check if lists have mismatching size + list_bad = [list_2d[0][:-1, :]] + list_2d[1:] + with pytest.raises(ValueError): + (data_list_bad, _, _, _) = _check_timeseries_input(list_bad) + + # List of 3-dimensional arrays + list_3d = [vector[:, np.newaxis, np.newaxis] + for _ in np.arange(10)] + (data_list_3d, n_TRs, + n_voxels, n_subjects) = _check_timeseries_input(list_3d) + assert n_TRs == 60 + assert n_voxels == 1 + assert n_subjects == 10 + + # 3-dimensional array + array_3d = np.dstack([vector[:, np.newaxis] + for _ in np.arange(10)]) + (data_array_3d, n_TRs, + n_voxels, n_subjects) = _check_timeseries_input(array_3d) + assert n_TRs == 60 + assert n_voxels == 1 + assert n_subjects == 10 + + # Check that 4-dimensional input array throws error + array_4d = array_3d[..., np.newaxis] + with pytest.raises(ValueError): + (data_array_4d, _, _, _) = _check_timeseries_input(array_4d) + + # Check they're the same + for pair in combinations([data_list_1d, data_array_2d, + data_list_2d, data_list_3d, + data_array_3d], 2): + assert np.array_equal(pair[0], pair[1]) + + # List of multivoxel arrays + matrix = np.random.randn(60, 30) + list_mv = [matrix + for _ in np.arange(10)] + (data_list_mv, n_TRs, + n_voxels, n_subjects) = _check_timeseries_input(list_mv) + assert n_TRs == 60 + assert n_voxels == 30 + assert n_subjects == 10 + + # 3-dimensional array with multiple voxels + array_mv = np.dstack([matrix for _ in np.arange(10)]) + (data_array_mv, n_TRs, + n_voxels, n_subjects) = _check_timeseries_input(array_mv) + assert n_TRs == 60 + assert n_voxels == 30 + assert n_subjects == 10 + + assert np.array_equal(data_list_mv, data_array_mv)