# Copyright 2017 Intel Corporation

#

# Licensed under the Apache License, Version 2.0 (the "License");

# you may not use this file except in compliance with the License.

# You may obtain a copy of the License at

#

# http://www.apache.org/licenses/LICENSE-2.0

#

# Unless required by applicable law or agreed to in writing, software

# distributed under the License is distributed on an "AS IS" BASIS,

# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.

# See the License for the specific language governing permissions and

# limitations under the License.

"""Intersubject correlation (ISC) analysis

Functions for computing intersubject correlation (ISC) and related

analyses (e.g., intersubject funtional correlations; ISFC), as well

as statistical tests designed specifically for ISC analyses.

The implementation is based on the work in [Hasson2004]_, [Kauppi2014]_,

[Simony2016]_, [Chen2016]_, and [Nastase2019]_.

.. [Chen2016] "Untangling the relatedness among correlations, part I:

nonparametric approaches to inter-subject correlation analysis at the

group level.", G. Chen, Y. W. Shin, P. A. Taylor, D. R. Glen, R. C.

Reynolds, R. B. Israel, R. W. Cox, 2016, NeuroImage, 142, 248-259.

https://doi.org/10.1016/j.neuroimage.2016.05.023

.. [Hasson2004] "Intersubject synchronization of cortical activity

during natural vision.", U. Hasson, Y. Nir, I. Levy, G. Fuhrmann,

R. Malach, 2004, Science, 303, 1634-1640.

https://doi.org/10.1126/science.1089506

.. [Kauppi2014] "A versatile software package for inter-subject

correlation based analyses of fMRI.", J. P. Kauppi, J. Pajula,

J. Tohka, 2014, Frontiers in Neuroinformatics, 8, 2.

https://doi.org/10.3389/fninf.2014.00002

.. [Simony2016] "Dynamic reconfiguration of the default mode network

during narrative comprehension.", E. Simony, C. J. Honey, J. Chen, O.

Lositsky, Y. Yeshurun, A. Wiesel, U. Hasson, 2016, Nature Communications,

7, 12141. https://doi.org/10.1038/ncomms12141

.. [Nastase2019] "Measuring shared responses across subjects using

intersubject correlation." S. A. Nastase, V. Gazzola, U. Hasson,

C. Keysers, 2019, Social Cognitive and Affective Neuroscience, 14,

667-685. https://doi.org/10.1093/scan/nsz037

"""

# Authors: Sam Nastase, Christopher Baldassano, Qihong Lu,

# Mai Nguyen, and Mor Regev

# Princeton University, 2018

import math

import numpy as np

import logging

from scipy.spatial.distance import squareform

from itertools import combinations, permutations, product

from brainiak.fcma.util import compute_correlation

from brainiak.utils.utils import (

array_correlation, phase_randomize, p_from_null, _check_timeseries_input)

logger = logging.getLogger(__name__)

__all__ = [

"bootstrap_isc",

"compute_summary_statistic",

"isfc",

"isc",

"permutation_isc",

"phaseshift_isc",

"squareform_isfc",

"timeshift_isc",

]

MAX_RANDOM_SEED = 2**32 - 1

def isc(data, pairwise=False, summary_statistic=None, tolerate_nans=True):

"""Intersubject correlation

For each voxel or ROI, compute the Pearson correlation between each

subject's response time series and other subjects' response time series.

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 ISC values for N subjects

(leave-one-out) or N(N-1)/2 ISC values for each pair of N subjects,

corresponding to the upper triangle of the pairwise correlation matrix

(see scipy.spatial.distance.squareform). Alternatively, use either

'mean' or 'median' to compute summary statistic of ISCs (Fisher Z will

be applied if using mean). Input data should be a n_TRs by n_voxels by

n_subjects array (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. If only two

subjects are supplied, simply compute Pearson correlation (precludes

averaging in leave-one-out approach, and does not apply summary statistic).

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]_.

Parameters

----------

data : list or ndarray (n_TRs x n_voxels x n_subjects)

fMRI data for which to compute ISC

pairwise : bool, default: False

Whether to use pairwise (True) or leave-one-out (False) approach

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

ISC for each subject or pair (or summary statistic) per voxel

"""

# Check response time series input format

data, n_TRs, n_voxels, n_subjects = _check_timeseries_input(data)

# No summary statistic if only two subjects

if n_subjects == 2:

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)

# Compute correlation for only two participants

if n_subjects == 2:

# Compute correlation for each corresponding voxel

iscs_stack = array_correlation(data[..., 0],

data[..., 1])[np.newaxis, :]

# Compute pairwise ISCs using voxel loop and corrcoef for speed

