isc_tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "view-in-github"
   },
   "source": [
    "<a href=\"https://colab.research.google.com/github/snastase/isc-tutorial/blob/master/isc_tutorial.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "fW_aiFKcIF1b"
   },
   "source": [
    "<h1>Intersubject correlation (ISC) tutorial</h1>\n",
    "This tutorial jupyter notebook accompanies the manuscript \"Measuring shared responses across subjects using intersubject correlation\" by Nastase, Gazzola, Hasson, and Keysers. The goal of the tutorial is to introduce basic intersubject correlation (ISC) analyses ([Hasson et al., 2004](https://doi.org/10.1126/science.1089506), [2010](https://doi.org/10.1016/j.tics.2009.10.011)) and subsequent statistical tests as implemented in Python using the Brain Imaging Analysis Kit ([BrainIAK](http://brainiak.org/)). Click *Open in playground* to interactively run and edit cells (you may need to sign into a Google account), and use *File* > *Save a copy in Drive...* or *Save a copy in GitHub...* to save your changes. To execute a code cell, click in that cell and press _Shift_ + _Enter_ (or _Shift_ + _Return_ on a Mac). The first time you run a code cell in playground mode, you may receive a warning and a prompt to reset runtimes; click _Run anyway_ followed by _Yes_.\n",
    "\n",
    "---\n",
    "\n",
    "Author: Samuel A. Nastase"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "_vlNUGYOPmce"
   },
   "source": [
    "## Getting started\n",
    "First, we'll need to [install BrainIAK](http://brainiak.org/docs/installation.html) and its requirements in the cloud instance—this may take a few minutes. For this tutorial, we'll install BrainIAK directly from the [GitHub repository](https://github.com/brainiak/brainiak) to get the most recent features. Note that if you are running this notebook locally, the installation process may vary and you should not execute this code cell. If you want to install BrainIAK locally, run `pip install git+https://github.com/snastase/brainiak.git` in a separate code cell or follow the online [installation instructions](http://brainiak.org/docs/installation.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "dZrDDlk8O27m"
   },
   "outputs": [],
   "source": [
    "# Install BrainIAK requirements in Google Colab Linux cloud instance\n",
    "!apt install build-essential libgomp1 libmpich-dev mpich python3-dev \\\n",
    "             python3-pip python3-venv\n",
    "\n",
    "# Install most recent BrainIAK from GitHub in the cloud instance\n",
    "!pip install pip==9.0.1\n",
    "!pip install git+https://github.com/brainiak/brainiak.git"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "BOfJ_JKy4hbX"
   },
   "source": [
    "Next, we'll import the relevant functions from BrainIAK"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "qIcccVRr22vX"
   },
   "outputs": [],
   "source": [
    "from brainiak.isc import (isc, isfc, bootstrap_isc, permutation_isc,\n",
    "                          timeshift_isc, phaseshift_isc,\n",
    "                          compute_summary_statistic)\n",
    "from brainiak.io import load_boolean_mask, load_images\n",
    "from brainiak.image import mask_images, MaskedMultiSubjectData"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "cwtApadEfo9F"
   },
   "source": [
    "If unable to install BrainIAK, we can use basic ISC functionality without the full BrainIAK package. We'll download `isc_standalone.py` from the [GitHub repository](https://github.com/snastase/isc-tutorial) for this tutorial and load the necessary modules locally. If you've cloned the `isc-tutorial` GitHub repository locally, this step is not necessary (as the local directory for the repository already contains `isc_standalone.py`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "zbPTnYXeIm_m"
   },
   "outputs": [],
   "source": [
    "# Download the standalone module if not using BrainIAK\n",
    "from urllib.request import urlretrieve\n",
    "urlretrieve('https://github.com/snastase/isc-tutorial/'\n",
    "            'raw/master/isc_standalone.py', 'isc_standalone.py');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "G9c3pQSZq-Kq"
   },
   "source": [
    "If you're using the `isc_standalone.py` module instead of BrainIAK, import the relevant functions from `isc_standalone`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "mTYTe5Kaq5-t"
   },
   "outputs": [],
   "source": [
    "from isc_standalone import (isc, isfc, bootstrap_isc, permutation_isc,\n",
    "                            timeshift_isc, phaseshift_isc,\n",
    "                            compute_summary_statistic, load_images,\n",
    "                            load_boolean_mask, mask_images,\n",
    "                            MaskedMultiSubjectData)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "-OEKROQZK64a"
   },
   "source": [
    "Finally, we'll load several other useful Python modules."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "9VRpzXF0K3ll"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from scipy.stats import norm, pearsonr, zscore\n",
    "from scipy.spatial.distance import squareform\n",
    "from statsmodels.stats.multitest import multipletests\n",
    "import nibabel as nib"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "ZkZ2xqvfHBhC"
   },
   "source": [
    "## Example data\n",
    "We'll create a simple simulated dataset for quickly applying ISC analyses, then later apply the analyses to a real fMRI dataset where participants listened to a spoken narrative ([Pie Man](https://themoth.org/stories/pie-man) by Jim O'Grady). Our simulated data will have 1,000 voxels in total comprising 10 \"networks\" and 300 time points (or TRs)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "Lvt3Fk_3Osjx"
   },
   "outputs": [],
   "source": [
    "# Set parameters for toy time series data\n",
    "n_subjects = 20\n",
    "n_TRs = 300\n",
    "n_voxels = 1000\n",
    "\n",
    "# Create simple simulated data with high intersubject correlation\n",
    "def simulated_timeseries(n_subjects, n_TRs, n_voxels=1, noise=1):\n",
    "  signal = np.random.randn(n_TRs, n_voxels // 100)\n",
    "  data = [zscore(np.repeat(signal, 100, axis=1) +\n",
    "                 np.random.randn(n_TRs, n_voxels) * noise,\n",
    "                 axis=0)\n",
    "          for subject in np.arange(n_subjects)]\n",
    "  return data\n",
    "\n",
    "# List of subject datasets\n",
    "data = simulated_timeseries(n_subjects, n_TRs, n_voxels=n_voxels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "AxGPkSaLhbP_"
   },
   "outputs": [],
   "source": [
    "# Inspect the shape of one of our simulated datasets\n",
    "print(f\"Simulated data shape first subject: {data[0].shape} \"\n",
    "      f\"\\ni.e., {data[0].shape[0]} time points and {data[0].shape[1]} voxels\")\n",
    "\n",
    "# Create a simple visualization of the data\n",
    "plt.matshow(data[0], cmap='RdYlBu_r', vmin=-3, vmax=3)\n",
    "plt.grid(False)\n",
    "plt.xlabel('voxels')\n",
    "plt.ylabel('time points');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Q9yhDQYxKJc1"
   },
   "source": [
    "## ISC analysis\n",
    "Let's start very simple by computing the ISC for a single voxel (or ROI) across only two participants. This should give us a simple Pearson correlation value (and should match other implementations of Pearson correlation). Note that when you call the `isc` function with `verbose=True` (the default), it outputs some warnings describing what it infers about the input data. If these don't match your assumptions, your input data may be organized improperly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "L4K_t2VxOvmN"
   },
   "outputs": [],
   "source": [
    "# Get the time series for a single voxel in two subjects\n",
    "subject_a = data[0][:, 0]\n",
    "subject_b = data[1][:, 0]\n",
    "\n",
    "# Check the shape of these mini-datasets\n",
    "print(f\"Subject A, first voxel, shape = {subject_a.shape} \"\n",
    "      f\"\\nSubject B, first voxel, shape = {subject_b.shape}\")\n",
    "\n",
    "# Combine these into a list\n",
    "both_subjects = [subject_a, subject_b]\n",
    "\n",
    "# Compute the ISC for this voxel across the two subjects\n",
    "iscs = isc(both_subjects, pairwise=True)\n",
    "print(f\"ISC for first voxel across subjects A and B = {iscs[0]}\")\n",
    "\n",
    "# NB: even for a single voxel, the output ISC is shaped to \n",
    "# to accommodate an n_ISCs x n_voxels matrix\n",
    "print(f\"ISC output shape = {iscs.shape}\"\n",
    "      f\"\\ni.e., {iscs.shape[0]} ISC value(s) by {iscs.shape[0]} voxel(s)\")\n",
    "\n",
    "# Check that ISC output matches of other correlation functions in python\n",
    "numpy_corrcoef = np.corrcoef(subject_a, subject_b)[0, 1]\n",
    "\n",
    "scipy_pearsonr = pearsonr(subject_a, subject_b)[0]\n",
    "\n",
    "print(f\"BrainIAK ISC = {iscs[0]:.6f}\"\n",
    "      f\"\\nNumpy's correlation = {numpy_corrcoef:.6f}\"\n",
    "      f\"\\nScipy's correlation = {scipy_pearsonr:.6f}\")\n",
    "assert np.isclose(iscs, numpy_corrcoef) and np.isclose(iscs, scipy_pearsonr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "zFXMgzm4rBHF"
   },
   "source": [
    "BrainIAK uses Python's logging functionality. To see non-critical messages while running ISC analyses, we temporarily can set the logging level to 'INFO'."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "HKOydSmCrUWH"
   },
   "outputs": [],
   "source": [
    "# Import logging module and set level to INFO\n",
    "import logging\n",
    "logging.basicConfig()\n",
    "logging.getLogger().setLevel(logging.INFO)\n",
    "\n",
    "# Re-run the previous ISC analyses to see logged info\n",
    "iscs = isc(both_subjects, pairwise=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "eJ35OZZXdNPv"
   },
   "outputs": [],
   "source": [
    "# Set logging back to default level of WARNING\n",
    "logging.getLogger().setLevel(logging.WARNING)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "OUIwWPb_1eiJ"
   },
   "source": [
    "When there are three or more subjects, we can compute ISCs using either the pairwise approach (`pairwise=True`), where we compute ISCs between each pair of subjects, or the leave-one-out (`pairwise=False`) approach, where we compute ISCs between each subject and the average time series of other subjects."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "k92U8CYz9ylK"
   },
   "source": [
    "### Pairwise approach\n",
    "Now we'll run the full-scale ISC analysis across all voxels and subjects using the pairwise approach. For a given voxel, the correlations between each pair of subjects are represented in a vector of length\n",
    "```\n",
    "n_subjects * (n_subjects - 1) / 2\n",
    "```\n",
    "or 190 pairs for 20 subjects. This vector of pairs corresponds to the off-diagonal values of a symmetric subjects-by-subjects correlation matrix.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "83ZNUVljnRe-"
   },
   "outputs": [],
   "source": [
    "# Pairwise approach across all subjects and voxels\n",
    "iscs = isc(data, pairwise=True)\n",
    "\n",
    "# Check shape of output ISC values\n",
    "print(f\"ISC values shape = {iscs.shape} \\ni.e., {iscs.shape[0]} \"\n",
    "      f\"pairs and {iscs.shape[1]} voxels\"\n",
    "      f\"\\nMinimum ISC = {np.amin(iscs):.3f}; \"\n",
    "      f\"maximum ISC = {np.amax(iscs):.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "S_SnbdKMHjyI"
   },
   "source": [
    "For a given voxel, we can convert the vector of pairs to the full correlation matrix for visualization. In the simulated dataset, all subjects were designed to have high ISCs; however, we can add noise to some of the subjects and then visualize the ISC matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "epCVUtkWHtW4"
   },
   "outputs": [],
   "source": [
    "# Visualize the correlation matrix for one voxel\n",
    "isc_matrix = squareform(iscs[:, 0])\n",
    "np.fill_diagonal(isc_matrix, 1)\n",
    "sns.heatmap(isc_matrix, cmap=\"RdYlBu_r\", vmin=-1, vmax=1, square=True, \n",
    "            xticklabels=range(1, 21), yticklabels=range(1, 21))\n",
    "plt.xlabel('subjects')\n",
    "plt.ylabel('subjects')\n",
    "plt.show()\n",
    "\n",
    "# Create noisier data\n",
    "noisy_data = np.dstack((np.dstack((\n",
    "    simulated_timeseries(n_subjects // 2, n_TRs,\n",
    "                         n_voxels=n_voxels, noise=1))),\n",
    "                        np.dstack((\n",
    "    simulated_timeseries(n_subjects // 2, n_TRs,\n",
    "                         n_voxels=n_voxels, noise=5)))))\n",
    "\n",
    "# Recompute ISC and visualize data with noisy subjects\n",
    "noisy_iscs = isc(noisy_data, pairwise=True)\n",
    "isc_matrix = squareform(noisy_iscs[:, 0])\n",
    "np.fill_diagonal(isc_matrix, 1)\n",
    "sns.heatmap(isc_matrix, cmap=\"RdYlBu_r\", vmin=-1, vmax=1, square=True, \n",
    "            xticklabels=range(1, 21), yticklabels=range(1, 21))\n",
    "plt.xlabel('subjects')\n",
    "plt.ylabel('subjects')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "-6w0xFXwHOGK"
   },
   "source": [
    "### Leave-one-out approach\n",
    "Instead of computing ISCs between each pair of subjects, for each subject we can compute the ISC between that subject and the average of all other subjects. Notice that the observed ISC values are typically higher in the leave-one-out approach due to computing correlations between the left-out subject and the cleaner averaged time series from the remaining subjects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "gv9CrA12l9Fs"
   },
   "outputs": [],
   "source": [
    "# Leave-one-out approach\n",
    "iscs = isc(data, pairwise=False)\n",
    "\n",
    "# Check shape of output ISC values\n",
    "print(f\"ISC values shape = {iscs.shape} \\ni.e., {iscs.shape[0]} \"\n",
    "      f\"left-out subjects and {iscs.shape[1]} voxel(s)\"\n",
    "      f\"\\nMinimum ISC = {np.amin(iscs):.3f}; \"\n",
    "      f\"maximum ISC = {np.amax(iscs):.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "MeTDDaC3QVIA"
   },
   "source": [
    "### Input types\n",
    "Currently, we're a submitting a list of numpy arrays to BrainIAK's `isc` function where each item in the list is a subject's response time course over some number of voxels. Alternatively, we could stack subjects along the 3rd dimension (`np.dstack`) into a single 3-dimensional numpy array and submit this to the `isc` function. If the `isc` function receives a single numpy array, it will assume that the last dimension indexes subjects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "1mcYGPCGQR3T"
   },
   "outputs": [],
   "source": [
    "# Input a list of subjects (same as before)\n",
    "iscs = isc(data, pairwise=False)\n",
    "\n",
    "# Stack subjects in 3rd-dimension and recompute ISC\n",
    "data_stack = np.dstack(data)\n",
    "print(f\"Stacked data shape = {data_stack.shape}\"\n",
    "      f\"\\ni.e., {data_stack.shape[0]} time points, {data_stack.shape[1]} \"\n",
    "      f\"voxels, and {data_stack.shape[2]} subjects\")\n",
    "\n",
    "# Input stacked numpy array\n",
    "iscs_from_stack = isc(data_stack, pairwise=False)\n",
    "\n",
    "# Make sure the ISC outputs are the same\n",
    "assert np.array_equal(iscs, iscs_from_stack)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "jX7-4nkAVyaV"
   },
   "source": [
    "### Summary statistics\n",
    "Rather than returning ISC values for each pair of subject (in the pairwise approach) or each left-out subject (in the leave-one-out approach), we can use the `summary_statistic` argument to output either the mean or median across the values. Note that by default `summary_statistic=False`. If we request the mean ISC value, the `isc` function will internally apply the Fisher *z*-transformation (`np.arctanh`) prior to computing the mean, then apply the inverse Fisher *z*-transformation (`np.tanh`) to the mean value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "eSpl72nsW82m"
   },
   "outputs": [],
   "source": [
    "# Compute mean leave-one-out ISC\n",
    "iscs = isc(data, pairwise=False, summary_statistic='mean')\n",
    "\n",
    "print(f\"ISC values shape = {iscs.shape} \\ni.e., the mean value across \"\n",
    "      f\"left-out subjects for {iscs.shape[0]} voxel(s)\"\n",
    "      f\"\\nMean ISC for first voxel = {iscs[0]:.3f}\")\n",
    "\n",
    "# Compute median leave-one-out ISC\n",
    "iscs = isc(data, pairwise=False, summary_statistic='median')\n",
    "\n",
    "print(f\"ISC values shape = {iscs.shape} \\ni.e., the median value across \"\n",
    "      f\"left-out subjects for {iscs.shape[0]} voxel(s)\"\n",
    "      f\"\\nMedian ISC for first voxel = {iscs[0]:.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "M980ge4ILSRk"
   },
   "source": [
    "## Statistical tests\n",
    "BrainIAK provides several nonparametric statistical tests for ISC analysis. Nonparametric tests are preferred due to the inherent correlation structure across ISC values—each subject contributes to the ISC of other subjects, violating assumptions of independence required for standard parametric tests (e.g., *t*-test, ANOVA). The nonparametric statistical tests discussed below return the actual observed ISC values, *p*-values, and the resampling distribution (the bootstrap hypothesis test also returns confidence intervals around the observed ISC statistic). For expediency, we only use 200 resampling iterations here, but 1,000 or more iterations are generally recommended."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "HfjWbhcGbBvu"
   },
   "source": [
    "### Phase randomization\n",
    "One approach for statistically assessing ISCs is to randomize the phase of time series across subjects prior to computing ISCs (e.g., [Lerner et al., 2011](https://doi.org/10.1523/jneurosci.3684-10.2011); [Simony et al., 2016](https://doi.org/10.1038/ncomms12141)). This method requires recomputing ISC at each iteration of the randomization test, and is therefore slow. In the pairwise approach, we phase-randomize each subject prior to computing ISCs; however, in the leave-one-out approach, we only phase-randomize the left-out subject prior to computing ISC. At each iteration of the phase randomization test, the same random phase shift is used across all voxels to preserve the spatial autocorrelation of typical fMRI data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "wqLNWogsbwu6"
   },
   "outputs": [],
   "source": [
    "# Phase randomization using pairwise approach (takes a couple minutes)\n",
    "observed, p, distribution = phaseshift_isc(data, pairwise=True,\n",
    "                                           summary_statistic='median',\n",
    "                                           n_shifts=200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "Y2IFdfPLe3rJ"
   },
   "outputs": [],
   "source": [
    "# Inspect shape of null distribution\n",
    "print(f\"Null distribution shape = {distribution.shape}\"\n",
    "      f\"\\ni.e., {distribution.shape[0]} randomizations \"\n",
    "      f\"and {distribution.shape[1]} voxels\")\n",
    "\n",
    "# Get actual ISC value and p-value for first voxel\n",
    "print(f\"Actual observed ISC value for first voxel = {observed[0]:.3f},\"\n",
    "      f\"\\np-value from randomization test = {p[0]:.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Cz7upBYscQNw"
   },
   "source": [
    "### Circular time-shift randomization\n",
    "A conceptually similar nonparametric approach is to circularly shift the response time series across subjects by random offsets ([Kauppi et al., 2014](https://doi.org/10.3389/fninf.2014.00002)). Time points that would be shifted beyond the end of the time series are wrapped around to the beginning of the time series."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "NIHrG5XucsWL"
   },
   "outputs": [],
   "source": [
    "# Circular time-shift using pairwise approach (takes a couple minutes)\n",
    "observed, p, distribution = timeshift_isc(data, pairwise=True,\n",
    "                                          summary_statistic='median',\n",
    "                                          n_shifts=200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "LXNQl6IbudTj"
   },
   "outputs": [],
   "source": [
    "# Inspect shape of null distribution\n",
    "print(f\"Null distribution shape = {distribution.shape}\"\n",
    "      f\"\\ni.e., {distribution.shape[0]} randomizations \"\n",
    "      f\"and {distribution.shape[1]} voxels\")\n",
    "\n",
    "# Get actual ISC value and p-value for first voxel\n",
    "print(f\"Actual observed ISC value for first voxel = {observed[0]:.3f},\"\n",
    "      f\"\\np-value from randomization test = {p[0]:.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "dxWZVZNpdD4C"
   },
   "source": [
    "### Bootstrap hypothesis test\n",
    "We can also perform group-level statistical tests that operate directly on the observed ISC values and do not require recomputing ISCs. For one-sample tests, we can resample subjects with replacement to construct a bootstrap distribution around our observed ISC statistic ([Chen et al., 2016](https://doi.org/10.1016/j.neuroimage.2016.05.023)). We can compute confidence intervals around the test statistic using the `ci_percentile` option (default 95%). Hypothesis test is performed by shifting the bootstrap distribution to zero. Note that when constructing the bootstrap distribution using the pairwise approach, subjects (i.e., rows and columns in the subject-by-subject correlation matrix) are sampled with replacement, not pairs (which would disrupt the correlation structure among pairs)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "M71PFb1Cb1Gd"
   },
   "outputs": [],
   "source": [
    "# Compute ISCs and then run bootstrap hypothesis test on ISCs\n",
    "iscs = isc(data, pairwise=True, summary_statistic=None)\n",
    "observed, ci, p, distribution = bootstrap_isc(iscs, pairwise=True,\n",
    "                                              ci_percentile=95,\n",
    "                                              summary_statistic='median',\n",
    "                                              n_bootstraps=200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "_uyla8jvv6uv"
   },
   "outputs": [],
   "source": [
    "# Inspect shape of null distribution\n",
    "print(f\"Null distribution shape = {distribution.shape}\"\n",
    "      f\"\\ni.e., {distribution.shape[0]} bootstraps \"\n",
    "      f\"and {distribution.shape[1]} voxels\")\n",
    "\n",
    "# Get actual ISC value and p-value for first voxel\n",
    "print(f\"Actual observed ISC value for first voxel = {observed[0]:.3f},\"\n",
    "      f\"\\np-value from bootstrap hypothesis test = {p[0]:.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "XCj_y0rGybFH"
   },
   "source": [
    "### Permutation test\n",
    "We can use a permutation test to statistically evaluate one- or two-sample tests ([Chen et al., 2016](https://doi.org/10.1016/j.neuroimage.2016.05.023)). In the case of a one-sample test, we use a sign-flipping (-1, +1) approach applied to the observed ISC. For a two-sample test, we supply a `group_assignment` list containing the group labels for each subject. The order of the group assignment list must match the order in which the subjects are supplied to the `isc` function. At each iteration, we randomly reassign the group labels, then compute the test statistic. In the one-sample test, there are `2**n_subjects` number of possible permutations, while in the the two-sample test, there are `n_subjects!` number of possible permutations. In both cases, if the requested number of permutations equals or exceeds the exhaustive list of permutations, an exact test is performed using all possible permutations. However, in most cases the number of subjects will yield an a prohibitively large number of permutations, in which case a Monte Carlo approximate permutation test is used instead of an exact test."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "IRDEs9bFLmSq"
   },
   "outputs": [],
   "source": [
    "# Compute ISCs and then run one-sample permutation test on ISCs\n",
    "iscs = isc(data, pairwise=True, summary_statistic=None)\n",
    "observed, p, distribution = permutation_isc(iscs, pairwise=True,\n",
    "                                            summary_statistic='median',\n",
    "                                            n_permutations=200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "yWCy7RW6vefV"
   },
   "outputs": [],
   "source": [
    "# Inspect shape of null distribution\n",
    "print(f\"Null distribution shape = {distribution.shape}\"\n",
    "      f\"\\ni.e., {distribution.shape[0]} permutations \"\n",
    "      f\"and {distribution.shape[1]} voxels\")\n",
    "\n",
    "# Get actual ISC value and p-value for first voxel\n",
    "print(f\"Actual observed ISC value for first voxel = {observed[0][0]:.3f},\"\n",
    "      f\"\\np-value from permutation test = {p[0]:.3f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "8pZSNtOW32E_"
   },
   "outputs": [],
   "source": [
    "# Note that with few subjects, an exact test is performed\n",
    "data_n6 = data[:6]\n",
    "iscs = isc(data_n6, pairwise=True, summary_statistic=None)\n",
    "observed, p, distribution = permutation_isc(iscs, pairwise=True,\n",
    "                                            summary_statistic='median',\n",
    "                                            n_permutations=200)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "vEyVRgV84kXR"
   },
   "source": [
    "If we have two groups we expect to have different ISC values, we must supply a `group_assignment` list. In the case of two groups, we compute the difference between the `summary_statistic` for each group. In the pairwise approach, we compute differences between the `summary_statistic` for within-group correlations, and ignore the between group correlations in the full subject-by-subject correlation matrix containing both groups. Furthermore, permutations are applied to subjects (i.e., rows and columns in the subject-by-subject correlation matrix) and not to pairs (which would disrupt the correlation structure among pairs). We'll construct a dataset where one group of subjects are noisier than the others."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "bmriCr0b5UQg"
   },
   "outputs": [],
   "source": [
    "# Create data with noisy subset of subjects\n",
    "noisy_data = np.dstack((np.dstack((\n",
    "    simulated_timeseries(n_subjects // 2, n_TRs,\n",
    "                         n_voxels=n_voxels, noise=1))),\n",
    "                        np.dstack((\n",
    "    simulated_timeseries(n_subjects // 2, n_TRs,\n",
    "                         n_voxels=n_voxels, noise=5)))))\n",
    "\n",
    "# Create group_assignment variable with group labels\n",
    "group_assignment = [1]*10 + [2]*10\n",
    "print(f\"Group assignments: \\n{group_assignment}\")\n",
    "\n",
    "# Compute ISCs and then run two-sample permutation test on ISCs\n",
    "iscs = isc(noisy_data, pairwise=True, summary_statistic=None)\n",
    "observed, p, distribution = permutation_isc(iscs,\n",
    "                                            group_assignment=group_assignment,\n",
    "                                            pairwise=True,\n",
    "                                            summary_statistic='median',\n",
    "                                            n_permutations=200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "BtHnvJ026l53"
   },
   "outputs": [],
   "source": [
    "# Inspect shape of null distribution\n",
    "print(f\"Null distribution shape = {distribution.shape}\"\n",
    "      f\"\\ni.e., {distribution.shape[0]} permutations \"\n",
    "      f\"and {distribution.shape[1]} voxels\")\n",
    "\n",
    "# Get actual ISC value and p-value for first voxel\n",
    "print(f\"Actual observed group difference in ISC values \"\n",
    "      f\"for first voxel = {observed[0]:.3f},\"\n",
    "      f\"\\np-value from permutation test = {p[0]:.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "IW7-B0A6MsX_"
   },
   "source": [
    "## Correcting for multiple tests\n",
    "Evaluating the statistical significance of an ISC analysis across many voxels will result in many false positives unless we somehow control for the large number of statistical tests. Here we'll use two simple approaches for correcting for multiple tests. In the first approach, we'll account for multiple tests by controlling the expected proportion of false positives or false discovery rate (FDR; [Benjamini & Hochberg, 1995](https://www.jstor.org/stable/2346101); [Benjamini & Yekutieli, 2001](https://www.jstor.org/stable/2674075); [Genovese et al., 2002](https://doi.org/10.1006/nimg.2001.1037)). In the second approach, we'll control the family-wise error rate (FWER) by constructing a null distribution from the maximum ISC value across all voxels at each iteration of a randomization test ([Nichols & Holmes, 2002](https://doi.org/10.1002/hbm.1058))."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Dn7s0Rq1hASj"
   },
   "source": [
    "First, we'll create a dataset where half the voxels are very consistent across subjects and the other half are very noisy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "6thjueBb_fQ3"
   },
   "outputs": [],
   "source": [
    "# Create data with where half of voxels are noisy\n",
    "noisy_data = np.hstack((np.dstack((\n",
    "    simulated_timeseries(n_subjects, n_TRs,\n",
    "                         n_voxels=n_voxels // 2, noise=1))),\n",
    "                        np.dstack((\n",
    "    simulated_timeseries(n_subjects, n_TRs,\n",
    "                         n_voxels=n_voxels // 2, noise=9)))))\n",
    "\n",
    "# Visualize data for first subject where half of voxels are noisy\n",
    "plt.matshow(noisy_data[..., 0], cmap='RdYlBu_r', vmin=-3, vmax=3)\n",
    "plt.grid(False)\n",
    "plt.xlabel('voxels')\n",
    "plt.ylabel('time points');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "mhQsh3fRi324"
   },
   "source": [
    "After visualizing these data, we'll compute ISCs and use a one-sample two-sided bootstrap hypothesis test, which yields *p*-values and a null distribution. For this example, we'll use a more realistic number of bootstrap samples (1,000)—this may take a couple minutes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "3nUC8PN0_I58"
   },
   "outputs": [],
   "source": [
    "# Compute ISCs and then run bootstrap hypothesis test on ISCs\n",
    "# using a realistic number of permutations (takes a few minutes)\n",
    "iscs = isc(noisy_data, pairwise=True, summary_statistic=None)\n",
    "observed, ci, p, distribution = bootstrap_isc(iscs, pairwise=True,\n",
    "                                              ci_percentile=95,\n",
    "                                              summary_statistic='median',\n",
    "                                              n_bootstraps=1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "73nFdL7EiArS"
   },
   "source": [
    "### Controlling FDR\n",
    "To control FDR, we'll use the the `multipletests` function from the StatsModels Python package. This returns an array of *q*-values, which are typically interpreted as FDR-corrected *p*-values. By thresholding uncorrected and corrected *p*- and *q*-values, we can determine how many voxels survived correction for multiple tests."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "0hTAf46V-PtG"
   },
   "outputs": [],
   "source": [
    "# Get q-values (i.e., FDR-controlled p-values) using statsmodels\n",
    "q = multipletests(p, method='fdr_by')[1]\n",
    "\n",
    "# We can also convert these q-values to z-values\n",
    "z = np.abs(norm.ppf(q))\n",
    "\n",
    "# Also get significant voxels with and without correction\n",
    "corrected = q[np.newaxis, :] < .05\n",
    "uncorrected = p[np.newaxis, :] < .05\n",
    "\n",
    "# Count significant voxels before and after correction\n",
    "print(f'{np.sum(uncorrected)} \"significant\" voxels before correction for '\n",
    "      f\"multiple tests; {np.sum(corrected)} significant voxels after \"\n",
    "      f\"controlling FDR at .05\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "KLjGo2Nvjrs6"
   },
   "source": [
    "Finally, we can visualize the voxel time series for an example subject, the ISC values across subjects, and which voxels are considered significant before and after controlling FDR at .05. Note that before correction even some of the noisy voxels are considered to have significant ISC; however, after correction, the number of significant noisy voxels is reduced."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "D_CVjOeUMwfH"
   },
   "outputs": [],
   "source": [
    "# Set up grid of subplots for visualizing voxel values and significance\n",
    "fig, (ax0, ax1, ax2, ax3) = plt.subplots(nrows=4, figsize=(12, 8),\n",
    "                                         sharex=True,\n",
    "                                         gridspec_kw={'height_ratios':\n",
    "                                                      [300, 190, 20, 20]})\n",
    "\n",
    "# Visualize data for first subject where half of voxels are noisy\n",
    "ax0.matshow(noisy_data[..., 0], cmap='RdYlBu_r', vmin=-3, vmax=3)\n",
    "ax0.grid(False)\n",
    "ax0.set_ylabel('time points')\n",
    "ax0.set_title('response time series for example subject', y=1)\n",
    "\n",
    "# Visualize ISC values across all pairs of subjects\n",
    "ax1.matshow(iscs, cmap='RdYlBu_r', vmin=-1, vmax=1)\n",
    "ax1.grid(False)\n",
    "ax1.set_ylabel('pairs of subjects')\n",
    "ax1.set_title('ISC values for all pairs of subjects', y=1)\n",
    "\n",
    "# Visualize uncorrected and corrected significant voxels\n",
    "ax2.matshow(np.repeat(uncorrected, 20, axis=0),\n",
    "            cmap='viridis',vmin=0, vmax=1)\n",
    "ax2.grid(False)\n",
    "ax2.set_yticks([])\n",
    "ax2.set_title('uncorrrected \"significant\" voxels (yellow)')\n",
    "\n",
    "ax3.matshow(np.repeat(corrected, 20, axis=0),\n",
    "            cmap='viridis',vmin=0, vmax=1)\n",
    "ax3.grid(False)\n",
    "ax3.set_xlabel('voxels')\n",
    "ax3.xaxis.tick_bottom()\n",
    "ax3.set_yticks([])\n",
    "ax3.set_title('FDR-corrrected significant voxels (yellow)')\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "54SD2b00lCS7"
   },
   "source": [
    "### Controlling FWER\n",
    "To strictly control the FWER, one method is to construct a null distribution of maximum ISC statistics across all voxels. First we'll use the `permutation_isc` function to run a one-sample two-sided permutation test using a sign-flipping procedure, which returns *p*-values and a null distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "I_I9WZrdoaJq"
   },
   "outputs": [],
   "source": [
    "# Compute ISCs and then run two-sample permutation test on ISCs\n",
    "iscs = isc(noisy_data, pairwise=True, summary_statistic=None)\n",
    "observed, p, distribution = permutation_isc(iscs, pairwise=True,\n",
    "                                            summary_statistic='mean',\n",
    "                                            n_permutations=1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "hKndHU0upTmB"
   },
   "source": [
    "Next, we'll write a simple function that takes a null distribution with multiple voxels, and aggregates the maximum ISC value across all voxels for each null sample."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "VzEhE_2ulpXw"
   },
   "outputs": [],
   "source": [
    "# Loop through null distribution and get maximum value across voxels\n",
    "def get_maxima(distribution):\n",
    "  max_distribution = []\n",
    "  for i in distribution:\n",
    "    max_isc = np.amax(i)\n",
    "    max_distribution.append(max_isc)\n",
    "  max_distribution = np.array(max_distribution)\n",
    "  return max_distribution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "5MMF8fNQpoDt"
   },
   "source": [
    "After we create a null distribution of maximum statistics, any voxel with an ISC value in the top 5% of distribution can be considered significant. Here we compute *p*-values from the null distribution of maximum statistics using a two-sided test."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "n6Ph-EV-qJG_"
   },
   "outputs": [],
   "source": [
    "# Create null distribution of maximum ISCs across all voxels\n",
    "max_distribution = get_maxima(distribution)\n",
    "\n",
    "# Broadcast our max distribution across all 1000 voxels\n",
    "max_distribution = np.repeat(max_distribution[:, np.newaxis], 1000, axis=1)\n",
    "\n",
    "# Get the summary statistic (median) for our actual ISC values\n",
    "# since we set summary_statistic=None above\n",
    "observed = np.median(iscs, axis=0)[np.newaxis, :]\n",
    "\n",
    "# Evaluate whether observed ISCs land in the tail of the max distribution\n",
    "p_max = ((np.sum(np.abs(max_distribution) >= np.abs(observed), axis=0) + 1) /\n",
    "          float((len(max_distribution) + 1)))[np.newaxis, :]\n",
    "\n",
    "# Get p-values less than .05 (corrected for multiple tests)\n",
    "corrected = p_max < .05"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "b0dZ2nSAuWFU"
   },
   "source": [
    "As with the FDR approach,  we can visualize the data and the voxels marked as significant before and after correction for multiple tests. This method of correction for multiple tests is considerably more conservative."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "-dY6rNoOnrz6"
   },
   "outputs": [],
   "source": [
    "# Set up grid of subplots for visualizing voxel values and significance\n",
    "fig, (ax0, ax1, ax2, ax3) = plt.subplots(nrows=4, figsize=(12, 8),\n",
    "                                         sharex=True,\n",
    "                                         gridspec_kw={'height_ratios':\n",
    "                                                      [300, 190, 20, 20]})\n",
    "\n",
    "# Visualize data for first subject where half of voxels are noisy\n",
    "ax0.matshow(noisy_data[..., 0], cmap='RdYlBu_r', vmin=-3, vmax=3)\n",
    "ax0.grid(False)\n",
    "ax0.set_ylabel('time points')\n",
    "ax0.set_title('response time series for example subject', y=1)\n",
    "\n",
    "# Visualize ISC values across all pairs of subjects\n",
    "ax1.matshow(iscs, cmap='RdYlBu_r', vmin=-1, vmax=1)\n",
    "ax1.grid(False)\n",
    "ax1.set_ylabel('pairs of subjects')\n",
    "ax1.set_title('ISC values for all pairs of subjects', y=1)\n",
    "\n",
    "# Visualize uncorrected and corrected significant voxels\n",
    "ax2.matshow(np.repeat(uncorrected, 20, axis=0),\n",
    "            cmap='viridis',vmin=0, vmax=1)\n",
    "ax2.grid(False)\n",
    "ax2.set_yticks([])\n",
    "ax2.set_title('uncorrrected \"significant\" voxels (yellow)')\n",
    "\n",
    "ax3.matshow(np.repeat(corrected, 20, axis=0),\n",
    "            cmap='viridis',vmin=0, vmax=1)\n",
    "ax3.grid(False)\n",
    "ax3.set_xlabel('voxels')\n",
    "ax3.xaxis.tick_bottom()\n",
    "ax3.set_yticks([])\n",
    "ax3.set_title('FWER-corrrected significant voxels (yellow)')\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "5ZxUKhmDkw35"
   },
   "source": [
    "Note that there are many other ways to correct for multiple tests, such as using cluster-extent thresholding; but these methods are beyond the scope of this tutorial."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "ZW7iXQvPOtiy"
   },
   "source": [
    "## ISFC analysis\n",
    "Rather than computing ISCs for corresponding voxels across participants, we can instead compute ISCs between all voxels to measure functional integration (i.e., connectivity). This method is called intersubject functional correlation (ISFC) analysis ([Simony et al., 2016](https://doi.org/10.1038/ncomms12141)). Using the `vectorize_isfcs` option, we can either return a tuple containing the condensed off-diagonal ISFC values and the diagonal ISC values or the square (redundant) ISFC values. If `vectorize_isfcs=True` (the default), the first array in the tuple contains the off-diagonal ISFC values for each pair of voxels as condensed by `scipy.spatial.distance.squareform` and is shaped `n_subjects` (or `n_pairs`) by `n_connections` where\n",
    "```\n",
    "n_connections = n_voxels * (n_voxels - 1) / 2\n",
    "```\n",
    "The second array in the tuple is the diagonal values shaped `n_subjects` (or `n_pairs`) by `n_voxels`. If `vectorize_isfcs=False`, we get a 3-dimensional array containing the square (redundant) ISFC and ISC values, shaped `n_subjects` (or `n_pairs`) by `n_voxels` by `n_voxels`.  If a `summary_statistic` is supplied, or only two subjects are input, the singleton first dimension is removed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "zd2lahKfO40H"
   },
   "outputs": [],
   "source": [
    "# Compute ISFCs using leave-one-out approach\n",
    "isfcs, iscs = isfc(data, pairwise=False, vectorize_isfcs=True)\n",
    "\n",
    "# Check shape of output ISFC values\n",
    "print(f\"ISFC output shape = {isfcs.shape}\\ni.e., {isfcs.shape[0]} \"\n",
    "      f\"left-out subjects by {isfcs.shape[1]} connections (i.e., voxel pairs)\"\n",
    "      f\"\\nISCs output shape = {iscs.shape}\\ni.e., {iscs.shape[0]} \"\n",
    "      f\"left-out subjects by {iscs.shape[1]} voxels\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "MLNhWr54wlPP"
   },
   "source": [
    "Alternatively, we can retain the (redundant) structure of the ISFC matrices using `vectorize_isfcs=False` to yield a 3-dimensional array of shape `n_subjects` by `n_voxels` by `n_voxels`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "fbjL1mfbww6z"
   },
   "outputs": [],
   "source": [
    "# Compute ISFCs using leave-one-out approach\n",
    "isfcs = isfc(data, pairwise=False, vectorize_isfcs=False)\n",
    "\n",
    "# Check shape of output ISFC values\n",
    "print(f\"ISFC output shape = {isfcs.shape}\\ni.e., {isfcs.shape[0]} \"\n",
    "      f\"left-out subjects by {isfcs.shape[1]} voxels by {isfcs.shape[2]} \"\n",
    "      \"voxels\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "tViO_TIlxWEc"
   },
   "source": [
    "We can also supply a `summary_statistic` to collapse the ISFC values over left-out subjects or pairs of subjects:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "6IwDWoDnxjp9"
   },
   "outputs": [],
   "source": [
    "# Compute ISFCs using leave-one-out approach with mean\n",
    "isfcs, iscs = isfc(data, pairwise=False, summary_statistic='mean',\n",
    "                   vectorize_isfcs=True)\n",
    "\n",
    "# Check shape of output ISFC values\n",
    "print(f\"Mean ISFC output shape = {isfcs.shape}\\ni.e., {isfcs.shape[0]} \"\n",
    "      f\"connections (i.e., voxel pairs)\"\n",
    "      f\"\\nMean ISC output shape = {iscs.shape}\\ni.e., {iscs.shape[0]} \"\n",
    "      \"voxels\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Lt5AC1sF7pf_"
   },
   "source": [
    "We can use the `brainiak.isc.squareform_isfc` convenience function to convert between the condensed representation of ISFCs (with ISCs) and the square (redundant) representation of ISFCs. This function mimics `scipy.spatial.distance.squareform`, but retains the diagonal ISC values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "sdneE-BR7nT9"
   },
   "outputs": [],
   "source": [
    "from brainiak.isc import squareform_isfc\n",
    "\n",
    "# Start with square (redundant) ISFCs and check shape\n",
    "isfcs_sq = isfc(data, pairwise=False, vectorize_isfcs=False)\n",
    "print(f\"Square (redundant) ISFCs shape: {isfcs_sq.shape}\")\n",
    "\n",
    "# Convert these directly to condensed ISFCs (and ISCs)\n",
    "isfcs_c, iscs = squareform_isfc(isfcs_sq)\n",
    "print(f\"Condensed ISFCs shape: {isfcs_c.shape}, \"\n",
    "      f\"ISCs shape: {iscs.shape}\")\n",
    "\n",
    "# Convert these directly back to redundant ISFCs\n",
    "isfcs_r = squareform_isfc(isfcs_c, iscs)\n",
    "print(f\"Converted redundant ISFCs shape: {isfcs_r.shape}\")\n",
    "\n",
    "# Check that they are identical to the original square ISFCs\n",
    "assert np.array_equal(isfcs_sq, isfcs_r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "cmF52nSk-oTn"
   },
   "source": [
    "Let's confirm that the diagonal of the ISFC matrix represents each voxel correlated with itself across subjects—the conventional ISC described above. We can see that the conventional ISC analysis is in fact a subset of the ISFC analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "wthaT68i-5jy"
   },
   "outputs": [],
   "source": [
    "# Get ISC values directly from ISFC matrix\n",
    "isfcs, iscs = isfc(data, pairwise=False, vectorize_isfcs=True)\n",
    "\n",
    "# Check that these are the same as conventional ISCs\n",
    "assert np.allclose(iscs, isc(data))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "ZCsZs_xfx-t3"
   },
   "source": [
    "Finally, we can visualize the matrix of mean (or median) ISFC values. If we used `vectorize_isfcs=True`, we'll first need to apply `squareform_isfc` the ISFC (and ISC) values. The diagonal blocks represent the 10 artificial \"networks\" in our simulated data; the 100 voxels in each network are highly correlated with each other and largely uncorrelated with voxels in other networks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "WIbf8t-ichv3"
   },
   "outputs": [],
   "source": [
    "# Recompute mean ISFCs\n",
    "isfcs, iscs = isfc(data, pairwise=False, summary_statistic='mean',\n",
    "                   vectorize_isfcs=True)\n",
    "\n",
    "# Convert these to a square representation\n",
    "isfcs = squareform_isfc(isfcs, iscs)\n",
    "\n",
    "# Visual mean ISFC matrix\n",
    "plt.matshow(isfcs, cmap=\"RdYlBu_r\", vmin=-1, vmax=1)\n",
    "plt.grid(False)\n",
    "plt.xticks(np.arange(0, 1001, 100)[1:], np.arange(100, 1001, 100),\n",
    "           rotation=45)\n",
    "plt.gca().xaxis.tick_top()\n",
    "plt.gca().xaxis.set_label_position('top')\n",
    "plt.yticks(np.arange(0, 1001, 100)[1:], np.arange(100, 1001, 100))\n",
    "plt.xlabel('voxels')\n",
    "plt.ylabel('voxels')\n",
    "ax = plt.gca()\n",
    "plt.colorbar(fraction=0.046, pad=0.04);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "PvGplkL506iG"
   },
   "source": [
    "We can inject some structure into our simulated data to yield a more realistic ISFC matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "OHj-L6w718qK"
   },
   "outputs": [],
   "source": [
    "# Create more structured simulated data with 7 \"networks\";\n",
    "# don't worry about the details\n",
    "from scipy.ndimage import gaussian_filter1d\n",
    "\n",
    "def structured_timeseries(n_subjects, n_TRs, n_voxels=1000, noise=1):\n",
    "  signals = np.random.randn(n_TRs, 3)\n",
    "  networks = np.column_stack((signals + np.random.randn(n_TRs, 3) * noise,\n",
    "                              signals[:, 0] + np.random.randn(n_TRs) * noise,\n",
    "                              signals[:, 0] + np.random.randn(n_TRs) * noise,\n",
    "                              -signals[:, 2] + np.random.randn(n_TRs) * noise,\n",
    "                              signals[:, 2] + np.random.randn(n_TRs) * noise))\n",
    "  networks = networks[:, [0, 3, 4, 5, 1, 2, 6]]\n",
    "  six = np.random.randint(n_voxels // 20, n_voxels // 6, 6)\n",
    "  seven = np.append(six, (n_voxels - np.sum(six)))\n",
    "  voxels = np.column_stack([np.tile(network[:, np.newaxis], (1, extent))\n",
    "                            for network, extent in zip(networks.T, seven)])\n",
    "  areas = [0] + sorted(np.random.randint(0, 1000, 16))\n",
    "  areas = np.diff(areas).tolist() + [(1000 - areas[-1])]\n",
    "  noise_sources = np.random.randn(n_TRs, 7)\n",
    "  structured_noise = np.column_stack([np.tile(\n",
    "      (noise_sources[:, np.random.choice(range(7))] *\n",
    "       np.random.choice([-1, 1, 1, 1]))[:, np.newaxis],\n",
    "                                              (1, extent)) \n",
    "                                      for extent in areas])\n",
    "  voxels = gaussian_filter1d(voxels, 8.0, axis=0)\n",
    "  structured_noise = gaussian_filter1d(structured_noise, 8.0, axis=0)\n",
    "  data = []\n",
    "  for s in np.arange(n_subjects):\n",
    "    data.append(voxels + structured_noise * noise * .2 +\n",
    "                np.random.randn(n_TRs, n_voxels) * noise * 1.35)\n",
    "\n",
    "  data = np.dstack(data)\n",
    "  return data\n",
    "\n",
    "structured_data = structured_timeseries(n_subjects, n_TRs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "b3VBL6M4EyTx"
   },
   "source": [
    "Now, we can recompute mean ISFCs using the leave-one-out approach and visualize the resulting ISFC matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "e_mY31oc5MbX"
   },
   "outputs": [],
   "source": [
    "# Compute ISFCs using leave-one-out approach with mean\n",
    "isfcs, iscs = isfc(structured_data, pairwise=False, summary_statistic='mean',\n",
    "                   vectorize_isfcs=True)\n",
    "\n",
    "# Convert these to a square representation\n",
    "isfcs = squareform_isfc(isfcs, iscs)\n",
    "\n",
    "# Visual mean ISFC matrix\n",
    "plt.matshow(isfcs, cmap=\"RdYlBu_r\", vmin=-.3, vmax=.3)\n",
    "plt.grid(False)\n",
    "plt.xticks(np.arange(0, 1001, 100)[1:], np.arange(100, 1001, 100),\n",
    "           rotation=45)\n",
    "plt.gca().xaxis.tick_top()\n",
    "plt.gca().xaxis.set_label_position('top')\n",
    "plt.yticks(np.arange(0, 1001, 100)[1:], np.arange(100, 1001, 100))\n",
    "plt.xlabel('voxels')\n",
    "plt.ylabel('voxels')\n",
    "ax = plt.gca()\n",
    "plt.colorbar(fraction=0.046, pad=0.04);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "llM5hK-Fd51Y"
   },
   "source": [
    "## Real fMRI data\n",
    "Next, we'll  download a publicly available fMRI dataset and run an ISC analysis. This dataset comprises fMRI data for 20 subjects listening to the spoken story [Pie Man](https://themoth.org/stories/pie-man) by Jim O'Grady (archived on the [Princeton DataSpace](https://dataspace.princeton.edu/jspui/handle/88435/dsp01dz010s83s)). Note that we use 20 subjects to minimize computational demands for this tutorial and recommend larger sample sizes for publication. The gzipped data archive file is ~1.5 GB in size, and may take a couple minutes to download and unzip. The functional data were acquired with 3 x 3 x 4 mm voxels and 1.5 s TRs. Data were preprocessed using [fMRIPrep](https://fmriprep.readthedocs.io/en/stable/) ([Esteban et al., 2018](https://doi.org/10.1038/s41592-018-0235-4)), including spatial normalization to MNI space (the T1-weighted [ICBM 2009c Nonlinear Asymmetric template](http://nist.mni.mcgill.ca/?p=904)). The data were then smoothed to 6 mm FWHM using [AFNI](https://afni.nimh.nih.gov/)'s [3dBlurToFWHM](https://afni.nimh.nih.gov/pub/dist/doc/program_help/3dBlurToFWHM.html) ([Cox, 1996](https://doi.org/10.1006/cbmr.1996.0014)). The following confound variables were regressed out using [3dTproject](https://afni.nimh.nih.gov/pub/dist/doc/program_help/3dTproject.html): six head motion parameters (and their first derivatives), framewise displacement, six prinicipal components from an anatomical mask of cerebrospinal fluid (CSF) and white matter, sine/cosine bases for high-pass filtering (cutoff: 0.00714 Hz; 140 s), as well as a linear and quadratic trends. The anatomical template and a brain mask (i.e., excluding skull) are supplied as well. These have been resampled to match resolution of the functional images."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "ozcGwfP4vtRN"
   },
   "outputs": [],
   "source": [
    "# Download data tarball from Princeton DataSpace\n",
    "from urllib.request import urlretrieve\n",
    "urlretrieve('https://dataspace.princeton.edu/jspui/bitstream/'\n",
    "            '88435/dsp01dz010s83s/6/pieman-isc-tutorial.tgz',\n",
    "            'pieman-isc-tutorial.tgz');\n",
    "!tar -xvzf pieman-isc-tutorial.tgz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "_XOoFmG_Fv2F"
   },
   "source": [
    "### Loading MRI data\n",
    "We'll use NiBabel as well as BrainIAK's `io` and `image` functionality to load the functional data and apply a brain mask."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "SmnvH1VUYkv7"
   },
   "outputs": [],
   "source": [
    "# Import functions helpful for managing file paths\n",
    "from glob import glob\n",
    "from os.path import join\n",
    "\n",
    "data_dir = 'pieman-isc-tutorial'\n",
    "\n",
    "# Filenames for MRI data; gzipped NIfTI images (.nii.gz)\n",
    "func_fns = glob(join(data_dir, ('sub-*_task-pieman_space-MNI152NLin2009cAsym'\n",
    "                                '_desc-tproject_bold.nii.gz')))\n",
    "mask_fn = join(data_dir, 'MNI152NLin2009cAsym_desc-brain_mask.nii.gz')\n",
    "mni_fn = join(data_dir, 'MNI152NLin2009cAsym_desc-brain_T1w.nii.gz')\n",
    "\n",
    "# Load a NIfTI of the brain mask as a reference Nifti1Image\n",
    "ref_nii = nib.load(mask_fn)\n",
    "\n",
    "# Load functional images and masks using brainiak.io\n",
    "func_imgs = load_images(func_fns)\n",
    "mask_img = load_boolean_mask(mask_fn)\n",
    "\n",
    "# Get coordinates of mask voxels in original image\n",
    "mask_coords = np.where(mask_img)\n",
    "\n",
    "# Apply the brain mask using brainiak.image\n",
    "masked_imgs = mask_images(func_imgs, mask_img)\n",
    "\n",
    "# Collate data into a single TR x voxel x subject array\n",
    "orig_data = MaskedMultiSubjectData.from_masked_images(masked_imgs,\n",
    "                                                      len(func_fns))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "t5LQcCZyLKA0"
   },
   "source": [
    "The data from each subject are stacked along the third dimension, yielding a `n_TRs` by `n_voxels` by `n_subjects` array. The functional acquisition originally included 13 s of music and 2 s of silence prepended to the story stimulus and an additional 13 s of silence after the story (450 s or 300 TRs in total). These segments as well as the first 12 s (8 TRs) after story onset can be discarded to minimize stimulus onset/offset effects. We may also opt to z-score the time series for each voxel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "SyMpEL5t1wM5"
   },
   "outputs": [],
   "source": [
    "print(f\"Original fMRI data shape: {orig_data.shape} \"\n",
    "      f\"\\ni.e., {orig_data.shape[0]} time points, {orig_data.shape[1]} voxels, \"\n",
    "      f\"{orig_data.shape[2]} subjects\")\n",
    "\n",
    "# Trim off non-story TRs and 12 s post-onset\n",
    "data = orig_data[18:-8, ...]\n",
    "\n",
    "print(f\"Trimmed fMRI data shape: {data.shape} \"\n",
    "      f\"\\ni.e., {data.shape[0]} time points, {data.shape[1]} voxels, \"\n",
    "      f\"{data.shape[2]} subjects\")\n",
    "\n",
    "# Z-score time series for each voxel\n",
    "data = zscore(data, axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "5dQSyCOxceRn"
   },
   "source": [
    "### ISC analysis\n",
    "Next, we'll run a leave-one-out ISC analysis on the preprocessed fMRI data, including all voxels in the brain mask—this may take a few minutes. Note that some voxels with no variance over time for one or more subjects were included in the brain mask due to the limited field of view during EPI acquisition and susceptibility artifacts (signal dropout). This yields NaN (not a number) values. Here, we'll use set `tolerate_nans` to `0.8` to ensure that, when computing the averaging time series for *N*–1 subjects in the leave-one-out approach, only voxels with >= 80% of subjects have non-NaN values are included."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "8uwaP0ph2Ie2"
   },
   "outputs": [],
   "source": [
    "# Leave-one-out approach\n",
    "iscs = isc(data, pairwise=False, tolerate_nans=.8)\n",
    "\n",
    "\n",
    "# Check shape of output ISC values\n",
    "print(f\"ISC values shape = {iscs.shape} \\ni.e., {iscs.shape[0]} \"\n",
    "      f\"left-out subjects and {iscs.shape[1]} voxel(s)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "5Dn56Iog8rE5"
   },
   "source": [
    "Since we didn't supply a `summary_statistic` in the `isc` call, we get an ISC value for each left-out subject (we'll need ISCs for each subject in the subsequent statistical test). If we want to preliminarily inspect the mean (or median) ISCs, we can apply the `brainiak.isc.compute_summary_statistic` function afterward. Note that if we specifed a `summary_statistic` in the `isc` call, the `isc` function would simply use `compute_summary_statistic` internally."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "4qkbnT-18qHh"
   },
   "outputs": [],
   "source": [
    "# Compute mean ISC (with Fisher transformation)\n",
    "mean_iscs = compute_summary_statistic(iscs, summary_statistic='mean', axis=0)\n",
    "\n",
    "print(f\"ISC values shape = {mean_iscs.shape} \\ni.e., {mean_iscs.shape[0]} \"\n",
    "      f\"mean value across left-out subjects and {iscs.shape[1]} voxel(s)\"\n",
    "      f\"\\nMinimum mean ISC across voxels = {np.nanmin(mean_iscs):.3f}; \"\n",
    "      f\"maximum mean ISC across voxels = {np.nanmax(mean_iscs):.3f}\")\n",
    "\n",
    "\n",
    "# Compute median ISC\n",
    "median_iscs = compute_summary_statistic(iscs, summary_statistic='median',\n",
    "                                        axis=0)\n",
    "\n",
    "print(f\"ISC values shape = {median_iscs.shape} \\ni.e., {median_iscs.shape[0]} \"\n",
    "      f\"median value across left-out subjects and {iscs.shape[1]} voxel(s)\"\n",
    "      f\"\\nMinimum median ISC across voxels = {np.nanmin(median_iscs):.3f}; \"\n",
    "      f\"maximum median ISC across voxels = {np.nanmax(median_iscs):.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "kcFv1DMkXsiV"
   },
   "source": [
    "### Statistical testing\n",
    "To test whether the observed ISCs are significantly greater than zero, we'll perform a bootstrap hypothesis test ([Chen et al., 2016](https://doi.org/10.1016/j.neuroimage.2016.05.023)). This may also take a couple minutes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "WzIRdU6OeLjv"
   },
   "outputs": [],
   "source": [
    "# Run bootstrap hypothesis test on ISCs\n",
    "observed, ci, p, distribution = bootstrap_isc(iscs, pairwise=False,\n",
    "                                              ci_percentile=95,\n",
    "                                              summary_statistic='median',\n",
    "                                              n_bootstraps=1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "xPq-KbyO8BBb"
   },
   "source": [
    "Before we correct for multiple tests, we should exclude any voxels with NaNs. To do this, we'll extract the non-NaN voxels, run the correction for multiple tests, then reinsert the non-NaN voxels into the full mask."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "aO21dj9rOriD"
   },
   "outputs": [],
   "source": [
    "# Get number of NaN voxels\n",
    "n_nans = np.sum(np.isnan(observed))\n",
    "print(f\"{n_nans} voxels out of {observed.shape[0]} are NaNs \"\n",
    "      f\"({n_nans / observed.shape[0] * 100:.2f}%)\")\n",
    "\n",
    "# Get voxels without NaNs\n",
    "nonnan_mask = ~np.isnan(observed)\n",
    "nonnan_coords = np.where(nonnan_mask)\n",
    "\n",
    "# Mask both the ISC and p-value map to exclude NaNs\n",
    "nonnan_isc = observed[nonnan_mask]\n",
    "nonnan_p = p[nonnan_mask]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "qKIFZxx9Aea7"
   },
   "source": [
    "Now we'll apply the `multipletests` function from StatsModels to the *p*-values from the bootstrap hypothesis test to control the false discovery rate (FDR) at 0.05 across all voxels. This yields a map of *q*-values. We can then threshold our ISC image based on the FDR-adjusted *q*-values, which are derived from the entire image, rather than the uncorrected *p*-values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "Jk0mZui0elF9"
   },
   "outputs": [],
   "source": [
    "# Get FDR-controlled q-values\n",
    "nonnan_q = multipletests(nonnan_p, method='fdr_by')[1]\n",
    "threshold = .05\n",
    "print(f\"{np.sum(nonnan_q < threshold)} significant voxels \"\n",
    "      f\"controlling FDR at {threshold}\")\n",
    "\n",
    "# Threshold ISCs according FDR-controlled threshold\n",
    "nonnan_isc[nonnan_q >= threshold] = np.nan\n",
    "\n",
    "# Reinsert thresholded ISCs back into whole brain image\n",
    "isc_thresh = np.full(observed.shape, np.nan)\n",
    "isc_thresh[nonnan_coords] = nonnan_isc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "fydiRkFyZLCw"
   },
   "source": [
    "### Visualizing results\n",
    "Finally, to visualize the significant ISC values, we must first reformat the 2-dimensional masked array into a 3-dimensional NIfTI image. We'll use an arbitrary reference NIfTI image `ref_nii` (the brain mask) to assign affine and header information correctly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "0NTEM7jCgRCO"
   },
   "outputs": [],
   "source": [
    "# Create empty 3D image and populate\n",
    "# with thresholded ISC values\n",
    "isc_img = np.full(ref_nii.shape, np.nan)\n",
    "isc_img[mask_coords] = isc_thresh\n",
    "\n",
    "# Convert to NIfTI image\n",
    "isc_nii = nib.Nifti1Image(isc_img, ref_nii.affine, ref_nii.header)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "2Ew86aC-ac13"
   },
   "source": [
    "We'll use `nilearn.plotting.plot_stat_map` to plot two views of the NIfTI image. We'll set the maximum ISC value for the colorbar at 0.5 and use a divergent colormap called `RdYlBu_r`. This yields maps where significant voxels are colored according to the median ISC value across left-out subjects. Statistical significance was assessed by a \n",
    "nonparametric bootstrap hypothesis test resampling left-out subjects and corrected for multiple tests by controlling FDR at .05.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "J0duUdiOacGI"
   },
   "outputs": [],
   "source": [
    "# Install Nilearn for plotting functionality\n",
    "!pip install nilearn\n",
    "from nilearn.plotting import plot_stat_map\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "R7glFCGhwDvW"
   },
   "outputs": [],
   "source": [
    "# Plot slices at coordinates -61, -20, 8\n",
    "plot_stat_map(\n",
    "    isc_nii,\n",
    "    cmap='RdYlBu_r',\n",
    "    vmax=.5,\n",
    "    cut_coords=(-61, -20, 8))\n",
    "\n",
    "# Plot slices at coordinates 0, -65, 40\n",
    "plot_stat_map(\n",
    "    isc_nii,\n",
    "    cmap='RdYlBu_r',\n",
    "    vmax=.5,\n",
    "    cut_coords=(0, -65, 40))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "5w1mXQa1ez3t"
   },
   "source": [
    "Some significant voxels have fairly low ISCs, so we can also use the `threshold` option in `plot_stat_map` to also exclude voxels with, e.g., ISC < .1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "HodP2a95ey_K"
   },
   "outputs": [],
   "source": [
    "# Plot slices at coordinates -61, -20, 8\n",
    "plot_stat_map(\n",
    "    isc_nii,\n",
    "    cmap='RdYlBu_r',\n",
    "    vmax=.5,\n",
    "    threshold=.1,\n",
    "    cut_coords=(-61, -20, 8))\n",
    "\n",
    "# Plot slices at coordinates 0, -65, 40\n",
    "plot_stat_map(\n",
    "    isc_nii,\n",
    "    cmap='RdYlBu_r',\n",
    "    vmax=.5,\n",
    "    threshold=.1,\n",
    "    cut_coords=(0, -65, 40))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "qCR8Qn_FfevP"
   },
   "source": [
    "Note that here we are only analyzing responses to the fully intact Pie Man story, which yields high ISCs in both low-level auditory areas and higher-level areas processing narrative content. However, if we scrambled the presentation of the story so as to disrupt the temporally contiguous narrative, we would expect to see high ISCs only in low-level auditory areas ([Hasson et al., 2008](https://doi.org/10.1523/jneurosci.5487-07.2008); [Lerner et al., 2011](https://doi.org/10.1523/jneurosci.3684-10.2011)). To finish, we can also use `nib.save` to save the NIfTI image for use with other neuroimaging data analysis and visualization programs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "g2fSRy6RF1iB"
   },
   "outputs": [],
   "source": [
    "# Save final ISC NIfTI image as .nii\n",
    "isc_fn = 'isc_thresh_pieman_n20.nii.gz'\n",
    "nib.save(isc_nii, isc_fn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "gLLGJ1HCFH3E"
   },
   "source": [
    "## References and suggested reading\n",
    "\n",
    "* Benjamini, Y., & Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society: Series B (Methodological), 289–300. https://www.jstor.org/stable/2346101\n",
    "\n",
    "* Benjamini, Y., & Yekutieli, D. (2001). The control of the false discovery rate in multiple testing under dependency. *Annals of Statistics*, *29*(4), 1165–1188. https://www.jstor.org/stable/2674075\n",
    "\n",
    "* Chen, G., Shin, Y. W., Taylor, P. A., Glen, D. R., Reynolds, R. C., Israel, R. B., & Cox, R. W. (2016). Untangling the relatedness among correlations, part I: nonparametric approaches to inter-subject correlation analysis at the group level. *NeuroImage*, *142*, 248–259. https://doi.org/10.1016/j.neuroimage.2016.05.023\n",
    "\n",
    "* Chen, G., Taylor, P. A., Shin, Y. W., Reynolds, R. C., & Cox, R. W. (2017). Untangling the relatedness among correlations, part II: inter-subject correlation group analysis through linear mixed-effects modeling. *NeuroImage*, *147*, 825–840. https://doi.org/10.1016/j.neuroimage.2016.08.029\n",
    "\n",
    "* Cox, R. W. (1996). AFNI: software for analysis and visualization of functional magnetic resonance neuroimages. *Computers and Biomedical Research*, *29*(3), 162–173. https://doi.org/10.1006/cbmr.1996.0014\n",
    "\n",
    "* Esteban, O., Markiewicz, C., Blair, R. W., Moodie, C., Isik, A. I., Erramuzpe, A., Kent, J. D., Goncalves, M., DuPre, E., Snyder, M., Oya, H., Ghosh, S., Wright, J., Durnez, J., Poldrack, R., & Gorgolewski, K. J. (2018). fMRIPrep: a robust preprocessing pipeline for functional MRI. *Nature Methods*. https://doi.org/10.1038/s41592-018-0235-4\n",
    "\n",
    "* Genovese, C. R., Lazar, N. A., & Nichols, T. (2002). Thresholding of statistical maps in functional neuroimaging using the false discovery rate. *NeuroImage*, *15*(4), 870–878. https://doi.org/10.1006/nimg.2001.1037\n",
    "\n",
    "* Hasson, U., Ghazanfar, A. A., Galantucci, B., Garrod, S., & Keysers, C. (2012). Brain-to-brain coupling: a mechanism for creating and sharing a social world. *Trends in Cognitive Sciences*, *16*(2), 114–121. https://doi.org/10.1016/j.tics.2011.12.007\n",
    "\n",
    "* Hasson, U., Malach, R., & Heeger, D. J. (2010). Reliability of cortical activity during natural stimulation. *Trends in Cognitive Sciences*, *14*(1), 40–48. https://doi.org/10.1016/j.tics.2009.10.011\n",
    "\n",
    "* Hasson, U., Nir, Y., Levy, I., Fuhrmann, G., & Malach, R. (2004). Intersubject synchronization of cortical activity during natural vision. *Science*, *303*(5664), 1634–1640. https://doi.org/10.1126/science.1089506\n",
    "\n",
    "* Hasson, U., Yang, E., Vallines, I., Heeger, D. J., & Rubin, N. (2008). A hierarchy of temporal receptive windows in human cortex. *Journal of Neuroscience*, *28*(10), 2539–2550. https://doi.org/10.1523/jneurosci.5487-07.2008\n",
    "\n",
    "* Kauppi, J. P., Pajula, J., & Tohka, J. (2014). A versatile software package for inter-subject correlation based analyses of fMRI. *Frontiers in Neuroinformatics*, *8*, 2. https://doi.org/10.3389/fninf.2014.00002\n",
    "\n",
    "* Lerner, Y., Honey, C. J., Silbert, L. J., & Hasson, U. (2011). Topographic mapping of a hierarchy of temporal receptive windows using a narrated story. *Journal of Neuroscience*, *31*(8), 2906–2915. https://doi.org/10.1523/jneurosci.3684-10.2011\n",
    "\n",
    "* Nichols, T. E., & Holmes, A. P. (2002). Nonparametric permutation tests for functional neuroimaging: a primer with examples. *Human Brain Mapping*, *15*(1), 1–25. https://doi.org/10.1002/hbm.1058\n",
    "\n",
    "* Silbert, L. J., Honey, C. J., Simony, E., Poeppel, D., & Hasson, U. (2014). Coupled neural systems underlie the production and comprehension of naturalistic narrative speech. *Proceedings of the National Academy of Sciences of the United States of America*, *111*(43), E4687–E4696. https://doi.org/10.1073/pnas.1323812111\n",
    "\n",
    "* Simony, E., Honey, C. J., Chen, J., Lositsky, O., Yeshurun, Y., Wiesel, A., & Hasson, U. (2016). Dynamic reconfiguration of the default mode network during narrative comprehension. *Nature Communications*, *7*, 12141. https://doi.org/10.1038/ncomms12141\n",
    "\n",
    "* Stephens, G. J., Silbert, L. J., & Hasson, U. (2010). Speaker–listener neural coupling underlies successful communication. *Proceedings of the National Academy of Sciences of the United States of America*, *107*(32), 14425–14430. https://doi.org/10.1073/pnas.1008662107"
   ]
  }
 ],
 "metadata": {
  "accelerator": "GPU",
  "colab": {
   "collapsed_sections": [],
   "include_colab_link": true,
   "name": "isc_tutorial.ipynb",
   "provenance": [],
   "toc_visible": true,
   "version": "0.3.2"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}