Synchrony & Alignment#

Written by Luke Chang and based on a tutorial on ISC

Synchrony between individuals happens at several levels from behavior to brain activity (Nummenmaa et al., 2018, Nastase et al., 2019). To an observer, synchrony during interaction or joint motion (Hale et al., 2020) can reflect prosocial qualities such rapport (Miles et al., 2009) or affiliation (Hove and Risen, 2009). During physically arousing ritualistic experiences, observers may selectively synchronize their heart rates while observing related, but not unrelated individuals participating in the ritual (Konvalinka et al., 2011) and the degree of affective synchrony can predict the perceived social connection between two individuals (Cheong et al., 2020). Synchrony of brain activity is associated with, among other things, shared psychological perspectives toward a stimulus (Lahnakoski et al., 2014, Yeshurun et al., 2017, Yeshurun et al., 2017, Kang & Wheatley, 2017, Cheong et al., 2020) and friendship (Parkinson et al., 2018), and may also be disturbed in psychiatric conditions ranging from developmental conditions such as autism (Hasson et al., 2009; Salmi et al., 2013, Byrge et al., 2015) to more acute conditions such as first-episode psychosis (Mäntylä et al., 2018). Thus, measures of synchrony can offer a simple window to many psychological processes.

In brain imaging, synchrony of brain activity is most commonly measured using intersubject correlations (ISC; Hasson et al., 2004). As the name implies, this method calculates linear correlations between participants and derives summary statistics from these correlations to measure the level of similarity of brain activity. Overall, the brain activity measured with fMRI during naturalistic stimulation conditions can be thought to consist of four main sources: (1) stimulus-driven brain activity between individuals shared by all or most of the participants, (2) individual or idiosyncratic activity elicited by the stimulus, (3) intrinsic activity that is not time-locked to the stimulus, and (4) noise from various sources. The idea behind ISC is to identify brain activity that is shared by a large proportion of the participants (category 1). Thus, this method evaluates how much of an individual’s brain activity is explained by this shared component. By contrast, if smaller groups of participants, e.g. pairs of friends within the study, share similar individual activity patterns (category 2), it may be better captured by the dyadic values in the pairwise matrices using techniques such as IS-RSA (Parkinson et al., 2018, Chen et al., 2020, Finn et al., 2020). Generally, the third category of activity is not readily detected by ISC, but can be revealed using alternative techniques that don’t necessarily assume time locking such as hidden markov models (Chang et al., 2021). Using inventive experimental designs, e.g. during verbal recall of previously experienced stimuli (Chen et al., 2017), it is still possible to extract shared brain activity patterns by temporally reorganizing the data, even when the original experiences of the participants were out of sync. The optimal choice of analysis will depend on the research question and the type of shared activity patterns that are of particular interest. Most commonly, ISCs are calculated locally within each voxel (or region), but this approach has also been extended to functional connectivity (Simony et al., 2016). Intersubject correlations give a summary statistic of synchrony over long periods of time, usually over an entire imaging session. However, the level of synchrony may change considerably from one moment to the next with the experimental condition and can be investigated using time-varying measures of synchrony.

In this tutorial, we will cover:

  • how to compute pairwise ISC

  • how to perform hypothesis tests with ISC data.org/features/notebooks/timecorr.html) tutorial)

  • how to perform dynamic ISC

  • how to apply isc to any feature space

  • how to apply isc to a low dimensional feature space

Computing Intersubject Correlation#

Pairwise vs Average Response#

Generally, ISCs are calculated using one of two main approaches. First, one calculates pairwise correlations between all participant pairs to build a full intersubject correlation matrix. The second approach, uses the average activity timecourse of the other participants as a model for each individual left out participant. This procedure produces individual, rather than pairwise, spatial maps of similarity reflecting how similarly, or “typically”, each person’s brain activates. These maps lend themselves to similar analyses as one might perform with first level results of a traditional general linear model (GLM) analysis. However, some of the individual variability is lost with the latter approach, and thus the ISC values are typically of much higher magnitude compared to the pairwise matrices.

