2-temporal_decoding.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multivariate statistics (decoding / MVPA) on MEG/EEG"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Authors : Alexandre Gramfort, Richard Höchenberger\n",
    "\n",
    "See more info on decoding on this page: https://mne.tools/stable/auto_tutorials/machine-learning/plot_sensors_decoding.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, load the mne package:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import mne"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We set the log-level 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": [
    "## Access the raw data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we import the `sample` dataset. It will be downloaded automatically (approx. 2 GB)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import os\n",
    "from mne.datasets import sample\n",
    "data_path = sample.data_path()\n",
    "\n",
    "raw_fname = os.path.join(data_path, 'MEG/sample/sample_audvis_filt-0-40_raw.fif')\n",
    "raw = mne.io.read_raw_fif(raw_fname, preload=True)\n",
    "raw"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "High-pass filter the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "raw.filter(l_freq=1, h_freq=None, verbose=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "print(raw.info)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define epochs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're only interested in `auditory left` and `auditory right`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First extract events:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "events = mne.find_events(raw, stim_channel='STI 014', verbose=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Look at the design in a graphical way. We're only interested in `auditory left` and `auditory right`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "event_id = {'aud_l': 1, 'aud_r': 2}  # event trigger and conditions\n",
    "fig = mne.viz.plot_events(events, sfreq=raw.info['sfreq'],\n",
    "                          first_samp=raw.first_samp, event_id=event_id)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define epochs parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "tmin = -0.1  # start of each epoch (in seconds)\n",
    "tmax = 0.4   # end of each epoch\n",
    "baseline = None  # no baseline correction, data were high-pass filtered\n",
    "\n",
    "reject = dict(eeg=80e-6, eog=40e-6)\n",
    "picks = mne.pick_types(raw.info, eeg=True, meg=True,\n",
    "                       eog=True, stim=False, exclude='bads')\n",
    "\n",
    "epochs = mne.Epochs(raw, events, event_id, tmin, tmax, proj=True,\n",
    "                    picks=picks, baseline=baseline,\n",
    "                    reject=reject, preload=True)  # with preload\n",
    "\n",
    "print(epochs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Look at the ERF and contrast between left and right stimulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "evoked_left = epochs['aud_l'].average()\n",
    "evoked_right = epochs['aud_r'].average()\n",
    "evoked_contrast = mne.combine_evoked([evoked_left, evoked_right],\n",
    "                                     weights=[1, -1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = evoked_left.plot()\n",
    "fig = evoked_right.plot()\n",
    "fig = evoked_contrast.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot some topographies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vmin, vmax = -4, 4  # Colorbar range\n",
    "fig = evoked_left.plot_topomap(ch_type='eeg', contours=0, vmin=vmin, vmax=vmax)\n",
    "fig = evoked_right.plot_topomap(ch_type='eeg', contours=0, vmin=vmin, vmax=vmax)\n",
    "fig = evoked_contrast.plot_topomap(ch_type='eeg', contours=0, vmin=None, vmax=None)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Now let's see if we can classify single trials."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To keep chance level at 50% accuracy, we first equalize the number of epochs in each condition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "epochs.equalize_event_counts(event_id)\n",
    "print(epochs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A classifier takes as input an `X` and returns `y` (0 or 1). Here `X` will be the data at one time point on all gradiometers (hence the term multivariate). We want to train our model to discriminate between the  `auditory left` and the `auditory right` trials.\n",
    "\n",
    "We work with all sensors jointly and try to find a discriminative pattern between the two conditions to predict the class.\n",
    "\n",
    "For classification we will use the scikit-learn package (http://scikit-learn.org/) and MNE functions.\n",
    "\n",
    "Let's first create the response vector, `y`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "y = np.empty(len(epochs.events), dtype=int)\n",
    "idx_auditory_left = epochs.events[:, 2] == event_id['aud_l']\n",
    "idx_auditory_right = epochs.events[:, 2] == event_id['aud_r']\n",
    "y[idx_auditory_left] = 0\n",
    "y[idx_auditory_right] = 1\n",
    "\n",
    "y.size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, the input matrix, `X`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = epochs.copy().pick_types(meg='grad').get_data()\n",
    "X.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "XX = X.reshape(108, -1)\n",
    "XX.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from sklearn.model_selection import StratifiedKFold, cross_val_score\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.pipeline import make_pipeline\n",
    "\n",
    "logreg = LogisticRegression()\n",
    "cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)\n",
    "clf = make_pipeline(StandardScaler(), logreg)\n",
    "\n",
    "scores = cross_val_score(clf, XX, y, cv=cv, scoring='roc_auc')\n",
    "\n",
    "roc_auc_mean = np.mean(scores)\n",
    "roc_auc_std = np.std(scores)\n",
    "\n",
    "print(f'CV scores: {scores}')\n",
    "print(f'Mean ROC AUC = {roc_auc_mean:.3f} (std: {roc_auc_std:.3f})')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## We can do this more simply using the `mne.decoding` module! Let's do some spatio-temporal decoding."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from sklearn.pipeline import make_pipeline\n",
    "from mne.decoding import Scaler, Vectorizer, cross_val_multiscore\n",
    "\n",
    "epochs_decoding = epochs.copy().pick_types(meg='grad')\n",
    "\n",
    "cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)\n",
    "clf = make_pipeline(Scaler(epochs_decoding.info),\n",
    "                    Vectorizer(),\n",
    "                    logreg)\n",
    "\n",
    "X = epochs_decoding.get_data()\n",
    "y = epochs_decoding.events[:, 2]\n",
    "\n",
    "scores = cross_val_multiscore(clf, X, y, cv=cv, scoring='roc_auc')\n",
    "\n",
    "roc_auc_mean = np.mean(scores)\n",
    "roc_auc_std = np.std(scores)\n",
    "\n",
    "print(f'CV scores: {scores}')\n",
    "print(f'Mean ROC AUC = {roc_auc_mean:.3f} (std: {roc_auc_std:.3f})')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    <b>EXERCISES:</b>\n",
    "     <ul>\n",
    "      <li>Why do you get different results from above? what did we change in the model?</li>\n",
    "      <li>How does the choice of cross-validation affect the results? Hint: Change the random_state</li>\n",
    "      <li>Try a different cross-validtion object like scikit-learn StratifiedShuffleSplit</li>\n",
    "      <li>Which sensor types give the best classification scores? EEG, MEG gradiometers, MEG magnetometers?</li>\n",
    "    </ul>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Decoding over time: Comparisons at every single time point."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the previous examples, we have trained a classifier to discriminate between experimentel conditions by using the spatio-temporal patterns of **entire trials**. Consequently, the classifier was (hopefully!) able to predict which activation pattern belonged to which condition. \n",
    "\n",
    "However, an interesting neuroscientific is: **Exactly *when* do the brain signals for two conditions differ?**\n",
    "\n",
    "We can try to answer this question by fitting a classifier **at every single time point.** If the classifier can successfully discriminate between the two conditions, we can conclude that the spatial activation patterns measured by the MEG or EEG sensors differed **at this time point**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import StandardScaler\n",
    "from mne.decoding import SlidingEstimator\n",
    "\n",
    "X = epochs_decoding.get_data()\n",
    "y = epochs_decoding.events[:, 2]\n",
    "\n",
    "\n",
    "clf = make_pipeline(StandardScaler(),\n",
    "                    logreg)\n",
    "\n",
    "time_decod = SlidingEstimator(clf, scoring='roc_auc', n_jobs=1, verbose=True)\n",
    "scores = cross_val_multiscore(time_decod, X, y, cv=5, n_jobs=1)\n",
    "\n",
    "# Mean scores across cross-validation splits, for each time point.\n",
    "mean_scores = np.mean(scores, axis=0)\n",
    "\n",
    "# Mean score across all time points.\n",
    "mean_across_all_times = np.mean(scores)\n",
    "print(f'Mean CV score across all time points: {mean_across_all_times:.3f}')\n",
    "\n",
    "# Plot\n",
    "import matplotlib.pyplot as plt\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(epochs.times, mean_scores, label='score')\n",
    "ax.axhline(0.5, color='k', linestyle='--', label='chance')\n",
    "ax.set_xlabel('Time (s)')\n",
    "ax.set_ylabel('Mean ROC AUC')\n",
    "ax.legend()\n",
    "ax.axvline(0, color='k', linestyle='-')\n",
    "ax.set_title('Sensor space decoding')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scores.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "scrolled": true
   },
   "source": [
    "For more details see: https://martinos.org/mne/stable/auto_tutorials/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/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    <b>EXERCISES:</b>\n",
    "     <ul>\n",
    "      <li>Plot the decoding score over time for the different channel types.</li>\n",
    "      <li>Do a decoding over time on the SPM `face` dataset to see if you can classify `face` vs. `scrambled face`.</li>\n",
    "         <li>Do a generalization over time analysis as explained in the <a href=\"https://mne.tools/stable/auto_tutorials/machine-learning/plot_sensors_decoding.html#temporal-generalization\">documentation on decoding</a>.\n",
    "</li>\n",
    "    </ul>\n",
    "</div>\n",
    "\n",
    "Hints:\n",
    "\n",
    "- Access the `face` dataset via:\n",
    "\n",
    "    ```\n",
    "    from mne.datasets import spm_face\n",
    "    data_path = spm_face.data_path()\n",
    "\n",
    "    raw_fname = os.path.join(data_path, 'MEG/spm/SPM_CTF_MEG_example_faces1_3D.ds')\n",
    "    raw = mne.io.read_raw_ctf(raw_fname, preload=True)\n",
    "    ```\n",
    "\n",
    "- The event IDs are:\n",
    "\n",
    "    ```\n",
    "    event_ids = {\"faces\": 1, \"scrambled\": 2}\n",
    "    ```\n",
    "\n",
    "See this online example for additional hints: https://mne.tools/dev/auto_examples/datasets/spm_faces_dataset.html"
   ]
  }
 ],
 "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.8.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}