elif pairwise:

# Swap axes for np.corrcoef

data = np.swapaxes(data, 2, 0)

# Loop through voxels

voxel_iscs = []

for v in np.arange(data.shape[1]):

voxel_data = data[:, v, :]

# Correlation matrix for all pairs of subjects (triangle)

iscs = squareform(np.corrcoef(voxel_data), checks=False)

voxel_iscs.append(iscs)

iscs_stack = np.column_stack(voxel_iscs)

# Compute leave-one-out ISCs

elif not pairwise:

# Loop through left-out subjects

iscs_stack = []

for s in np.arange(n_subjects):

# Correlation between left-out subject and mean of others

iscs_stack.append(array_correlation(

data[..., s],

mean(np.delete(data, s, axis=2), axis=2)))

iscs_stack = np.array(iscs_stack)

# 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:

iscs = compute_summary_statistic(iscs,

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, targets=None, pairwise=False, summary_statistic=None,

vectorize_isfcs=True, tolerate_nans=True):

"""Intersubject functional correlation (ISFC)

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)

or a list where each item is a n_TRs by n_voxels ndarray per subject.

Multiple input ndarrays must be the same shape. If a 2D array is supplied,

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). 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]_.

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), 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

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 : 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)

# 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),

return_nans=True)

isfcs = (isfcs + isfcs.T) / 2

isfcs = isfcs[..., np.newaxis]

summary_statistic = None

logger.info("Only two subjects! Computing ISFC between them.")

# Compute all pairwise ISFCs (only for symmetric approach)

elif pairwise:

isfcs = []

for pair in combinations(np.arange(n_subjects), 2):

isfc_pair = compute_correlation(np.ascontiguousarray(

data[..., pair[0]].T),

np.ascontiguousarray(

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)

# Compute ISFCs using leave-one-out approach

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(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 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=0)

# Throw away first dimension if singleton

if isfcs.shape[0] == 1:

isfcs = isfcs[0]

# 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):

"""Checks ISC inputs for statistical tests

Input ISCs should be n_subjects (leave-one-out approach) or

n_pairs (pairwise approach) by n_voxels or n_ROIs array or a 1D

array (or list) of ISC values for a single voxel or ROI. This

function is only intended to be used internally by other

functions in this module (e.g., bootstrap_isc, permutation_isc).

Parameters

----------

iscs : ndarray or list

ISC values

Returns

-------

iscs : ndarray

Array of ISC values

n_subjects : int

Number of subjects

n_voxels : int

Number of voxels (or ROIs)

"""

# Standardize structure of input data

if type(iscs) == list:

iscs = np.array(iscs)[:, np.newaxis]

elif isinstance(iscs, np.ndarray):

if iscs.ndim == 1:

iscs = iscs[:, np.newaxis]

# Check if incoming pairwise matrix is vectorized triangle

if pairwise:

try:

test_square = squareform(iscs[:, 0], force='tomatrix')

n_subjects = test_square.shape[0]

except ValueError:

raise ValueError("For pairwise input, ISCs must be the "

"vectorized triangle of a square matrix.")

elif not pairwise:

n_subjects = iscs.shape[0]

# Infer subjects, voxels and print for user to check

n_voxels = iscs.shape[1]

logger.info("Assuming {0} subjects with and {1} "

"voxel(s) or ROI(s) in bootstrap ISC test.".format(n_subjects,

n_voxels))

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

Computes either the 'mean' or 'median' across a set of ISCs. In the

case of the mean, ISC values are first Fisher Z transformed (arctanh),

averaged, then inverse Fisher Z transformed (tanh).

The implementation is based on the work in [SilverDunlap1987]_.

.. [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

Parameters

----------

iscs : list or ndarray

ISC values

summary_statistic : str, default: 'mean'

Summary statistic, 'mean' or 'median'

axis : None or int or tuple of ints, optional

Axis or axes along which the means are computed. The default is to

compute the mean of the flattened array.

Returns

-------

statistic : float or ndarray

Summary statistic of ISC values

"""

if summary_statistic not in ('mean', 'median'):

raise ValueError("Summary statistic must be 'mean' or 'median'")

# Compute summary statistic

if summary_statistic == 'mean':

statistic = np.tanh(np.nanmean(np.arctanh(iscs), axis=axis))

elif summary_statistic == 'median':

statistic = np.nanmedian(iscs, axis=axis)

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, side='right',

random_state=None):