Summary Statistic#

The next step after computing either the pairwise or average similarity across participants is to summarize the overall level of synchrony across participants. A straightforward approach is to compute the mean correlation. However, to make the correlation coefficients more normally distributed across the range of values, the Fisher’s Z transformation (inverse hyperbolic tangent) is typically applied before computing the mean correlation. The Fisher’s Z transformation mainly affects the higher end of absolute correlation values, eventually stretching the correlation coefficient 1 to infinity. However, with the typical range of values of pairwise ISCs, the effects of this transformation are relatively small reaching ~10% at the higher end of the scale of r=0.5. More recently, it has been suggested that computing the median, particularly when using the pairwise approach, provides a more accurate summary of the correlation values (Chen et al., 2016).

Hypothesis tests#

Performing hypothesis tests that appropriately account for the false positive rate can be tricky with ISC because of the dependence between the pairwise correlation values and the inflated number of variables in the pairwise correlation matrices. Though there have been proposals to use mixed-effects models for a parametric solution (Chen et al., 2017), we generally recommend using non-parametric statistics when evaluating the significance of these correlations.

There are two general non-parametric approaches to performing hypothesis tests with ISC. The first method is a permutation or randomization method, achieved by creating surrogate data and repeating the same analysis many times to build an empirical null distribution (e.g. 5-10k iterations). However, to meet the exchangeability assumption, it is important to consider the temporal dependence structure. Surrogate data can be created by circularly shifting the timecourses of the participants, or scrambling the phases of the Fourier transform of the signals and transforming these scrambled signals back to the time domain (Theiler et al., 1992, Lancaster et al., 2018). Various blockwise scrambling techniques have also been applied and autoregressive models have been proposed to create artificial data for statistical inference . These approaches have the benefit that, when properly designed, they retain important characteristics of the original signal, such as the frequency content and autocorrelation, while removing temporal synchrony in the data.

To illustrate these permutation based approaches, the animation below depicts the process of creating these null distributions and compares these to a similar distribution built based on real resting state data of the same duration in the same participants recorded just prior to the movie data in the same imaging study. Resting state is an ideal condition for demonstrating the true null distribution of no synchrony intersubject correlation as it involves no external synchronizing factors apart from the repeating noise of the scanner gradients, which are generally not of interest to us. Thus, any correlations in the resting state data arise by chance. As can be seen, the null distributions based on the surrogate data follow the distribution of resting state ISCs well as the number of iterations increases. However, the distributions are sometimes considered too liberal.

null_distribution

The second approach employs a subject-wise bootstrap on the pairwise similarity matrices. Essentially, participants are randomly sampled with replacement and then a new similarity matrix is computed with these resampled participants. As a consequence of the resampling procedure, sometimes the same subjects are sampled multiple times, which introduces correlation values of 1 off the diagonal. Summarizing the ISC using the median can minimize the impact of these outliers. These values are then shifted by the real summary statistics to produce an approximately zero-centered distribution. Python implementations in Brainiak and nltools convert these values to NaNs by default so that they are not included in the overall ISC summary statistic. If you would like to learn more about resampling methods, we encourage you to read our brief introduction available on the dartbrains course.

Getting Started#

Before getting started with this tutorial, we need to make sure you have the necessary software installed and data downloaded.

Software#

This tutorial requires the following Python packages to be installed. See the Software Installation tutorial for more information.

  • seaborn

  • matplotlib

  • numpy

  • scipy

  • pandas

  • nltools

  • nilearn

  • sklearn

  • datalad

Let’s now load the modules we will be using for this tutorial.

%matplotlib inline

