2-Time_domain_evoked_responses.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sensor space analysis: Time domain analysis (ERP/ERF) in MNE-Python\n",
    "\n",
    "`\n",
    "Authors:\n",
    "Britta Westner, Alexandre Gramfort, Denis A. Engemann\n",
    "`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup\n",
    "\n",
    "We already know this part from last time:\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import os\n",
    "import numpy as np\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import mne\n",
    "\n",
    "mne.set_log_level('error')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set the path to the data and load the previously saved epochs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should have downloaded the `ds000117-practical` folder. We have to let Python know, where to find this folder on your disk. You will have to adjust the path below to reflect your computer and path structure!\n",
    "You can print the whole path and check the directory to double check it's correct."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Change the following path to where the folder ds000117 is on your disk\n",
    "data_path = os.path.expanduser(\"~/Documents/teaching/practical_meeg_2022_data/ds000117\")\n",
    "\n",
    "epochs_fname = os.path.join(data_path,\n",
    "    'derivatives/meg_derivatives/sub-01/ses-meg/meg/sub-01_ses-meg_task-facerecognition_run-01_proc-sss-epo.fif')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "epochs = mne.read_epochs(epochs_fname, proj=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "epochs.info"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wait, why are the projections \"on\"? We did not apply them before saving - can you find out why the seem to be applied now?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Explore the different conditions in the data\n",
    "\n",
    "We have seen in the last session, that the data has different events, which mark which stimulus was shown to the participants. The event/trigger structure is as follows:\n",
    "- 5, 6, 7: famous faces\n",
    "- 13, 14, 15: unfamiliar faces\n",
    "- 17, 18, 19: scrambled faces\n",
    "\n",
    "We have epoched the data according to these events previously. Let's see how we can use these event codes to work with the conditions!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first look at how we registered the events:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "epochs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now query the data based on the event IDs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "epochs['long']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "epochs['scrambled', 'famous']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "epochs['famous']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now for example use this for plotting the epochs based on condition:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "epochs['famous'].plot_image(picks='MEG2143', sigma=1.);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    <b>EXERCISE</b>:\n",
    "     <ul>\n",
    "      <li>Can you make the same plot for the other two conditions?</li>\n",
    "      <li>Is there anything you notice?</li>   \n",
    "    </ul>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Average the epochs to get ERF/ERP "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "evoked = epochs.average()\n",
    "evoked.apply_proj();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "evoked.info"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at our evoked field and evoked potential!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "evoked.plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Separate the conditions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at the ERP at channel EEG065 for all three conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remember how we can index the epochs:\n",
    "for condition in ['famous', 'unfamiliar', 'scrambled']:\n",
    "    epochs[condition].average().plot(picks='EEG065');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To more easily handle our evoked object, let's make a _list of evoked objects_ for the different conditions we are interested in. Then we can also make it easier for ourselved and compare the ERPs in one plot!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "conditions = ['famous', 'unfamiliar', 'scrambled']\n",
    "evokeds_list = []\n",
    "for k, cond in enumerate(conditions):\n",
    "   evokeds_list.append(epochs[cond].average().apply_proj().crop(-0.5, 1))\n",
    "   evokeds_list[k].comment = cond  # we want to update the name of the condition for plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mne.viz.plot_compare_evokeds(evokeds_list, picks='EEG065');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But what if we want to see _all_ channels?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib qt\n",
    "mne.viz.plot_compare_evokeds(evokeds_list, picks='meg', axes='topo');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    <b>EXERCISE</b>:\n",
    "     <ul>\n",
    "      <li>Try clicking on the plot, what happens?</li>\n",
    "      <li>Can you make the same plots for magnetometers and gradiometers?</li>\n",
    "    </ul>\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot standard deviation across trials\n",
    "\n",
    "It can be helpful to visualize evoked data with confidence intervals or standard deviations. Let's see how!\n",
    "\n",
    "For standard deviations across _trials_, there is no ready-to-use MNE-Python function. The reason: Evoked objects do not store the information to compute this metric (i.e., single trials) anymore!\n",
    "\n",
    "But don't worry: we can plot this ourselves using matplotlib - and gain some knowledge about matplotlib and writing functions on the way!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's write our own function to compute the standard deviations!\n",
    "# Since we have to use this operation for each condition, it's handy to write a function for it:\n",
    "\n",
    "# Everything we need within the function, we have to pass to the function.\n",
    "# In this case, we need the epochs object, and we want to pick this object for the right condition\n",
    "# and for channels, so we also pass condition and picks.\n",
    "def compute_sd(epochs, condition, picks):\n",
    "    \"\"\" Get CI from epochs object for condition.\"\"\"\n",
    "\n",
    "    # this gets the data for the condition and channel as an array:\n",
    "    epochs_data = epochs[condition].get_data(picks=picks)\n",
    "\n",
    "    # compute the mean first\n",
    "    epochs_mean = np.mean(epochs_data, axis=0)\n",
    "\n",
    "    # compute the standard deviation\n",
    "    epochs_std = np.std(epochs_data, axis=0)\n",
    "\n",
    "    # add above and below the mean:\n",
    "    sd_upper = np.squeeze(epochs_mean + epochs_std)\n",
    "    sd_lower = np.squeeze(epochs_mean - epochs_std)\n",
    "\n",
    "    # and return what we need\n",
    "    return sd_lower, sd_upper    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we have cropped the evoked, so we also have to crop the epochs!\n",
    "epochs_crop = epochs.copy().crop(-0.5, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "picks = 'EEG065'  # we plot one channel\n",
    "\n",
    "# let's first add the lines at zero that MNE-Python has\n",
    "plt.axhline(y=0., color='black', linewidth=0.75)\n",
    "plt.axvline(x=0., color='black', linewidth=0.75, linestyle='--')\n",
    "\n",
    "for ii, condition in enumerate(conditions):\n",
    "\n",
    "    # first let's compute the standard deviation with our function\n",
    "    sd_lower, sd_upper = compute_sd(epochs_crop, condition, picks)\n",
    "\n",
    "    # get the evoked data from the evokeds list\n",
    "    evoked_data = evokeds_list[ii].get_data(picks=picks)\n",
    "    times = evokeds_list[ii].times\n",
    "\n",
    "    # plot the evoked\n",
    "    plt.plot(times, np.squeeze(evoked_data), linewidth=2, label=condition)\n",
    "\n",
    "    # plot the standard deviation and add a label once\n",
    "    plt.fill_between(times, sd_lower, sd_upper, alpha=0.3,\n",
    "                     label=('standard dev.' if ii==2 else None))\n",
    "\n",
    "\n",
    "\n",
    "# labels and legend \n",
    "plt.ylabel('V')  # note the difference in scaling to MNE-Python!\n",
    "plt.xlabel('Time (s)')\n",
    "plt.legend(loc='lower left', framealpha=0.3)\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot topographies\n",
    "\n",
    "Next, we want to plot the topographies at the peak. We want to plot this for all channel types _and_ all conditions. To not wind up with 9 individual plots, let's see how to combine plots into one figure!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "# per channel type, we have different color limits for the topography:\n",
    "color_lims = [(-150, 150), (0, 45), (-4.5, 4.5)]\n",
    "channel_types = ['mag', 'grad', 'eeg']\n",
    "sup_titles = ['Magnetometers', 'Gradiometers', 'EEG']\n",
    "\n",
    "# we loop across the channel types first:\n",
    "for ch_type, c_lim, s_title in zip(channel_types, color_lims, sup_titles):\n",
    "\n",
    "    # Here we build our figure. Note that we ask for 4 subplots! That is because our colorbar\n",
    "    # needs its own axis.\n",
    "    fig, axes = plt.subplots(1, 4, figsize=(7, 4), gridspec_kw={'width_ratios': [3, 3, 3, 0.5]})\n",
    "\n",
    "    # Now we can get the peak time across the conditions for this channel type.\n",
    "    # We search for the peak time for up to 150 ms after stimulus onset. We will plot \n",
    "    # the topography for 100 ms around this peak - see the average parameter in the \n",
    "    # plotting function call!\n",
    "    _, peak_time = evoked.get_peak(ch_type=ch_type, tmax=0.15)\n",
    "\n",
    "    for ii, (ax, cond) in enumerate(zip(axes, conditions)):\n",
    "        evokeds_list[ii].plot_topomap(\n",
    "            ch_type=ch_type, times=peak_time, average=0.1,\n",
    "            vlim=c_lim, \n",
    "            colorbar=(True if ii==2 else False),  # only one colorbar\n",
    "            axes=(axes[2::] if ii == 2 else ax),  # last plot needs 2 axes\n",
    "            show=False)\n",
    "\n",
    "        ax.set_title(cond)\n",
    "\n",
    "    fig.suptitle('%s, %.3f s' % (s_title, peak_time))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Write evoked data to disk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    " # note that the file for evoked ends with -ave.fif\n",
    "evoked_fname = epochs_fname.replace('-epo.fif', '-ave.fif')\n",
    "evoked_fname"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "evoked.save(evoked_fname, overwrite=True) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also write multiple conditions into one file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "evoked_list_fname = evoked_fname.replace('-ave.fif', '_list-ave.fif')\n",
    "mne.write_evokeds(evoked_list_fname, evokeds_list, overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reading evoked from disk\n",
    "\n",
    "You can read the data back from the stored `.fif` file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "evokeds_list = mne.read_evokeds(evoked_fname, baseline=(None, 0), proj=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or give the explicit name of the averaged condition when reading in the list:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "evoked_famous = mne.read_evokeds(evoked_list_fname, condition=\"famous\",\n",
    "                           baseline=(None, 0), proj=True)\n",
    "evoked_famous"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Remark:** Did you notice that you can apply some preprocessing on reading the evokeds from disk?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Contrast the conditions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# now we can use the condition \"face\", which consists of famous and unfamiliar\n",
    "evoked_face = epochs['face'].average().apply_proj()\n",
    "evoked_scrambled = epochs['scrambled'].average().apply_proj()\n",
    "\n",
    "# adapt the naming\n",
    "evoked_face.comment = 'face'\n",
    "evoked_scrambled.comment = 'scrambled'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "contrast = mne.combine_evoked([evoked_face, evoked_scrambled], [0.5, -0.5])\n",
    "\n",
    "contrast.comment = 'face - scrambled'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that this combines evokeds by taking into account the number of averaged epochs (to scale the noise variance). Let's check out the number of epochs we averaged over:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(evoked_face.nave) \n",
    "print(evoked_scrambled.nave) \n",
    "print(contrast.nave) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(contrast)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot the contrast together with the ERPs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mne.viz.plot_compare_evokeds([evoked_face, evoked_scrambled, contrast], picks='EEG065');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    <b>EXERCISE</b>:\n",
    "     <ul>\n",
    "      <li>Can you plot the difference between famous and unfamiliar faces for an occipital MEG sensor?</li>\n",
    "    </ul>\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot a 3D field map\n",
    "\n",
    "We can also plot the topography in 3D. Since MNE-Python also plots the head surface here, we need a `-trans.fif` file, that tells us the _transform_ between the MEG device coordinate system and the head surface coordinate system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get the path of the transform file\n",
    "# usually, this needs to be created by manual coregistration, but we have it all ready!\n",
    "extra_path = os.path.expanduser(\"~/Documents/teaching/practical_meeg_2022_data/extra_data_mne/\")\n",
    "trans_fname = os.path.join(extra_path, 'source_recon/sub-01/sub-01-trans.fif')\n",
    "\n",
    "# we also need the Freesurfer subjects direcory:\n",
    "subjects_dir = os.path.join(extra_path, 'freesurfer')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's double check it exists:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ls $trans_fname"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# first make the map\n",
    "maps = mne.make_field_map(evoked_face, trans=trans_fname, ch_type='meg',\n",
    "                          subject='sub-01', subjects_dir=subjects_dir)\n",
    "evoked_face.plot_field(maps, time=0.1)"
   ]
  }
 ],
 "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"
   }
  },
  "widgets": {
   "state": {},
   "version": "1.1.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}