"""One-sample group-level bootstrap hypothesis test for ISCs

For ISCs from one more voxels or ROIs, resample subjects with replacement

to construct a bootstrap distribution. Input is a list or ndarray of

ISCs for a single voxel/ROI, or an ISCs-by-voxels ndarray. ISC values

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 correlation matrix

(see scipy.stats.distance.squareform). Shifts bootstrap distribution by

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 (FPRs) for one-sample tests

in the pairwise approach. Note that the bootstrap hypothesis test may not

strictly control FPRs in the leave-one-out approach.

The implementation is based on the work in [Chen2016]_ and

[HallWilson1991]_.

.. [HallWilson1991] "Two guidelines for bootstrap hypothesis testing.",

P. Hall, S. R., Wilson, 1991, Biometrics, 757-762.

https://doi.org/10.2307/2532163

Parameters

----------

iscs : list or ndarray, ISCs by voxels array

ISC values for one or more voxels

pairwise : bool, default: False

Indicator of pairwise or leave-one-out, should match ISCs structure

summary_statistic : str, default: 'median'

Summary statistic, either 'median' (default) or 'mean'

n_bootstraps : int, default: 1000

Number of bootstrap samples (subject-level with replacement)

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

Returns

-------

observed : float, median (or mean) ISC value

Summary statistic for actual ISCs

ci : tuple, bootstrap confidence intervals

Confidence intervals generated from bootstrap distribution

p : float, p-value

p-value based on bootstrap hypothesis test

distribution : ndarray, n_bootstraps by voxels

Bootstrap distribution

"""

# Standardize structure of input data

iscs, n_subjects, n_voxels = _check_isc_input(iscs, pairwise=pairwise)

# Check for valid summary statistic

if summary_statistic not in ('mean', 'median'):

raise ValueError("Summary statistic must be 'mean' or 'median'")

# Compute summary statistic for observed ISCs

observed = compute_summary_statistic(iscs,

summary_statistic=summary_statistic,

axis=0)

# Set up an empty list to build our bootstrap distribution

distribution = []

# Loop through n bootstrap iterations and populate distribution

for i in np.arange(n_bootstraps):

# 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)

# Randomly sample subject IDs with replacement

subject_sample = sorted(prng.choice(np.arange(n_subjects),

size=n_subjects))

# Squareform and shuffle rows/columns of pairwise ISC matrix to

# to retain correlation structure among ISCs, then get triangle

if pairwise:

# Loop through voxels

isc_sample = []

for voxel_iscs in iscs.T:

# Square the triangle and fill diagonal

voxel_iscs = squareform(voxel_iscs, force='tomatrix')

np.fill_diagonal(voxel_iscs, 1)

# Shuffle square correlation matrix and get triangle

voxel_sample = voxel_iscs[subject_sample, :][:, subject_sample]

voxel_sample = squareform(voxel_sample, checks=False)

# Censor off-diagonal 1s for same-subject pairs

voxel_sample[voxel_sample == 1.] = np.nan

isc_sample.append(voxel_sample)

isc_sample = np.column_stack(isc_sample)

# Get simple bootstrap sample if not pairwise

elif not pairwise:

isc_sample = iscs[subject_sample, :]

# Compute summary statistic for bootstrap ISCs per voxel

# (alternatively could construct distribution for all voxels

# then compute statistics, but larger memory footprint)

distribution.append(compute_summary_statistic(

isc_sample,

summary_statistic=summary_statistic,

axis=0))

# Update random state for next iteration

random_state = np.random.RandomState(prng.randint(0, MAX_RANDOM_SEED,

dtype=np.int64))

# Convert distribution to numpy array

distribution = np.array(distribution)

# Compute CIs of median from bootstrap distribution (default: 95%)

ci = (np.percentile(distribution, (100 - ci_percentile)/2, axis=0),

np.percentile(distribution, ci_percentile + (100 - ci_percentile)/2,

axis=0))

# Shift bootstrap distribution to 0 for hypothesis test

shifted = distribution - observed

# Get p-value for actual median from shifted distribution

p = p_from_null(observed, shifted,

side=side, exact=False,

axis=0)

return observed, ci, p, distribution

def _check_group_assignment(group_assignment, n_subjects):

if type(group_assignment) == list:

pass

elif type(group_assignment) == np.ndarray:

group_assignment = group_assignment.tolist()

else:

logger.info("No group assignment provided, "

"performing one-sample test.")

if group_assignment and len(group_assignment) != n_subjects:

raise ValueError("Group assignments ({0}) "

"do not match number of subjects ({1})!".format(

len(group_assignment), n_subjects))

return group_assignment

def _get_group_parameters(group_assignment, n_subjects, pairwise=False):

# Set up dictionary to contain group info

group_parameters = {'group_assignment': group_assignment,

'n_subjects': n_subjects,

'group_labels': None, 'groups': None,

'sorter': None, 'unsorter': None,

'group_matrix': None, 'group_selector': None}