import os
from tqdm.notebook import tqdm
import glob
import numpy as np
from numpy.fft import fft, ifft, fftfreq
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib import gridspec
from matplotlib.animation import FuncAnimation
import seaborn as sns
from nltools.data import Brain_Data, Adjacency
from nltools.mask import expand_mask, roi_to_brain, create_sphere
from nltools.stats import (isc, isc_group, isfc, isps, fdr, threshold, phase_randomize, circle_shift, _butter_bandpass_filter, _phase_mean_angle, _phase_vector_length, zscore, align)
from nilearn.plotting import view_img_on_surf, view_img, plot_stat_map
from sklearn.metrics import pairwise_distances
from sklearn.utils import check_random_state
from scipy.stats import ttest_1samp, pearsonr
from scipy.signal import hilbert
from IPython.display import HTML
from fnl_tools.utils import rec_to_time
import datalad.api as dl
import warnings

warnings.simplefilter('ignore')

Data#

This tutorial will be using the Friday Night Lights dataset and will require downloading the Average ROI csv files.

We have already extracted the data for you to make this easier and have written out the average activity within each ROI into a separate csv file for each participant. If you would like to get practice doing this yourself, here is the code we used. Note that we are working with the hdf5 files as they load much faster than the nifti images, but either one will work for this example.

mask = Brain_Data('http://neurovault.org/media/images/2099/Neurosynth%20Parcellation_0.nii.gz')

file_list = glob.glob(os.path.join(download_path, 'ds003521', 'derivatives', 'denoised', 'smoothed', '*hdf5'))
for f in tqdm(file_list):
    sub_dat = Brain_Data(f)
    sub_roi = sub_dat.extract_roi(mask)
    sub_roi = pd.DataFrame(sub_roi).T
    output_name = f"{os.path.basename(f).split('_bold.hdf5')[0]}_average_roi.csv"
    sub_roi.to_csv(os.path.join(download_path, 'ds003521', 'derivatives', 'denoised', 'smoothed', output_name), index=False)

For now, we will be downloading the csv files directly from github, but you should also be able to do this using datalad, once some of the OpenNeuro syncing has completed.

You will want to change data_dir to wherever you have installed the Friday Night Lights datalad repository. We will initialize a datalad dataset instance and get the files we need for this tutorial. If you’ve already downloaded everything, you can skip this cell. See the Download Data Tutorial for more information about how to install and use datalad.

Here is example code:


data_dir = '/Users/lukechang/Downloads/FNL'

# If dataset hasn't been installed, clone from GIN repository
if not os.path.exists(data_dir):
    dl.clone(source='https://github.com/OpenNeuroDatasets/ds003521.git', path=data_dir)

# Initialize dataset
ds = dl.Dataset(data_dir)

# Get Cropped & Denoised CSV Files
result = ds.get(glob.glob(os.path.join(data_dir, 'derivatives', 'denoised', 'smoothed', f'*average_roi.csv')))

For this tutorial, you can also access the data directly from github.

sub_list = ['sub-sid000705',
 'sub-sid000708',
 'sub-sid000707',
 'sub-sid000678',
 'sub-sid000677',
 'sub-sid000475',
 'sub-sid000375',
 'sub-sid000681',
 'sub-sid000683',
 'sub-sid000498',
 'sub-sid000697',
 'sub-sid000698',
 'sub-sid000695',
 'sub-sid000692',
 'sub-sid000216',
 'sub-sid000570',
 'sub-sid000669',
 'sub-sid000688',
 'sub-sid000668',
 'sub-sid000573',
 'sub-sid000667',
 'sub-sid000671',
 'sub-sid000217',
 'sub-sid000529',
 'sub-sid000679',
 'sub-sid000476',
 'sub-sid000704',
 'sub-sid000706',
 'sub-sid000710',
 'sub-sid000680',
 'sub-sid000682',
 'sub-sid000494',
 'sub-sid000694',
 'sub-sid000699']

    
