7-Bonus-Temporal_Decoding.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Decoding MEG/EEG data\n",
    "### -- A bonus episode to PracticalMEEG 2022 --\n",
    "\n",
    "`Authors: Britta Westner, Alexandre Gramfort`\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see if we can do some _single trial_ decoding of our data (dataset #1)! \n",
    "\n",
    "For that, we will use the MNE-Python decoding module, which makes use of the extensive machine learning library and toolbox [scikit-learn](https://scikit-learn.org/stable/).\n",
    "\n",
    "`\n",
    "Reference:\n",
    "Scikit-learn: Machine Learning in Python,\n",
    "Pedregosa et al., JMLR 12, pp. 2825-2830, 2011.\n",
    "`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import os\n",
    "import matplotlib.pyplot as plt\n",
    "import mne\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load the epochs we saved"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mne.set_log_level('error')\n",
    "\n",
    "# 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)\n",
    "epochs.crop(-0.3, 0.4)\n",
    "epochs.pick_types('grad')  # let's only keep the gradiometers this time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compute the ERFs and look at the contrast between _face_ and _scrambled_ responses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "evoked_face = epochs['face'].average()\n",
    "evoked_scrambled = epochs['scrambled'].average()\n",
    "evoked_contrast = mne.combine_evoked([evoked_face, evoked_scrambled],\n",
    "                                     [0.5, -0.5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's check the signal ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "evoked_face.plot();\n",
    "evoked_scrambled.plot();\n",
    "evoked_contrast.plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... and plot some topographies:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vmax = 55\n",
    "evoked_face.plot_topomap(contours=0, vmin=0, vmax=vmax);\n",
    "evoked_scrambled.plot_topomap(contours=0, vmin=0, vmax=vmax);\n",
    "\n",
    "vmax = 35\n",
    "evoked_contrast.plot_topomap(contours=0, vmin=0, vmax=vmax);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Let's prepare our data for classification"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we equalize the trials per condition, so 50% is actually our chance level:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "epochs.equalize_event_counts(['face', 'scrambled'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For classification, we first need a _response_ vector. In our case, that is a vector that describes for each trial, whether the trial was \n",
    "a _face_ or a _scrambled_ trial. \n",
    "\n",
    "We set _face_ to 1, and _scrambled_ to 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# initialize with zeros and then set faces to 1\n",
    "y = np.zeros(len(epochs.events), dtype=int)\n",
    "y[epochs.events[:, 2] < 17] = 1  # 1 means face\n",
    "\n",
    "y.size # check the length"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also need an array that contains our data (_observations x channels x time_):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = epochs.get_data()\n",
    "X.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Let's see if we can classify the data based on single trials"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we use a logistic regression model and import this from scikit learn:\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "\n",
    "# We then set the parameters (regularization value and solver):\n",
    "logreg = LogisticRegression(C=1e6, solver='liblinear')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's use the `mne.decoding` module to manage our classification:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.pipeline import make_pipeline\n",
    "from mne.decoding import Scaler, Vectorizer, cross_val_multiscore\n",
    "\n",
    "# We make a pipeline object that specifies how to handle the data:\n",
    "clf = make_pipeline(Scaler(epochs.info),\n",
    "                    Vectorizer(),\n",
    "                    logreg)\n",
    "\n",
    "# Now let's run the classification and score based on a 5-fold-cross-validation:\n",
    "scores = cross_val_multiscore(clf, X, y, cv=5, n_jobs=1)\n",
    "\n",
    "# Let's compute the mean scores across cross-validation splits:\n",
    "score = np.mean(scores, axis=0)\n",
    "print('Spatio-temporal: %0.1f%%' % (100 * score,))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Decoding over time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Often more interesting for electrophysiological data is to keep the time dimension and to slide\n",
    "the decoding model across time:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import StandardScaler\n",
    "from mne.decoding import SlidingEstimator\n",
    "\n",
    "# We prepare our pipeline again:\n",
    "clf = make_pipeline(StandardScaler(), logreg)\n",
    "\n",
    "# But now wrap this pipeline in a sliding estimator:\n",
    "time_decod = SlidingEstimator(clf, n_jobs=1, scoring='roc_auc', verbose=True)\n",
    "\n",
    "# Let's score it:\n",
    "scores = cross_val_multiscore(time_decod, X, y, cv=5, n_jobs=1)\n",
    "\n",
    "# Mean scores across cross-validation splits:\n",
    "scores = np.mean(scores, axis=0)\n",
    "\n",
    "# Let's see the dimesion of our scores:\n",
    "scores.shape\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We kept the time dimension - by fitting one model per time-point. Let's plot our scores!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# prepare plot\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "# plot scores and chance level\n",
    "ax.plot(epochs.times, scores, label='score')\n",
    "ax.axhline(.5, color='k', linestyle='--', label='chance')\n",
    "ax.axvline(.0, color='k', linestyle='-')  # mark time 0\n",
    "\n",
    "# set axis labels, legend, title\n",
    "ax.set_xlabel('Times')\n",
    "ax.set_ylabel('AUC')  # Area Under the Curve\n",
    "ax.legend()\n",
    "\n",
    "ax.set_title('Sensor space decoding')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "scrolled": true
   },
   "source": [
    "## Further reading\n",
    "\n",
    "For more details see: https://mne.tools/stable/auto_tutorials/machine-learning/plot_sensors_decoding.html\n",
    "\n",
    "and this book chapter:\n",
    "```\n",
    "Jean-Rémi King, Laura Gwilliams, Chris Holdgraf, Jona Sassenhagen, Alexandre Barachant, Denis Engemann, Eric Larson, Alexandre Gramfort. Encoding and Decoding Neuronal Dynamics: Methodological Framework to Uncover the Algorithms of Cognition. 2018. https://hal.archives-ouvertes.fr/hal-01848442/\n",
    "```"
   ]
  }
 ],
 "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
}