# Set up group selectors for two-group scenario

if group_assignment and len(np.unique(group_assignment)) == 2:

group_parameters['n_groups'] = 2

# Get group labels and counts

group_labels = np.unique(group_assignment)

groups = {group_labels[0]: group_assignment.count(group_labels[0]),

group_labels[1]: group_assignment.count(group_labels[1])}

# For two-sample pairwise approach set up selector from matrix

if pairwise:

# Sort the group_assignment variable if it came in shuffled

# so it's easier to build group assignment matrix

sorter = np.array(group_assignment).argsort()

unsorter = np.array(group_assignment).argsort().argsort()

# Populate a matrix with group assignments

upper_left = np.full((groups[group_labels[0]],

groups[group_labels[0]]),

group_labels[0])

upper_right = np.full((groups[group_labels[0]],

groups[group_labels[1]]),

np.nan)

lower_left = np.full((groups[group_labels[1]],

groups[group_labels[0]]),

np.nan)

lower_right = np.full((groups[group_labels[1]],

groups[group_labels[1]]),

group_labels[1])

group_matrix = np.vstack((np.hstack((upper_left, upper_right)),

np.hstack((lower_left, lower_right))))

np.fill_diagonal(group_matrix, np.nan)

group_parameters['group_matrix'] = group_matrix

# Unsort matrix and squareform to create selector

group_parameters['group_selector'] = squareform(

group_matrix[unsorter, :][:, unsorter],

checks=False)

group_parameters['sorter'] = sorter

group_parameters['unsorter'] = unsorter

# If leave-one-out approach, just user group assignment as selector

else:

group_parameters['group_selector'] = group_assignment

# Save these parameters for later

group_parameters['groups'] = groups

group_parameters['group_labels'] = group_labels

# Manage one-sample and incorrect group assignments

elif not group_assignment or len(np.unique(group_assignment)) == 1:

group_parameters['n_groups'] = 1

# If pairwise initialize matrix of ones for sign-flipping

if pairwise:

group_parameters['group_matrix'] = np.ones((

group_parameters['n_subjects'],

group_parameters['n_subjects']))

elif len(np.unique(group_assignment)) > 2:

raise ValueError("This test is not valid for more than "

"2 groups! (got {0})".format(

len(np.unique(group_assignment))))

else:

raise ValueError("Invalid group assignments!")

return group_parameters

def _permute_one_sample_iscs(iscs, group_parameters, i, pairwise=False,

summary_statistic='median', group_matrix=None,

exact_permutations=None, prng=None):

"""Applies one-sample permutations to ISC data

Input ISCs should be n_subjects (leave-one-out approach) or

n_pairs (pairwise approach) by n_voxels or n_ROIs array.

This function is only intended to be used internally by the

permutation_isc function in this module.

Parameters

----------

iscs : ndarray or list

ISC values

group_parameters : dict

Dictionary of group parameters

i : int

Permutation iteration

pairwise : bool, default: False

Indicator of pairwise or leave-one-out, should match ISCs variable

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

Initial random seed

Returns

-------

isc_sample : ndarray

Array of permuted ISC values

"""

# Randomized sign-flips

if exact_permutations:

sign_flipper = np.array(exact_permutations[i])

else:

sign_flipper = prng.choice([-1, 1],

size=group_parameters['n_subjects'],

replace=True)

# If pairwise, apply sign-flips by rows and columns

if pairwise:

matrix_flipped = (group_parameters['group_matrix'] * sign_flipper

* sign_flipper[

:, np.newaxis])

sign_flipper = squareform(matrix_flipped, checks=False)

# Apply flips along ISC axis (same across voxels)

isc_flipped = iscs * sign_flipper[:, np.newaxis]

# Get summary statistics on sign-flipped ISCs

isc_sample = compute_summary_statistic(

isc_flipped,

summary_statistic=summary_statistic,

axis=0)

return isc_sample

def _permute_two_sample_iscs(iscs, group_parameters, i, pairwise=False,

summary_statistic='median',

exact_permutations=None, prng=None):

"""Applies two-sample permutations to ISC data

Input ISCs should be n_subjects (leave-one-out approach) or

n_pairs (pairwise approach) by n_voxels or n_ROIs array.

This function is only intended to be used internally by the

permutation_isc function in this module.

Parameters

----------

iscs : ndarray or list

ISC values

group_parameters : dict

Dictionary of group parameters

i : int

Permutation iteration

pairwise : bool, default: False

Indicator of pairwise or leave-one-out, should match ISCs variable

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

Initial random seed