sub_timeseries = {}
for sub in tqdm(sub_list):
    sub_data = pd.read_csv(f'https://raw.githubusercontent.com/mind-tutorials/mind_book/master/data/fnl_average_roi/{sub}_task-movie_run-1_space-MNI152NLin2009cAsym_desc-preproc_trim_smooth6_denoised_average_roi.csv')
    sub_timeseries[sub] = sub_data
data_dir = '/Users/lukechang/Downloads/FNL'

# If dataset hasn't been installed, clone from GIN repository
if not os.path.exists(data_dir):
    dl.clone(source='https://github.com/OpenNeuroDatasets/ds003521.git', path=data_dir)

# Initialize dataset
ds = dl.Dataset(data_dir)

# Get Cropped & Denoised CSV Files
result = ds.get(glob.glob(os.path.join(data_dir, 'derivatives', 'denoised', 'smoothed', f'*average_roi.csv')))

# Load CSV files
file_list = glob.glob(os.path.join(data_dir, 'derivatives', 'denoised', 'smoothed', f'*average_roi.csv'))

sub_timeseries = {}
for f in tqdm(file_list):
    sub = os.path.basename(f).split('_')[0]
    sub_timeseries[sub] = pd.read_csv(f)
action summary:
  get (notneeded: 34)

Pairwise Synchrony#

What if we were interested in examining if two subjects synchronized within their dmPFC while watching the movie?

There are many ways to test this, but a very straightforward approach is to simply correlate the timeseries of average activity within the dmPFC to the average activity from a different participant.

Ok, now let’s download the k=50 whole brain meta-analytic parcellation of the neurosynth database (de la Vega, 2016) from neurovault. Each ROI is indicated with a unique integer. We can expand this mask into 50 separate binary masks with the expand_mask function. The dmPFC corresponds to an integer index of 2.

mask = Brain_Data('http://neurovault.org/media/images/2099/Neurosynth%20Parcellation_0.nii.gz')
mask_x = expand_mask(mask)

mask.plot()
../_images/f12c73fac293abf04b1fb964d375b8b83945744b8ce4392faa52a41f07c2a8df.png
sub_1 = 'sub-sid000705'
sub_2 = 'sub-sid000708'

roi = 2

mask_x[roi].plot()

r, p = pearsonr(sub_timeseries[sub_1].iloc[:, roi], sub_timeseries[sub_2].iloc[:, roi])

f,a = plt.subplots(figsize=(15,3))
a.plot(sub_timeseries[sub_1].iloc[:, roi])
a.plot(sub_timeseries[sub_2].iloc[:, roi])
a.legend([sub_1, sub_2])
a.set_ylabel('Average activity')
a.set_title(f'Average time series within dmPFC - r={r:.02}')
Text(0.5, 1.0, 'Average time series within dmPFC - r=0.11')
../_images/878d81102906072aed1e20308e262b5236d0250ffd66e7feffc90077d1c414a6.png ../_images/01eae95f08d264a60582ae7f0c2b3a75da9435078fb266533648a2a6588b7b32.png

As you can see this is very simple to assess overall synchrony. One thing that can be slightly more tricky is how to perform hypothesis tests to assess statistical significance. The p-values from the pearsonr scipy function assumes that each row of the data are independent, which is likely violated in this example as fMRI data tends to be autocorrelated. There are numerous ways to account for this using resampling statistics. We recommend using circle shifting or phase randomization, which are described in more detail below.

Intersubject correlation (ISC) Analysis#

Ok, now what if we were interested in assessing the overall level of synchrony across our entire sample? We can extend our pairwise synchrony concept to all participants and then summarizing the total pairwise synchrony using the mean or median.

let’s see how ISC works with real data. In this tutorial, we will be computing ISC on average activity within 50 ROIs using data from the Friday Night Lights dataset using the nltools Python toolbox. We also recommend checking out the Brainiak tutorial, and the Matlab ISC-Toolbox (Kauppi et al., 2014).

