5-group_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Group analysis with MNE\n",
    "\n",
    "\n",
    "The aim of this lecture is to show you how to do group analysis\n",
    "with MNE-Python.\n",
    "\n",
    "    Authors: Alexandre Gramfort <alexandre.gramfort@inria.fr>\n",
    "             Denis A. Engemann <denis.engemann@gmail.com>\n",
    "             Mainak Jas <mainakjas@gmail.com>\n",
    "             Hicham Janati <hicham.janati@inria.fr>\n",
    "\n",
    "    License: BSD (3-clause)\n",
    "   \n",
    "The code in this notebook is taken from https://github.com/mne-tools/mne-biomag-group-demo/tree/master/scripts/results/statistics\n",
    "\n",
    "You'll need to [download the data from here](https://www.dropbox.com/s/eh8yexs8d5fx3kp/ds000117-mne-group.zip?dl=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib qt\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import os\n",
    "import numpy as np\n",
    "import mne\n",
    "\n",
    "mne.set_log_level('WARNING')\n",
    "\n",
    "# Change the following path to where the folder ds000117-practical is on your disk\n",
    "data_path = os.path.expanduser(\"~/work/data/ds000117-mne-group/\")\n",
    "evokeds_path = os.path.join(data_path, \"evokeds\")\n",
    "stc_path = os.path.join(data_path, \"stcs\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Subjects\n",
    "\n",
    "we have here the evoked data and the source time courses (STCs) for\n",
    "all but 3 subjects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['sub002',\n",
       " 'sub003',\n",
       " 'sub004',\n",
       " 'sub006',\n",
       " 'sub007',\n",
       " 'sub008',\n",
       " 'sub009',\n",
       " 'sub010',\n",
       " 'sub011',\n",
       " 'sub012',\n",
       " 'sub013',\n",
       " 'sub014',\n",
       " 'sub015',\n",
       " 'sub017',\n",
       " 'sub018',\n",
       " 'sub019']"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "subjects = [\"sub%03d\" % i for i in range(1, 20) if i not in [1, 5, 16]]\n",
    "subjects"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Let's look at our different subjects"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<Evoked  |  'famous' (average, N=104), [-0.2, 2.9] sec, 376 ch, ~9.3 MB>,\n",
       " <Evoked  |  'scrambled' (average, N=84), [-0.2, 2.9] sec, 376 ch, ~9.3 MB>,\n",
       " <Evoked  |  'unfamiliar' (average, N=101), [-0.2, 2.9] sec, 376 ch, ~9.3 MB>,\n",
       " <Evoked  |  'contrast' (average, N=60), [-0.2, 2.9] sec, 376 ch, ~9.3 MB>,\n",
       " <Evoked  |  'faces' (average, N=205), [-0.2, 2.9] sec, 376 ch, ~9.3 MB>,\n",
       " <Evoked  |  'faces_eq' (average, N=84), [-0.2, 2.9] sec, 376 ch, ~9.3 MB>,\n",
       " <Evoked  |  'scrambled_eq' (average, N=84), [-0.2, 2.9] sec, 376 ch, ~9.3 MB>]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "subject = 'sub015'\n",
    "fname = os.path.join(evokeds_path, \"%s-highpass-ave.fif\" % subject)\n",
    "evokeds = mne.read_evokeds(fname, verbose=False)\n",
    "evokeds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "evokeds[0].plot(spatial_colors=True);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now loop over all subjects to see if there are any subject\n",
    "that look problematic."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "conditions = [\"famous\", \"scrambled\", \"faces\"]\n",
    "\n",
    "f, axes = plt.subplots(4, 4, figsize=(13, 9), sharex=True, sharey=True)\n",
    "\n",
    "for ax, subject in zip(axes.ravel(), subjects):\n",
    "    evokeds_dict = dict()\n",
    "    fname = os.path.join(evokeds_path, \"%s-highpass-ave.fif\" % subject)\n",
    "    evokeds = mne.read_evokeds(fname)\n",
    "    evokeds = [ev for ev in evokeds if ev.comment in conditions]\n",
    "    for condition, ev in zip(conditions, evokeds):\n",
    "        evokeds_dict[condition] = ev.crop(tmax=0.5)\n",
    "    mne.viz.plot_compare_evokeds(evokeds_dict, picks=\"grad\", show=False,\n",
    "                                 axes=ax, title=subject)\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    <b>EXERCISE</b>:\n",
    "     <ul>\n",
    "      <li>Compute the same type of GFP plots for EEG and MEG Magnetometers</li>\n",
    "      <li>Do you see any problematic subject?</li>\n",
    "    </ul>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    <b>EXERCISE</b>:\n",
    "     <ul>\n",
    "      <li>Read all the evoked data from all subjects and store the ones for conditions 'faces' or 'scrambled' in a Python list.</li>\n",
    "      <li>Do a GFP plots of the contrast between \"faces\" vs \"scrambled conditions. You will store the contrast of all subjects in a list.</li>\n",
    "    </ul>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[<Evoked  |  'scrambled' (average, N=238), [-0.2, 2.9] sec, 376 ch, ~9.3 MB>,\n",
       "  <Evoked  |  'faces' (average, N=465), [-0.2, 2.9] sec, 376 ch, ~9.3 MB>],\n",
       " [<Evoked  |  'scrambled' (average, N=211), [-0.2, 2.9] sec, 376 ch, ~9.3 MB>,\n",
       "  <Evoked  |  'faces' (average, N=470), [-0.2, 2.9] sec, 376 ch, ~9.3 MB>],\n",
       " [<Evoked  |  'scrambled' (average, N=216), [-0.2, 2.9] sec, 376 ch, ~9.3 MB>,\n",
       "  <Evoked  |  'faces' (average, N=436), [-0.2, 2.9] sec, 376 ch, ~9.3 MB>],\n",
       " [<Evoked  |  'scrambled' (average, N=181), [-0.2, 2.9] sec, 376 ch, ~9.3 MB>,\n",
       "  <Evoked  |  'faces' (average, N=357), [-0.2, 2.9] sec, 376 ch, ~9.3 MB>],\n",
       " [<Evoked  |  'scrambled' (average, N=286), [-0.2, 2.9] sec, 376 ch, ~9.3 MB>,\n",
       "  <Evoked  |  'faces' (average, N=581), [-0.2, 2.9] sec, 376 ch, ~9.3 MB>],\n",
       " [<Evoked  |  'scrambled' (average, N=214), [-0.2, 2.9] sec, 376 ch, ~9.3 MB>,\n",
       "  <Evoked  |  'faces' (average, N=444), [-0.2, 2.9] sec, 376 ch, ~9.3 MB>],\n",
       " [<Evoked  |  'scrambled' (average, N=128), [-0.2, 2.9] sec, 376 ch, ~9.3 MB>,\n",
       "  <Evoked  |  'faces' (average, N=296), [-0.2, 2.9] sec, 376 ch, ~9.3 MB>],\n",
       " [<Evoked  |  'scrambled' (average, N=177), [-0.2, 2.9] sec, 376 ch, ~9.3 MB>,\n",
       "  <Evoked  |  'faces' (average, N=353), [-0.2, 2.9] sec, 376 ch, ~9.3 MB>],\n",
       " [<Evoked  |  'scrambled' (average, N=176), [-0.2, 2.9] sec, 376 ch, ~9.3 MB>,\n",
       "  <Evoked  |  'faces' (average, N=352), [-0.2, 2.9] sec, 376 ch, ~9.3 MB>],\n",
       " [<Evoked  |  'scrambled' (average, N=125), [-0.2, 2.9] sec, 376 ch, ~9.3 MB>,\n",
       "  <Evoked  |  'faces' (average, N=269), [-0.2, 2.9] sec, 376 ch, ~9.3 MB>],\n",
       " [<Evoked  |  'scrambled' (average, N=196), [-0.2, 2.9] sec, 376 ch, ~9.3 MB>,\n",
       "  <Evoked  |  'faces' (average, N=363), [-0.2, 2.9] sec, 376 ch, ~9.3 MB>],\n",
       " [<Evoked  |  'scrambled' (average, N=174), [-0.2, 2.9] sec, 376 ch, ~9.3 MB>,\n",
       "  <Evoked  |  'faces' (average, N=327), [-0.2, 2.9] sec, 376 ch, ~9.3 MB>],\n",
       " [<Evoked  |  'scrambled' (average, N=84), [-0.2, 2.9] sec, 376 ch, ~9.3 MB>,\n",
       "  <Evoked  |  'faces' (average, N=205), [-0.2, 2.9] sec, 376 ch, ~9.3 MB>],\n",
       " [<Evoked  |  'scrambled' (average, N=127), [-0.2, 2.9] sec, 376 ch, ~9.3 MB>,\n",
       "  <Evoked  |  'faces' (average, N=299), [-0.2, 2.9] sec, 376 ch, ~9.3 MB>],\n",
       " [<Evoked  |  'scrambled' (average, N=182), [-0.2, 2.9] sec, 376 ch, ~9.3 MB>,\n",
       "  <Evoked  |  'faces' (average, N=401), [-0.2, 2.9] sec, 376 ch, ~9.3 MB>],\n",
       " [<Evoked  |  'scrambled' (average, N=243), [-0.2, 2.9] sec, 376 ch, ~9.3 MB>,\n",
       "  <Evoked  |  'faces' (average, N=457), [-0.2, 2.9] sec, 376 ch, ~9.3 MB>]]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# TODO\n",
    "evokeds_list = []\n",
    "\n",
    "for subject in subjects:\n",
    "    fname = os.path.join(evokeds_path, \"%s-highpass-ave.fif\" % subject)\n",
    "    evokeds = mne.read_evokeds(fname)\n",
    "    evokeds = [ev for ev in evokeds if ev.comment in ['faces', 'scrambled']]\n",
    "    evokeds_list.append(evokeds)\n",
    "\n",
    "evokeds_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO\n",
    "\n",
    "contrast_list = []\n",
    "f, axes = plt.subplots(4, 4, figsize=(13, 9), sharex=True, sharey=True)\n",
    "\n",
    "for ax, subject, evokeds in zip(axes.ravel(), subjects, evokeds_list):\n",
    "    contrast = mne.combine_evoked(evokeds, weights=[0.5, -0.5])\n",
    "    contrast.crop(None, 0.4)\n",
    "    contrast_list.append(contrast)\n",
    "    mne.viz.plot_compare_evokeds(contrast, picks=\"grad\", show=False,\n",
    "                                 axes=ax, title=subject)\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Look at grand averages\n",
    "\n",
    "Grand averages are obtained by averaging the sensor space data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<Evoked  |  '0.500 * scrambled + -0.500 * faces' (average, N=629.701280227596), [-0.2, 0.4] sec, 376 ch, ~7.7 MB>,\n",
       " <Evoked  |  '0.500 * scrambled + -0.500 * faces' (average, N=582.496328928047), [-0.2, 0.4] sec, 376 ch, ~7.7 MB>,\n",
       " <Evoked  |  '0.500 * scrambled + -0.500 * faces' (average, N=577.7668711656441), [-0.2, 0.4] sec, 376 ch, ~7.7 MB>,\n",
       " <Evoked  |  '0.500 * scrambled + -0.500 * faces' (average, N=480.42379182156134), [-0.2, 0.4] sec, 376 ch, ~7.7 MB>,\n",
       " <Evoked  |  '0.500 * scrambled + -0.500 * faces' (average, N=766.6251441753171), [-0.2, 0.4] sec, 376 ch, ~7.7 MB>,\n",
       " <Evoked  |  '0.500 * scrambled + -0.500 * faces' (average, N=577.6048632218844), [-0.2, 0.4] sec, 376 ch, ~7.7 MB>,\n",
       " <Evoked  |  '0.500 * scrambled + -0.500 * faces' (average, N=357.4339622641509), [-0.2, 0.4] sec, 376 ch, ~7.7 MB>,\n",
       " <Evoked  |  '0.500 * scrambled + -0.500 * faces' (average, N=471.5547169811321), [-0.2, 0.4] sec, 376 ch, ~7.7 MB>,\n",
       " <Evoked  |  '0.500 * scrambled + -0.500 * faces' (average, N=469.33333333333337), [-0.2, 0.4] sec, 376 ch, ~7.7 MB>,\n",
       " <Evoked  |  '0.500 * scrambled + -0.500 * faces' (average, N=341.3705583756345), [-0.2, 0.4] sec, 376 ch, ~7.7 MB>,\n",
       " <Evoked  |  '0.500 * scrambled + -0.500 * faces' (average, N=509.1091234347049), [-0.2, 0.4] sec, 376 ch, ~7.7 MB>,\n",
       " <Evoked  |  '0.500 * scrambled + -0.500 * faces' (average, N=454.2754491017964), [-0.2, 0.4] sec, 376 ch, ~7.7 MB>,\n",
       " <Evoked  |  '0.500 * scrambled + -0.500 * faces' (average, N=238.33910034602076), [-0.2, 0.4] sec, 376 ch, ~7.7 MB>,\n",
       " <Evoked  |  '0.500 * scrambled + -0.500 * faces' (average, N=356.55399061032864), [-0.2, 0.4] sec, 376 ch, ~7.7 MB>,\n",
       " <Evoked  |  '0.500 * scrambled + -0.500 * faces' (average, N=500.73413379073753), [-0.2, 0.4] sec, 376 ch, ~7.7 MB>,\n",
       " <Evoked  |  '0.500 * scrambled + -0.500 * faces' (average, N=634.5771428571428), [-0.2, 0.4] sec, 376 ch, ~7.7 MB>]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "contrast_list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    <b>EXERCISE</b>:\n",
    "     <ul>\n",
    "      <li>Using the mne.grand_average function see the average constract in your group.</li>\n",
    "      <li>Comment what you see.</li>\n",
    "    </ul>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "mne.grand_average?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "evoked_gave = mne.grand_average(contrast_list)\n",
    "evoked_gave.plot(spatial_colors=True);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Let's so some statistics\n",
    "\n",
    "We'll start with a single channel EEG065"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "stat_fun(H1): min=-1.568384 max=8.532331\n",
      "Running initial clustering\n",
      "Found 2 clusters\n",
      "Permuting 999 times...\n",
      "[                                                            ]   0.00%  |\n",
      "Computing cluster p-values\n",
      "Step-down-in-jumps iteration #1 found 2 clusters to exclude from subsequent iterations\n",
      "Permuting 999 times...\n",
      "[                                                            ]   0.00%  |\n",
      "Computing cluster p-values\n",
      "Step-down-in-jumps iteration #2 found 0 additional clusters to exclude from subsequent iterations\n",
      "Done.\n"
     ]
    }
   ],
   "source": [
    "from scipy import stats\n",
    "from mne.stats import permutation_cluster_1samp_test\n",
    "\n",
    "n_jobs = 2  # nb of parallel jobs\n",
    "\n",
    "channel = 'EEG065'\n",
    "idx = contrast.ch_names.index(channel)\n",
    "data = np.array([c.data[idx] for c in contrast_list])\n",
    "\n",
    "n_permutations = 1000  # number of permutations to run\n",
    "\n",
    "# set initial threshold\n",
    "p_initial = 0.001\n",
    "\n",
    "# set family-wise p-value\n",
    "p_thresh = 0.01\n",
    "\n",
    "connectivity = None\n",
    "tail = 0.  # for two sided test\n",
    "\n",
    "# set cluster threshold\n",
    "n_samples = len(data)\n",
    "threshold = -stats.t.ppf(p_initial / (1 + (tail == 0)), n_samples - 1)\n",
    "if np.sign(tail) < 0:\n",
    "    threshold = -threshold\n",
    "\n",
    "cluster_stats = permutation_cluster_1samp_test(\n",
    "    data, threshold=threshold, n_jobs=n_jobs, verbose=True, tail=tail,\n",
    "    step_down_p=0.05, connectivity=connectivity,\n",
    "    n_permutations=n_permutations, seed=42)\n",
    "\n",
    "T_obs, clusters, cluster_p_values, _ = cluster_stats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Visualize results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "times = 1e3 * contrast.times\n",
    "\n",
    "fig, axes = plt.subplots(2, sharex=True)\n",
    "ax = axes[0]\n",
    "ax.plot(times, 1e6 * data.mean(axis=0), label=\"ERP Contrast\")\n",
    "ax.set(title='Channel : ' + channel, ylabel=\"EEG (uV)\", ylim=[-1, 3])\n",
    "ax.legend()\n",
    "\n",
    "ax = axes[1]\n",
    "for i_c, c in enumerate(clusters):\n",
    "    c = c[0]\n",
    "    if cluster_p_values[i_c] < p_thresh:\n",
    "        h1 = ax.axvspan(times[c.start], times[c.stop - 1],\n",
    "                        color='r', alpha=0.3)\n",
    "hf = ax.plot(times, T_obs, 'g')\n",
    "ax.legend((h1,), (u'p < %s' % p_thresh,), loc='upper right', ncol=1)\n",
    "ax.set(xlabel=\"time (ms)\", ylabel=\"T-values\",\n",
    "       ylim=[-10., 10.], xlim=contrast.times[[0, -1]] * 1000)\n",
    "fig.tight_layout(pad=0.5)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    <b>EXERCISE</b>:\n",
    "     <ul>\n",
    "      <li>Do a cluster level stats in sensor space.</li>\n",
    "    </ul>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "stat_fun(H1): min=-13.252473 max=9.820945\n",
      "No disjoint connectivity sets found\n",
      "Running initial clustering\n",
      "Found 10 clusters\n",
      "Permuting 1023 times...\n",
      "[                                                            ]   0.00%  |"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-14-5ea33054f3bc>:24: RuntimeWarning: Provided stat_fun does not treat variables independently. Setting buffer_size to None.\n",
      "  check_disjoint=True, step_down_p=0.05, seed=42)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Computing cluster p-values\n",
      "Step-down-in-jumps iteration #1 found 1 cluster to exclude from subsequent iterations\n",
      "Permuting 1023 times...\n",
      "[                                                            ]   0.00%  |\n",
      "Computing cluster p-values\n",
      "Step-down-in-jumps iteration #2 found 0 additional clusters to exclude from subsequent iterations\n",
      "Done.\n",
      "Good clusters: [9]\n"
     ]
    }
   ],
   "source": [
    "# Assemble the data and run the cluster stats on channel data\n",
    "data = np.array([c.pick_types(meg=False, eeg=True).crop(None, 0.8).data\n",
    "                 for c in contrast_list])\n",
    "\n",
    "connectivity = None\n",
    "tail = 0.  # for two sided test\n",
    "\n",
    "# set cluster threshold\n",
    "p_thresh = 0.01 / (1 + (tail == 0))\n",
    "n_samples = len(data)\n",
    "threshold = -stats.t.ppf(p_thresh, n_samples - 1)\n",
    "if np.sign(tail) < 0:\n",
    "    threshold = -threshold\n",
    "\n",
    "# Make a triangulation between EEG channels locations to\n",
    "# use as connectivity for cluster level stat\n",
    "connectivity = mne.channels.find_ch_connectivity(contrast.info, 'eeg')[0]\n",
    "\n",
    "data = np.transpose(data, (0, 2, 1))  # transpose for clustering\n",
    "\n",
    "cluster_stats = permutation_cluster_1samp_test(\n",
    "    data, threshold=threshold, n_jobs=2, verbose=True, tail=1,\n",
    "    connectivity=connectivity, out_type='indices',\n",
    "    check_disjoint=True, step_down_p=0.05, seed=42)\n",
    "\n",
    "T_obs, clusters, p_values, _ = cluster_stats\n",
    "good_cluster_inds = np.where(p_values < 0.05)[0]\n",
    "\n",
    "print(\"Good clusters: %s\" % good_cluster_inds)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Visualize the spatio-temporal clusters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "from mne.viz import plot_topomap\n",
    "\n",
    "times = contrast.times * 1e3\n",
    "colors = 'r', 'steelblue'\n",
    "linestyles = '-', '--'\n",
    "\n",
    "pos = mne.find_layout(contrast.info).pos\n",
    "\n",
    "T_obs_max = 5.\n",
    "T_obs_min = -T_obs_max\n",
    "\n",
    "# loop over significant clusters\n",
    "for i_clu, clu_idx in enumerate(good_cluster_inds):\n",
    "\n",
    "    # unpack cluster information, get unique indices\n",
    "    time_inds, space_inds = np.squeeze(clusters[clu_idx])\n",
    "    ch_inds = np.unique(space_inds)\n",
    "    time_inds = np.unique(time_inds)\n",
    "\n",
    "    # get topography for T0 stat\n",
    "    T_obs_map = T_obs[time_inds, ...].mean(axis=0)\n",
    "\n",
    "    # get signals at significant sensors\n",
    "    signals = data[..., ch_inds].mean(axis=-1)\n",
    "    sig_times = times[time_inds]\n",
    "\n",
    "    # create spatial mask\n",
    "    mask = np.zeros((T_obs_map.shape[0], 1), dtype=bool)\n",
    "    mask[ch_inds, :] = True\n",
    "\n",
    "    # initialize figure\n",
    "    fig, ax_topo = plt.subplots(1, 1, figsize=(7, 2.))\n",
    "\n",
    "    # plot average test statistic and mark significant sensors\n",
    "    image, _ = plot_topomap(T_obs_map, pos, mask=mask, axes=ax_topo,\n",
    "                            vmin=T_obs_min, vmax=T_obs_max,\n",
    "                            show=False)\n",
    "\n",
    "    # advanced matplotlib for showing image with figure and colorbar\n",
    "    # in one plot\n",
    "    divider = make_axes_locatable(ax_topo)\n",
    "\n",
    "    # add axes for colorbar\n",
    "    ax_colorbar = divider.append_axes('right', size='5%', pad=0.05)\n",
    "    plt.colorbar(image, cax=ax_colorbar, format='%0.1f')\n",
    "    ax_topo.set_xlabel('Averaged t-map\\n({:0.1f} - {:0.1f} ms)'.format(\n",
    "        *sig_times[[0, -1]]\n",
    "    ))\n",
    "\n",
    "    # add new axis for time courses and plot time courses\n",
    "    ax_signals = divider.append_axes('right', size='300%', pad=1.2)\n",
    "    for signal, name, col, ls in zip(signals, ['Contrast'], colors,\n",
    "                                     linestyles):\n",
    "        ax_signals.plot(times, signal * 1e6, color=col,\n",
    "                        linestyle=ls, label=name)\n",
    "\n",
    "    # add information\n",
    "    ax_signals.axvline(0, color='k', linestyle=':', label='stimulus onset')\n",
    "    ax_signals.set_xlim([times[0], times[-1]])\n",
    "    ax_signals.set_xlabel('Time [ms]')\n",
    "    ax_signals.set_ylabel('Amplitude [uV]')\n",
    "\n",
    "    # plot significant time range\n",
    "    ymin, ymax = ax_signals.get_ylim()\n",
    "    ax_signals.fill_betweenx((ymin, ymax), sig_times[0], sig_times[-1],\n",
    "                             color='orange', alpha=0.3)\n",
    "    ax_signals.legend(loc='lower right')\n",
    "    title = 'Cluster #{0} (p < {1:0.3f})'.format(i_clu + 1, p_values[clu_idx])\n",
    "    ax_signals.set(ylim=[ymin, ymax], title=title)\n",
    "\n",
    "    # clean up viz\n",
    "    fig.tight_layout(pad=0.5, w_pad=0)\n",
    "    fig.subplots_adjust(bottom=.05)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Let's do now some source space analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Current object is not active, please select an active object.\n",
      "Current object is not active, please select an active object.\n"
     ]
    }
   ],
   "source": [
    "surfer_kwargs = dict(initial_time=0.170, subject=\"fsaverage\",\n",
    "                     background=\"white\", foreground='black',\n",
    "                     cortex=(\"gray\", -1, 6, True), smoothing_steps=10,\n",
    "                     clim=dict(kind='value', lims=[3, 6, 7]),\n",
    "                     time_viewer=True)\n",
    "\n",
    "def plot_stc(subject, condition, hemi=\"both\", views=\"ventral\"):\n",
    "    \"\"\"Plot source estimates.\"\"\"\n",
    "    stc_fname = os.path.join(stc_path, subject + \"-%s\" % condition)\n",
    "    stc = mne.read_source_estimate(stc_fname)\n",
    "\n",
    "    brain = stc.plot(hemi=hemi, views=views, **surfer_kwargs)\n",
    "    brain.add_text(0.01, 0.9, condition, condition, font_size=15)\n",
    "    brain.add_text(0.8, 0.9, subject, subject, font_size=15)\n",
    "    return brain\n",
    "\n",
    "\n",
    "brain = plot_stc('sub015', 'faces')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Current object is not active, please select an active object.\n",
      "Current object is not active, please select an active object.\n",
      "Current object is not active, please select an active object.\n"
     ]
    }
   ],
   "source": [
    "def plot_average_stc(condition, subjects=subjects, plot=False, hemi=\"both\",\n",
    "                     views=\"ventral\"):\n",
    "    \"\"\"Compute the average stc for a certain condition.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    condition: str\n",
    "        experimental condition (\"faces\", \"famous\" or \"scrambled\")\n",
    "    subjects: list of str\n",
    "        names of subjects considered in the average.\n",
    "    plot: bool\n",
    "        plot\n",
    "    hemi: str\n",
    "        hemisphere to plot (ignored if `plot` == False)\n",
    "    \"\"\"\n",
    "    stc_fnames = [os.path.join(stc_path, subject + \"-%s\" % condition)\n",
    "                  for subject in subjects]\n",
    "    stcs = [mne.read_source_estimate(stc_fname) for stc_fname in stc_fnames]\n",
    "    average_stc = np.mean(stcs)\n",
    "\n",
    "    if plot:\n",
    "        brain = average_stc.plot(hemi=hemi, views=views, **surfer_kwargs)\n",
    "        brain.add_text(0.01, 0.9, condition, condition, font_size=15)\n",
    "    return average_stc\n",
    "\n",
    "conditions = [\"famous\", \"faces\", \"scrambled\"]\n",
    "for condition in conditions:\n",
    "    plot_average_stc(condition, plot=True)"
   ]
  }
 ],
 "metadata": {
  "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
}