6-Group_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Group analysis with MNE-Python\n",
    "\n",
    "\n",
    "The aim of this lecture is to show you how to do group analysis\n",
    "with MNE-Python.\n",
    "\n",
    "    Authors: Britta Westner, Alexandre Gramfort, Denis A. Engemann, Mainak Jas, Hicham Janati\n",
    "\n",
    "    License: BSD (3-clause)\n",
    "   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import os\n",
    "import numpy as np\n",
    "import mne\n",
    "\n",
    "mne.set_log_level('error')\n",
    "\n",
    "# Change the following path to where the folder with the extra data is on your disk:\n",
    "extra_path = os.path.expanduser('~/Documents/teaching/practical_meeg_2022_data/extra_data_mne')\n",
    "evokeds_path = os.path.join(extra_path, 'group_analysis/')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's check if setting the path worked and if the data is there:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ls $evokeds_path"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Datasets\n",
    "\n",
    "We have here the evoked data for all but 1 participant. That one dataset had some issues with the triggers that would have needed fixing to be read into MNE-Python, so here we just discard it (dataset # 10)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "datasets = ['sub-%02d' % ii for ii in range(1, 17) if ii not in [10]]\n",
    "datasets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Let's first look at our different datasets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can look at the datasets one-by-one:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subject = 'sub-02'\n",
    "fname = os.path.join(evokeds_path, ('%s_list-ave.fif' % subject))\n",
    "evokeds = mne.read_evokeds(fname, verbose=False)\n",
    "evokeds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "evokeds[0].plot_joint()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also loop over all datasets to get an overview and to see if there are any datasets\n",
    "that look problematic."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ch_type = 'eeg'  # we pick the EEG channels here\n",
    "conditions = ['famous', 'unfamiliar', 'scrambled']\n",
    "\n",
    "f, axes = plt.subplots(4, 4, figsize=(13, 9), sharex=True, sharey=True)\n",
    "\n",
    "for ax, subject in zip(axes.ravel(), datasets):\n",
    "    evokeds_dict = dict()\n",
    "    fname = os.path.join(evokeds_path, ('%s_list-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(tmin=-0.3, tmax=0.6)\n",
    "    mne.viz.plot_compare_evokeds(evokeds_dict, picks=ch_type, 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 MEG Magnetometers</li>\n",
    "      <li>Do you see any problematic datasets?</li>\n",
    "    </ul>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read all data and compute contrast\n",
    "\n",
    "Next we will read in all datasets and append them to _one list_:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "evokeds_list = []\n",
    "\n",
    "for subject in datasets:\n",
    "    fname = os.path.join(evokeds_path, ('%s_list-ave.fif' % subject))\n",
    "    evokeds = mne.read_evokeds(fname)\n",
    "    evokeds = [ev for ev in evokeds if ev.comment in ['famous', 'scrambled']]\n",
    "    evokeds_list.append(evokeds)\n",
    "\n",
    "evokeds_list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From here, we can compute the contrast, e.g. between _famous_ and _scrambled_. We can also plot this contrast."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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(), datasets, evokeds_list):\n",
    "    contrast = mne.combine_evoked(evokeds, weights=[0.5, -0.5]).crop(-0.3, 0.5)\n",
    "    contrast.comment = 'contrast'\n",
    "    contrast_list.append(contrast)\n",
    "    mne.viz.plot_compare_evokeds(contrast, picks=ch_type, show=False,\n",
    "                                 axes=ax, title=subject)\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "contrast_list"
   ]
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mne.grand_average?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "evoked_gave = mne.grand_average(contrast_list)\n",
    "evoked_gave.plot_joint();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Let's do some statistics - cluster-permutation test on one channel\n",
    "\n",
    "We'll start with a single channel, EEG065. \n",
    "\n",
    "We will run a cluster-permutation test of our contrast against 0.\n",
    "\n",
    "Let's start with the setup."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# imports for the statistics\n",
    "\n",
    "import scipy as sp\n",
    "from mne.stats import permutation_cluster_1samp_test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's get our data into a numpy array so we can feed it to the stats function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "channel = 'EEG065'\n",
    "ch_idx = contrast.ch_names.index(channel)\n",
    "data = np.array([c.data[ch_idx] for c in contrast_list])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's set some parameters for our test:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_jobs = 2  # number of parallel jobs, you might want to set this to 1\n",
    "n_permutations = 5000  # number of permutations to run \n",
    "\n",
    "# specify adjacency of points in the data: \n",
    "# here, no special adjacency is needed because it's only one channel and MNE will now\n",
    "# just assume a regular grid, which is fine for our time dimension\n",
    "adjacency = None  \n",
    "\n",
    "tail = 0.  # we will run a two-sided test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Okay, with that out of the way, let's set our cluster thresholds:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set the p-value we want to conduct our test at\n",
    "p_value = 0.05\n",
    "\n",
    "# Now, let's compute the threshold for our t-value - we need to pass that to the \n",
    "# stats function.\n",
    "# Beware, that this changes for a one-sided test!\n",
    "deg_of_free = len(data) - 1\n",
    "threshold = sp.stats.t.ppf(1 - p_value/ (1 + (tail == 0)), deg_of_free)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have everything to run our statistics!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "cluster_stats = permutation_cluster_1samp_test(\n",
    "    data, \n",
    "    threshold=threshold, \n",
    "    n_jobs=n_jobs, \n",
    "    verbose=True, \n",
    "    tail=tail,\n",
    "    adjacency=adjacency,\n",
    "    n_permutations=n_permutations, \n",
    "    seed=42)\n",
    "\n",
    "T_obs, clusters, cluster_p_values, _ = cluster_stats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Visualize results from one-channel test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set up the figure:\n",
    "fig, axes = plt.subplots(2, sharex=True)\n",
    "\n",
    "# on first axis: plot the data (contrast averaged across datasets)\n",
    "ax = axes[0]\n",
    "ax.plot(contrast.times, data.mean(axis=0), label='ERP Contrast')\n",
    "ax.set(title='Channel : ' + channel, ylabel='fT/cm')\n",
    "ax.legend()\n",
    "\n",
    "# on second axis:\n",
    "ax = axes[1]\n",
    "# enumerate across the clusters\n",
    "for ii, cluster_ii in enumerate(clusters):\n",
    "    c = cluster_ii[0]\n",
    "\n",
    "    # check if the matching p-value is smaller than the threshold:\n",
    "    if cluster_p_values[ii] < p_value:\n",
    "        # if so, mark the stretch of time:\n",
    "        h1 = ax.axvspan(contrast.times[c[0]], \n",
    "                        contrast.times[c[-1] - 1],\n",
    "                        color='r', alpha=0.3)\n",
    "\n",
    "# plot the t-values \n",
    "hf = ax.plot(contrast.times, T_obs, 'g')\n",
    "\n",
    "# set legend, axes, and show\n",
    "ax.legend((h1,), (u'p < %s' % p_value,), loc='upper right', ncol=1)\n",
    "ax.set(xlabel='Time (ms)', ylabel='T-values',\n",
    "       ylim=[-10., 10.], xlim=contrast.times[[0, -1]])\n",
    "fig.tight_layout(pad=0.5)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cluster-permutation test across time and sensors\n",
    "\n",
    "Now let's not choose a sensor, but run the clustering approach across all _magnetometers_!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# here we can use a convenience function for spatio-temporal data:\n",
    "from mne.stats import spatio_temporal_cluster_1samp_test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# let's use magnetometers here:\n",
    "ch_type = 'mag'\n",
    "\n",
    "# We again need to make our data into an array \n",
    "data = np.array([c.copy().pick_types(meg=ch_type).data\n",
    "                 for c in contrast_list])\n",
    "\n",
    "data.shape\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The clustering function tells us:\n",
    "```\n",
    "X : array, shape (n_observations, p[, q], n_vertices)\n",
    "    The data to be clustered. The first dimension should correspond to the\n",
    "    difference between paired samples (observations) in two conditions.\n",
    "    The second, and optionally third, dimensions correspond to the\n",
    "    time or time-frequency data. And, the last dimension should be spatial.\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make sure spatial dimension is last:\n",
    "data = np.transpose(data, (0, 2, 1)) \n",
    "data.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now our parameters again. This time we need to pay close attentiont to the adjacency!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# tail and cluster thresholds\n",
    "tail = 0.  # for two sided test\n",
    "\n",
    "# set cluster threshold (see above for more detail)\n",
    "deg_of_free = len(data) - 1\n",
    "threshold = sp.stats.t.ppf(1 - p_value/ (1 + (tail == 0)), deg_of_free)\n",
    "\n",
    "# Make a triangulation between MEG sensor locations to\n",
    "# use as adjacency informatiom for cluster level stats:\n",
    "adjacency = mne.channels.find_ch_adjacency(contrast.info, ch_type)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "cluster_stats = spatio_temporal_cluster_1samp_test(\n",
    "    data, \n",
    "    threshold=threshold, \n",
    "    n_jobs=2, \n",
    "    verbose=True, \n",
    "    tail=tail,\n",
    "    adjacency=adjacency, \n",
    "    out_type='indices',\n",
    "    check_disjoint=True, \n",
    "    seed=42)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's make the output easier to use for plotting:\n",
    "T_obs, clusters, p_values, _ = cluster_stats\n",
    "good_cluster_inds = np.where(p_values < 0.05)[0]\n",
    "\n",
    "print(\"Good clusters: %i\" % len(good_cluster_inds))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Visualize the spatio-temporal clusters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "from mne.viz import plot_topomap\n",
    "\n",
    "# some setup for plotting\n",
    "colors = 'r', 'steelblue'\n",
    "linestyles = '-', '--'\n",
    "\n",
    "# find the relevant sensors\n",
    "pos = mne.find_layout(contrast.info, ch_type=ch_type).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 per cluster\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 stats and average across time\n",
    "    T_obs_map = T_obs[time_inds, ...].mean(axis=0)\n",
    "\n",
    "    # get signals at significant sensors and average across sensors\n",
    "    signals = data[..., ch_inds].mean(axis=-1)\n",
    "    sig_times = contrast.times[time_inds]\n",
    "\n",
    "    # create a 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",
    "    # how to mark sensors in the cluster\n",
    "    mask_params = dict(marker='.', markerfacecolor='k', markersize=2)\n",
    "\n",
    "    # plot average test statistic and mark significant sensors\n",
    "    sel = mne.pick_types(contrast.info, meg=ch_type)\n",
    "    info = mne.pick_info(contrast.info, sel)\n",
    "    image, _ = plot_topomap(T_obs_map, info, ch_type=ch_type, mask=mask, \n",
    "                            axes=ax_topo, sensors=False,\n",
    "                            mask_params=mask_params,\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, add colorbar \n",
    "    ax_colorbar = divider.append_axes('right', size='5%', pad=0.05)\n",
    "    plt.colorbar(image, cax=ax_colorbar, format='%0.1f')\n",
    "\n",
    "    # label the topoplot\n",
    "    ax_topo.set_xlabel('Averaged t-map\\n({:0.2f} - {:0.2f} ms)'.format(\n",
    "        *sig_times[[0, -1]]\n",
    "    ))\n",
    "\n",
    "    # add a 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(contrast.times, signal, color=col,\n",
    "                        linestyle=ls, label=name)\n",
    "\n",
    "    # mark stimulus onset\n",
    "    ax_signals.axvline(0, color='k', linestyle=':', label='stimulus onset')\n",
    "\n",
    "    # adjust and label axes\n",
    "    ax_signals.set_xlim([contrast.times[0], contrast.times[-1]])\n",
    "    ax_signals.set_xlabel('Time [s]')\n",
    "    ax_signals.set_ylabel('Amplitude')\n",
    "\n",
    "    # plot the significant time range in time course\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 the figure a little\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": [
    "#### Further reading and credit:\n",
    "\n",
    "The code in this notebook is heavily inspired by https://github.com/mne-tools/mne-biomag-group-demo/tree/master/scripts/results/statistics\n",
    "\n",
    "More about sensor level statistics can be found here: https://mne.tools/stable/auto_tutorials/stats-sensor-space/index.html\n",
    "\n",
    "If you want to do stats in source space, check this: https://mne.tools/stable/auto_tutorials/stats-source-space/index.html"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.10.8 ('mne_aix')",
   "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.10.8"
  },
  "vscode": {
   "interpreter": {
    "hash": "96118219f514f0f7c28e51c58abb2aa3e9c527721ceabff83c3dd194a32d9fc1"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}