Now, let’s pick a single ROI to demonstrate how to perform ISC analyses. We can pick any region, but let’s continue with the dmPFC (roi=2), feel free to play with different regions as you work with the code.

We can plot the mask and create a new pandas DataFrame that has the average vmPFC activity for each participant.

roi = 2

mask_x[roi].plot()

def get_subject_roi(data, roi):
    sub_rois = {}
    for sub in data:
        sub_rois[sub] = data[sub].iloc[:, roi]
    return pd.DataFrame(sub_rois)

sub_rois = get_subject_roi(sub_timeseries, roi)
sub_rois.head()
sub-sid000705 sub-sid000708 sub-sid000707 sub-sid000678 sub-sid000677 sub-sid000475 sub-sid000375 sub-sid000681 sub-sid000683 sub-sid000498 ... sub-sid000679 sub-sid000476 sub-sid000704 sub-sid000706 sub-sid000710 sub-sid000680 sub-sid000682 sub-sid000494 sub-sid000694 sub-sid000699
0 -1.257059 -1.858981 0.055510 -2.927054 0.001443 0.976214 -0.091158 1.122817 0.531115 -1.476508 ... -0.816042 -1.364593 0.133366 -2.057195 1.901717 -0.849037 -0.132114 -4.361229e+00 -0.863538 -1.691733
1 0.025426 -0.666974 -0.661008 -1.038320 -1.200573 0.360406 -0.165949 -0.136538 0.619254 -1.421292 ... -0.026740 -2.155036 0.347466 -2.337125 0.614552 -0.656935 -0.109674 -1.787996e-12 -0.656501 -1.883794
2 1.147395 -1.490594 -1.106923 -1.106748 -0.171776 0.517466 -0.011003 0.039086 -0.711524 -1.461454 ... 1.600135 -2.700619 0.747714 -1.796568 1.175404 -0.381551 0.379087 -3.886152e-12 0.714567 -1.965230
3 -0.312306 0.235479 -0.832777 -2.265747 2.064173 1.440031 0.581691 1.074168 -0.916229 -0.676008 ... 2.817513 -2.574495 1.161320 -1.159793 1.798948 -0.109981 0.521651 3.098448e-12 0.582420 -1.078426
4 -0.960072 0.649027 -0.094664 -2.235001 1.333921 1.734485 0.790967 0.677707 -0.524876 0.147219 ... 2.265140 -0.002140 0.662223 -0.648039 1.305884 0.865690 -0.304318 -5.992719e-13 0.707843 -1.072950

5 rows × 34 columns

../_images/878d81102906072aed1e20308e262b5236d0250ffd66e7feffc90077d1c414a6.png

Circle Shift Randomization#

To perform ISC we will be using the nltools.stats.isc function. There are three different methods implemented to perform hypothesis tests (i.e., circular shifting data, phase randomization, and subject-wise bootstrap). We will walk through an example of how to run each one.

The idea behind circular shifting the data is to generate random surrogate data that has the same autoregressive and temporal properties of the original data (Lancaster et al., 2018). This is fairly straightforward and involves randomly selecting a time point to become the new beginning of the timeseries and then concatenating the rest of the data at the end so that it has the same length as the original data. Of course, there will potentially be a sudden change in the data where the two parts were merged.

To demonstrate how this works and that the circular shifted data has the same spectral properties of the original data, we will plot one subject’s time series and shift it using the nltools.stats.circle_shift function. Next to both timeseries we plot the coefficients from a fast fourier transform.

sub_1 = 'sub-sid000705'
sampling_freq = .5

f,a = plt.subplots(nrows=2, ncols=2, figsize=(15, 5))
a[0,0].plot(sub_rois[sub], linewidth=2)
a[0,0].set_ylabel('Avg Activity', fontsize=16)
a[0,1].set_xlabel('Time (TR)', fontsize=18)
a[0,0].set_title('Observed Data', fontsize=16)

