1-Preprocessing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Preprocessing data 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 start out with loading the packages we need. These include `matplotlib` for plotting, `os` for path management, `numpy` for numerical computations, and of course `mne`.\n",
    "We also use matplotlib magic to ask for figure to be plotted inline. "
   ]
  },
  {
   "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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's double check your MNE-Python version. This should give back 1.2.2 or 1.2.3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mne.__version__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We set the log-level of MNE-Python to 'warning' so the output is less verbose:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mne.set_log_level('warning')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Help!\n",
    "\n",
    "Remember, if you need help just ask ... the machine!\n",
    "Let's see how to get the docstring information for a function - here, the function `pick_types`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mne.pick_types?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set the path to the data"
   ]
  },
  {
   "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",
    "raw_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_meg.fif')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(raw_fname)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use `bash` to verify the path is really there - if this gives an error, chances are you made a typo in the path!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ls $raw_fname"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Access and read the raw data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mne.io.read_raw_fif?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw = mne.io.read_raw_fif(raw_fname, preload=False)\n",
    "print(raw)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For general info on importing data you can check:\n",
    "- for MEG: https://mne.tools/stable/auto_tutorials/io/plot_10_reading_meg_data.html\n",
    "- for EEG: https://mne.tools/stable/auto_tutorials/io/plot_20_reading_eeg_data.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Understand your data file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Now let's look at the measurement info. It can give details about:\n",
    "\n",
    "   - sampling rate\n",
    "   - filtering parameters\n",
    "   - available channel types\n",
    "   - bad channels\n",
    "   - etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(raw.info)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    <b>Exercise</b>:\n",
    "     <ul>\n",
    "    <li>How many channels do you have for each type of sensors?</li>\n",
    "    <li>What is the sampling frequency?</li>\n",
    "    <li>Have the data been filtered?</li>\n",
    "    <li>What is the frequency of the line noise?</li>\n",
    "    <li>Is there any bad channel?</li>\n",
    "    </ul>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A closer look at the info dictionary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "raw.info is just a dictionary:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "isinstance(raw.info, dict)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So we can access its elements this way:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw.info['sfreq']  # Sampling frequency"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw.info['bads']  # list of marked bad channels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw.info['line_freq']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A closer look at the channels\n",
    "Next let's see what channels are present. It is available via the `raw.ch_names` attribute."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "type(raw.ch_names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw.ch_names[:10]  # this prints the first ten channels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can index it as a list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw.ch_names[42]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also query the channel type of a specific channel:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "channel_type = mne.io.pick.channel_type(raw.info, 75)\n",
    "print('Channel #75 is of type:', channel_type)  # print this out in a neat way\n",
    "\n",
    "channel_type = mne.io.pick.channel_type(raw.info, 320)\n",
    "print('Channel #320 is of type:', channel_type)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The info also contains all the details about the sensors (type, locations, coordinate frame etc.) in `chs`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(raw.info['chs'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "type(raw.info['chs'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw.info['chs'][0]  # check the first channel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw.info['chs'][330]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we know that there is EEG and MEG channels in the data, we can plot both separately:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw.plot_sensors(kind='topomap', ch_type='grad');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw.plot_sensors(kind='topomap', ch_type='eeg');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting channel types and re-referencing\n",
    "\n",
    "Some channels are wrongly defined as EEG in the file. \n",
    "Two of these are EOG (EEG061 and EEG062) and EEG063 is actually an ECG channel. EEG064 was recording but not connected to anything, so we'll make it `'misc'` instead. \n",
    "We will now set the channel types for those wrongly classified channels. This will be useful for automatic artifact rejection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw.set_channel_types?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw.set_channel_types({'EEG061': 'eog',  # actually EOG not EEG\n",
    "                       'EEG062': 'eog',  # actually EOG not EEG\n",
    "                       'EEG063': 'ecg',  # actually ECG not EEG\n",
    "                       'EEG064': 'misc'})  # EEG064 free-floating electrode\n",
    "\n",
    "# we also rename the EOG and ECG channels:\n",
    "raw.rename_channels({'EEG061': 'EOG061',\n",
    "                     'EEG062': 'EOG062',\n",
    "                     'EEG063': 'ECG063'})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw.info"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw.plot_sensors(kind='topomap', ch_type='eeg');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After we have fixed the channels, we can now compute an average reference for the EEG. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# For setting the reference, we have to load the data into memory:\n",
    "raw.load_data()\n",
    "print(raw.info['custom_ref_applied'])  # let's see if there is a reference applied"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# now let's re-reference\n",
    "raw.set_eeg_reference(ref_channels='average', projection=False)\n",
    "print(raw.info['projs'])  # not added as a projection\n",
    "print(raw.info['custom_ref_applied'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Accessing the data\n",
    "\n",
    "To access the data just use the `[]` syntax as to access any element of a list, dict etc. Note that `raw[]` returns two things: the data and the times array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "start, stop = 0, 10\n",
    "data, times = raw[:, start:stop]  # fetch all channels and the first 10 time points\n",
    "print(data.shape)\n",
    "print(times.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "times  # always starts at 0 by convention"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that `raw[]` returns both the data and the times array."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Resampling the data\n",
    "\n",
    "We will now change the sampling frequency of the data to speed up the computations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw.load_data()  # load data into memory\n",
    "raw.resample(300)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And let's remove unecessary channels - some empty stimulus channels, misc. channels, and HPI channels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw.drop_channels?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "to_drop = ['STI201', 'STI301', 'MISC201', 'MISC202', 'MISC203',\n",
    "           'MISC204', 'MISC205', 'MISC206', 'MISC301', 'MISC302',\n",
    "           'MISC303', 'MISC304', 'MISC305', 'MISC306', 'CHPI001',\n",
    "           'CHPI002', 'CHPI003', 'CHPI004', 'CHPI005', 'CHPI006',\n",
    "           'CHPI007', 'CHPI008', 'CHPI009']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw.drop_channels(to_drop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Filtering the data and plotting raw data\n",
    "\n",
    "We want to filter the data between 0 and 40 Hz using a linear-phase finite-impulse response (FIR) filter."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    <b>EXERCISE</b>:\n",
    "     <ul>\n",
    "      <li>Which parameters do we have to set to achieve this, based on the docstring of the `filter` method?</li>\n",
    "    </ul>\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw.filter?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To see what effect filtering has for our data, let's quickly look at our data first! For full functionality, we ask matplotlib to show the plot in a separate window."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib qt\n",
    "raw.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw.filter(0, 40)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we filtered our data, let's look at it again. Can you spot the difference?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    <b>EXERCISE</b>:\n",
    "     <ul>\n",
    "      <li> Which data changed more due to the filtering: EEG or MEG?</li>\n",
    "      <li> Can you find reasons why?</li>\n",
    "      <li>  Do you see any bad channels?</li>\n",
    "      <li>  Is there any characteristics you can see in the data?</li>\n",
    "       </ul>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For more information on visualizing of raw data, see here: \n",
    "https://mne.tools/0.16/auto_tutorials/plot_visualize_raw.html\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Look at the event structure of the data\n",
    "\n",
    "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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first look at which events are there:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "events = mne.find_events(raw, stim_channel='STI101', verbose=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    <b>Exercise</b>:\n",
    "     <ul>\n",
    "    <li>What is the type of the variable events?</li>\n",
    "    <li>What is the meaning of the 3 columns of events?</li>\n",
    "    <li>How many events of the value 5 are there?\n",
    "    </ul>\n",
    "</div>\n",
    "\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There was a time offset of 34.5 ms in the stimulus presentation. We need to correct the events accordingly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "delay = int(round(0.0345 * raw.info['sfreq']))\n",
    "events[:, 0] = events[:, 0] + delay"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's visualize the paradigm:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "events = events[events[:, 2] < 20]  # take only events with code less than 20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = mne.viz.plot_events(events, raw.info['sfreq']);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For event trigger and conditions we use a Python dictionary with keys that contain \"/\" for grouping sub-conditions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "event_id = {\n",
    "    'face/famous/first': 5,\n",
    "    'face/famous/immediate': 6,\n",
    "    'face/famous/long': 7,\n",
    "    'face/unfamiliar/first': 13,\n",
    "    'face/unfamiliar/immediate': 14,\n",
    "    'face/unfamiliar/long': 15,\n",
    "    'scrambled/first': 17,\n",
    "    'scrambled/immediate': 18,\n",
    "    'scrambled/long': 19,\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = mne.viz.plot_events(events, sfreq=raw.info['sfreq'],\n",
    "                          event_id=event_id);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now re-visit our raw data plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw.plot(event_id=event_id, events=events);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Epoch data and artifact rejection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define epochs parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmin = -0.5  # start of each epoch (500ms before the trigger)\n",
    "tmax = 2.0  # end of each epoch (2000ms after the trigger)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the baseline period:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "baseline = (-0.2, 0)  # means from 200ms before to stim onset (t = 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also pick channels now - MEG, EEG and EOG channels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "picks = mne.pick_types(raw.info, meg=True, eeg=True, eog=True,\n",
    "                       stim=False, exclude='bads')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The easiest (and maybe also most dangerous?) way to clean your data is to define peak-to-peak (amplitude range) rejection parameters for gradiometers, magnetometers and EOG."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">\n",
    "    <b>REMARK</b>:\n",
    "     <ul>\n",
    "    <li>The <a href=\"https://autoreject.github.io/\">autoreject project</a> aims to solve this problem of reject parameter setting. See the <a href=\"https://www.sciencedirect.com/science/article/pii/S1053811917305013\">paper</a> for more info.</li>\n",
    "    </ul>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reject = dict(grad=4000e-13, mag=4e-12, eog=150e-6)  # this can be highly data dependent"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can put all of this together and create epochs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "epochs = mne.Epochs(raw, events, event_id, tmin, tmax, proj=True,\n",
    "                    picks=picks, baseline=baseline,\n",
    "                    reject=reject)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(epochs)  # let's look at some details about the epochs object"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's explicitly drop the epochs we identified as _bad_ through the thresholds we identified above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "epochs.drop_bad()  # remove bad epochs based on reject"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "epochs.load_data()  # load data in memory"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A closer look at artifact rejection\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's have a closer look at the methods of the epochs object.\n",
    "Uncomment the line below and hit ``epochs.<TAB>``"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#epochs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See how epochs were dropped"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "epochs.plot_drop_log();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Wait a second, did we just loose half of our epochs due to EOG???\n",
    "\n",
    "We can probably do better. Let's use the PCA-based signal space projection (SSP) to regress out spatial patterns related to EOG and other offenders, ie., ECG.\n",
    "\n",
    "Here is the workflow, we'll first detect EOG artifacts and visualize their impact. Then we'll compute related spatial patterns to mitigate these artifacts."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# There is a function to create EOG epochs:\n",
    "eog_epochs = mne.preprocessing.create_eog_epochs(raw.copy().filter(1, None))\n",
    "eog_epochs.average().plot_joint();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see where those EOG segments show up in our raw data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib qt\n",
    "raw.plot(events=eog_epochs.events);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compute SSP projections based on the EOG:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "projs_eog, _ = mne.preprocessing.compute_proj_eog(\n",
    "    raw, n_mag=3, n_grad=3, n_eeg=3, average=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "projs_eog  # let's check what they look like"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "mne.viz.plot_projs_topomap(projs_eog, info=epochs.info);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the important question is how many components one should keep? Tip: some of them don't look like clear artifact patterns. \n",
    "\n",
    "The good news is that we don't need to decide __*right*__ now - as you could see the projectors are stored with the data, but inactive at the moment.\n",
    "\n",
    "BUT: let's repeat this procedure for the ECG, i.e. heart beat artifacts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# same business, same issue for ECG\n",
    "ecg_epochs = mne.preprocessing.create_ecg_epochs(raw.copy().filter(1, None))\n",
    "ecg_epochs.average().plot_joint()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that we also face contamination from the cardiac signal... we'll project that out as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "projs_ecg, _ = mne.preprocessing.compute_proj_ecg(\n",
    "    raw, n_mag=3, n_grad=3, n_eeg=3, average=True)\n",
    "mne.viz.plot_projs_topomap(projs_ecg, info=epochs.info);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Apply projections and visualize effect\n",
    "\n",
    "Now let's reverse our previous artifact rejection and apply the projections instead. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we remove the EOG from our rejection here:\n",
    "reject_no_eog = dict(mag=reject['mag'], grad=reject['grad']) \n",
    "\n",
    "epochs_clean = mne.Epochs(raw, events, event_id, tmin, tmax, proj=False,\n",
    "                          picks=picks, baseline=baseline,\n",
    "                          preload=False,\n",
    "                          reject=reject_no_eog)\n",
    "\n",
    "# and then we add the EOG and ECG projs (but we don't apply them yet!)\n",
    "epochs_clean.add_proj(projs_eog + projs_ecg)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at one frontal MEG channel before applying the projections:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "epochs_clean.plot_image(picks='MEG0123', sigma=1.);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's apply the projections to a copy and plot this channel again!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "epochs_proj = epochs_clean.copy().apply_proj()  # apply projs on a copy\n",
    "\n",
    "epochs_proj.plot_image(picks='MEG0123', sigma=1.);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We established earlier, that probably not all projections capture eye blinks and cardiac artifacts. So let's repeat this procedure but only project out the _first_ projection per channel type!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "epochs_clean.del_proj()\n",
    "epochs_clean.add_proj(projs_eog[::3] + projs_ecg[::3])  # only add some SSP projs\n",
    "epochs_proj = epochs_clean.copy().apply_proj()  # apply projs on a copy\n",
    "\n",
    "epochs_proj.plot_image(picks='MEG0123', sigma=1.);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "epochs_proj.info"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This way, we now keep all trials, but remove eye blinks and cardiac artifacts. We will now save the data with the SSP projections _unapplied_."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">\n",
    "    <b>REMARK</b>:\n",
    "     <ul>\n",
    "    <li>MNE keeps SSP projections inside the info and allows to apply them later.</li>\n",
    "    </ul>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Some thoughts on artifact rejection\n",
    "\n",
    "We now tackled the artifacts in this data set by computing SSP projections. There are many other ways to do artifact rejection:\n",
    "\n",
    "- mark artifacts by hand (visual inspection)\n",
    "- use thresholds (which failed on this dataset!)\n",
    "- use ICA\n",
    "- use an automated pipeline, e.g. the <a href=\"https://autoreject.github.io/\">autoreject project</a>\n",
    "- ...\n",
    "\n",
    "The best recommendation is: get to know your (raw) data!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Save Epochs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# let's overwrite\n",
    "epochs = epochs_clean"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The standard way is to save the epochs as a `.fif` file together with all the header data. Epochs are saved with the suffix `-epo.fif`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "epochs_fname = raw_fname.replace('_meg.fif', '-epo.fif')  # create the file name\n",
    "epochs_fname"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "epochs.save(epochs_fname, overwrite=True) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bonus: Visualizing epochs data\n",
    "\n",
    "See [this page](https://mne.tools/stable/auto_tutorials/epochs/20_visualize_epochs.html) for options on how to visualize epochs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We have already looked at the epochs in a stacked plot:\n",
    "\n",
    "epochs_proj.plot_image(picks='EEG065', sigma=1.);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also look at the epochs in a data browser window:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib qt\n",
    "epochs.plot();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "epochs.info"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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
}