fft_data = fft(sub_rois[sub_1])
freq = fftfreq(len(fft_data), 1/sampling_freq)
n_freq = int(np.floor(len(fft_data)/2))
a[0,1].plot(freq[:n_freq], np.abs(fft_data)[:n_freq], linewidth=2)
a[0,1].set_xlabel('Frequency (Hz)', fontsize=18)
a[0,1].set_ylabel('Amplitude', fontsize=18)
a[0,1].set_title('Power Spectrum', fontsize=18)

circle_shift_data = circle_shift(sub_rois[sub_1])
a[1,0].plot(circle_shift_data, linewidth=2, color='red')
a[1,0].set_ylabel('Avg Activity', fontsize=16)
a[1,0].set_xlabel('Time (TR)', fontsize=16)
a[1,0].set_title('Circle Shifted Data', fontsize=16)

fft_circle = fft(circle_shift_data)
a[1,1].plot(freq[:n_freq], np.abs(fft_circle)[:n_freq], linewidth=2, color='red')
a[1,1].set_xlabel('Frequency (Hz)', fontsize=18)
a[1,1].set_ylabel('Amplitude', fontsize=18)
a[1,1].set_title('Circle Shifted Power Spectrum', fontsize=18)

plt.tight_layout()
../_images/93d97a84d62411310f98a0910158e2924c58a1f7d4ed09f0ae977c299c1dd6ac.png

Now that you understand how the circular shifting method works with a single time series, let’s compute the ISC of the Friday Night Lights viewing data within the dmPFC roi with 5000 permutations. The function outputs a dictionary that contains the ISC values, the p-value, the 95% confidence intervals, and optionally returns the 5000 samples. All of the permutation and bootstraps are parallelized and will use as many cores as your computer has available.

stats_circle = isc(sub_rois, method='circle_shift', n_samples=5000, return_null=True)

print(f"ISC: {stats_circle['isc']:.02}, p = {stats_circle['p']:.03}")
ISC: 0.085, p = 0.0002
stats_circle.keys()
dict_keys(['isc', 'p', 'ci', 'null_distribution'])

Phase Randomization#

Next we will show how phase randomization works. This method projects the data into frequency space using a fast fourier transform and randomizes the phase and then projects the data back into the time domain (Theiler et al., 1992, Lancaster et al., 2018). If you would like to learn more about Fourier transforms, check out our tutorial on Dartbrains. Similar to the circular shifting method, this generates a random surrogate of the data, while maintaining a similar temporal and autoregressive structure as the original data. We will generate the same plots from the example above to illustrate how this works using the nltools.stats.phase_randomize function.

sub_1 = 'sub-sid000705'
sampling_freq = .5

f,a = plt.subplots(nrows=2, ncols=2, figsize=(15, 5))
a[0,0].plot(sub_rois[sub_1], linewidth=2)
a[0,0].set_ylabel('Avg Activity', fontsize=16)
a[0,1].set_xlabel('Time (TR)', fontsize=18)
a[0,0].set_title('Observed Data', fontsize=16)

fft_data = fft(sub_rois[sub_1])
freq = fftfreq(len(fft_data), 1/sampling_freq)
n_freq = int(np.floor(len(fft_data)/2))
a[0,1].plot(freq[:n_freq], np.abs(fft_data)[:n_freq], linewidth=2)
a[0,1].set_xlabel('Frequency (Hz)', fontsize=18)
a[0,1].set_ylabel('Amplitude', fontsize=18)
a[0,1].set_title('Power Spectrum', fontsize=18)

phase_random_data = phase_randomize(sub_rois[sub_1])
a[1,0].plot(phase_random_data, linewidth=2, color='red')
a[1,0].set_ylabel('Avg Activity', fontsize=16)
a[1,0].set_xlabel('Time (TR)', fontsize=16)
a[1,0].set_title('Phase Randomized Data', fontsize=16)

fft_phase = fft(phase_random_data)
a[1,1].plot(freq[:n_freq], np.abs(fft_phase)[:n_freq], linewidth=2, color='red')
a[1,1].set_xlabel('Frequency (Hz)', fontsize=18)
a[1,1].set_ylabel('Amplitude', fontsize=18)
a[1,1].set_title('Phase Randomized Power Spectrum', fontsize=18)

plt.tight_layout()
../_images/8ffb1bc1cd8c14e8bdaad8b42a2097bb9bd08172101c865407e785d8c845b78c.png
stats_phase = isc(sub_rois, method='phase_randomize', n_samples=5000, return_null=True)

print(f"ISC: {stats_phase['isc']:.02}, p = {stats_phase['p']:.03}")
ISC: 0.085, p = 0.0002

As you can see the ISC value is identical as the median of the pairwise correlations are identical. The p-values are also similar and likely reflect the limited precision of the possible p-values that can be computed using only 5,000 permutations. For greater precision, you will need to increase the number of permutations, but remember that this will also increase the computational time.

Subject-wise Bootstrapping#

The final approach we will illustrate is subject-wise bootstrapping of the pairwise similarity matrix. This approach is more conservative than the previously described methods and is almost an order of magnitude faster shuffling the similarity matrix rather than recomputing the pairwise similarity for the null distribution (Chen et al., 2016). Bootstrapping and permutation tests are different types of resampling statistics (see our resampling tutorial for a more in depth overview). Bootstrapping is typically used more for generating confidence intervals around an estimator, while permutation tests are used for performing hypothesis tests. However, p-values can also be computed using a bootstrap by subtracting the ISC from the null distribution and evaluating the percent of samples from the distribution that are smaller than the observed ISC (Hall et al., 1991).

Just like the examples above, we will illustrate what an example bootstrapped similarity matrix looks like. Notice how some subjects are repeatedly resampled, which means that there are multiple values of perfect correlations found off the diagonal. This should be accounted for by using the median summary statistic of the lower triangle. However, both Brainiak and nltools toolboxes convert these values to NaNs to minimize the impact of these outliers on the summary statistic.

def bootstrap_subject_matrix(similarity_matrix, random_state=None):
    '''This function shuffles subjects within a similarity matrix based on recommendation by Chen et al., 2016'''
    
    random_state = check_random_state(random_state)
    n_sub = similarity_matrix.shape[0]
    bootstrap_subject = sorted(random_state.choice(np.arange(n_sub), size=n_sub, replace=True))
    return similarity_matrix[bootstrap_subject, :][:, bootstrap_subject]


similarity = 1 - pairwise_distances(pd.DataFrame(sub_rois).T, metric='correlation')

f,a = plt.subplots(ncols=2, figsize=(12, 6), sharey=True)
sns.heatmap(similarity, square=True, cmap='RdBu_r', vmin=-1, vmax=1, xticklabels=False, yticklabels=False, ax=a[0])
a[0].set_ylabel('Subject', fontsize=18)
a[0].set_xlabel('Subject', fontsize=18)
a[0].set_title('Pairwise Similarity', fontsize=16)

sns.heatmap(bootstrap_subject_matrix(similarity), square=True, cmap='RdBu_r', vmin=-1, vmax=1, xticklabels=False, yticklabels=False, ax=a[1])
a[1].set_ylabel('Subject', fontsize=18)
a[1].set_xlabel('Subject', fontsize=18)
a[1].set_title('Bootstrapped Pairwise Similarity', fontsize=16)
Text(0.5, 1.0, 'Bootstrapped Pairwise Similarity')
../_images/103a4527078732052abc833e336ff74c3aedbfc0fb2a2726bf4664740fd04cd0.png
stats_boot = isc(sub_rois, method='bootstrap', n_samples=5000, return_null=True)

print(f"ISC: {stats_boot['isc']:.02}, p = {stats_boot['p']:.03}")
ISC: 0.085, p = 0.0002

The bootstrap procedure tends to run much faster than the permutation methods on our computer, which is one of the reasons that the authors of the toolbox recommend this approach, beyond it being more conservative (Chen et al., 2016).

Null Distributions#

Since each of our examples saved the null distribution, we can plot a histogram of the null distribution from each method including the confidence intervals. Notice that the circle shift and phase randomization methods produce highly similar null distributions and confidence intervals, while the bootstrap method has a wider and less symmetric distribution with the current number of iterations. However, the observed ISC of 0.085 (red line) exceeds all of the samples from the null distribution yielding a very small p-value. Notice that our observed null distributions using the surrogate data derived from the Sherlock data are very similar to the animation of null data presented above.

plt.figure(figsize=(12,6))
sns.distplot(stats_boot['null_distribution'] - stats_boot['isc'], kde=False, label='Bootstrap')
sns.distplot(stats_circle['null_distribution'], kde=False, label='Bootstrap')
sns.distplot(stats_phase['null_distribution'], kde=False, label='Bootstrap')
plt.ylabel('Frequency', fontsize=18)
plt.xlabel('ISC Values (correlation)', fontsize=18)
plt.title('ISC Null Distribution', fontsize=20)
plt.axvline(stats_boot['isc'], linestyle='-', color='red', linewidth=4)
plt.legend(['Bootstrap Null','Circle Shift Null', 'Phase Randomize Null','Observed ISC'], fontsize=18)

plt.axvline(stats_boot['ci'][0] - stats_boot['isc'], linestyle='--', color='blue')
plt.axvline(stats_boot['ci'][1] - stats_boot['isc'], linestyle='--', color='blue')

plt.axvline(stats_circle['ci'][0], linestyle='--', color='orange')
plt.axvline(stats_circle['ci'][1], linestyle='--', color='orange')

plt.axvline(stats_phase['ci'][0], linestyle='--', color='green')
plt.axvline(stats_phase['ci'][1], linestyle='--', color='green')
<matplotlib.lines.Line2D at 0x16c7b3e20>
../_images/bf7951ddb1d8ba867e163395faa2b5cd0b7d399444a003fbaf1b17dfa1552943.png

Whole brain ISC#

Now, let’s calculate ISC looping over each of the 50 ROIs from the whole-brain meta-analytic parcellation (de la Vega, 2016). Here we loop over each ROI and grab the column from each subject’s dataframe. We then run ISC on the combined subject’s ROI timeseries using the median method and compute a hypothesis test using the subject-wise bootstrap method with 5000 samples (Chen et al., 2016). Finally, we convert each correlation and p-value from each region back into a Brain_Data instance.

isc_r, isc_p = {}, {}
for roi in tqdm(range(len(mask_x))):
    stats = isc(get_subject_roi(sub_timeseries, roi), n_samples=5000, metric='median', method='bootstrap')
    isc_r[roi], isc_p[roi] = stats['isc'], stats['p']
isc_r_brain, isc_p_brain = roi_to_brain(pd.Series(isc_r), mask_x), roi_to_brain(pd.Series(isc_p), mask_x)

Now we can plot the ISC values to visualize which regions had a higher overall ISC.

isc_r_brain.plot(cmap='RdBu_r')
../_images/1cfeb4a02e3c4bcbd42e3079c46f53cfd375abe30dee37277a2a2b57e030d571.png
view_img(isc_r_brain.to_nifti())

We can threshold using bonferroni correction (p < 0.001 for k=50 parcellation). Alternatively, we can threshold using false discovery rate, by setting thr=fdr(isc_p_brain.data). In this example, FDR is more conservative than bonferroni.

view_img_on_surf(threshold(isc_r_brain, isc_p_brain, thr=.001).to_nifti())