3-decoding_csp_eeg.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Motor imagery decoding from EEG data using the Common Spatial Pattern (CSP)\n",
    "\n",
    "Authors: Martin Billinger\n",
    "\n",
    "Decoding of motor imagery applied to EEG data decomposed using CSP.\n",
    "Here the classifier is applied to features extracted on CSP filtered signals.\n",
    "\n",
    "See http://en.wikipedia.org/wiki/Common_spatial_pattern and [1]\n",
    "\n",
    "The EEGBCI dataset is documented in [2]\n",
    "The data set is available at PhysioNet [3]\n",
    "\n",
    "[1] Zoltan J. Koles. The quantitative extraction and topographic mapping\n",
    "    of the abnormal components in the clinical EEG. Electroencephalography\n",
    "    and Clinical Neurophysiology, 79(6):440--447, December 1991.\n",
    "\n",
    "[2] Schalk, G., McFarland, D.J., Hinterberger, T., Birbaumer, N.,\n",
    "    Wolpaw, J.R. (2004) BCI2000: A General-Purpose Brain-Computer Interface\n",
    "    (BCI) System. IEEE TBME 51(6):1034-1043\n",
    "\n",
    "[3] Goldberger AL, Amaral LAN, Glass L, Hausdorff JM, Ivanov PCh, Mark RG,\n",
    "    Mietus JE, Moody GB, Peng C-K, Stanley HE. (2000) PhysioBank,\n",
    "    PhysioToolkit, and PhysioNet: Components of a New Research Resource for\n",
    "    Complex Physiologic Signals. Circulation 101(23):e215-e220\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import mne\n",
    "from mne.channels import read_layout\n",
    "from mne.io import concatenate_raws, read_raw_edf\n",
    "from mne.datasets import eegbci\n",
    "from mne.decoding import CSP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set parameters and read data\n",
    "\n",
    "# avoid classification of evoked responses by using epochs that\n",
    "# start 1s after cue onset.\n",
    "tmin, tmax = -1., 4.\n",
    "subject = 1\n",
    "runs = [6, 10, 14]  # motor imagery: hands vs feet\n",
    "\n",
    "raw_fnames = eegbci.load_data(subject, runs)\n",
    "raw_fnames"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw_files = [read_raw_edf(f, preload=True) for f in raw_fnames]\n",
    "raw = concatenate_raws(raw_files)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set channel locations to plot topographies\n",
    "from mne.channels import make_standard_montage\n",
    "raw.rename_channels(lambda x: x.strip('.'))\n",
    "montage = make_standard_montage('standard_1005')\n",
    "raw.set_montage(montage, False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw.info"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw.annotations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib qt\n",
    "raw.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "# strip channel names of \".\" characters\n",
    "raw.rename_channels(lambda x: x.strip('.'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw.ch_names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Apply band-pass filter\n",
    "raw.filter(7., 30., method='iir')\n",
    "\n",
    "# events = find_events(raw, shortest_event=0, stim_channel='STI 014', verbose=True)\n",
    "events, _ = mne.events_from_annotations(raw)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "events[:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "event_id = dict(hands=2, feet=3)\n",
    "\n",
    "picks = mne.pick_types(raw.info, meg=False, eeg=True, stim=False, eog=False,\n",
    "                       exclude='bads')\n",
    "\n",
    "# Read epochs (train will be done only between 1 and 2s)\n",
    "# Testing will be done with a running classifier\n",
    "epochs = mne.Epochs(raw, events, event_id, tmin, tmax, proj=True, picks=picks,\n",
    "                    baseline=None, preload=True)\n",
    "epochs_train = epochs.copy().crop(tmin=1., tmax=2.)\n",
    "labels = epochs.events[:, -1] - 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "epochs['feet'].average().plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Classification with linear discrimant analysis\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.discriminant_analysis import LinearDiscriminantAnalysis\n",
    "from sklearn.model_selection import ShuffleSplit\n",
    "\n",
    "# Assemble a classifier\n",
    "lda = LinearDiscriminantAnalysis()\n",
    "csp = CSP(n_components=4, reg=None, log=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "epochs.get_data().shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define a monte-carlo cross-validation generator (reduce variance):\n",
    "cv = ShuffleSplit(10, test_size=0.2, random_state=42)\n",
    "scores = []\n",
    "X = epochs.get_data()\n",
    "X_train = epochs_train.get_data()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use scikit-learn Pipeline with cross_val_score function\n",
    "from sklearn.pipeline import Pipeline\n",
    "from sklearn.model_selection import cross_val_score\n",
    "clf = Pipeline([('CSP', csp), ('LDA', lda)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scores = cross_val_score(clf, X_train, labels, cv=cv, n_jobs=1)\n",
    "print(scores)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Printing the results\n",
    "class_balance = np.mean(labels == labels[0])\n",
    "class_balance = max(class_balance, 1. - class_balance)\n",
    "print(\"Classification accuracy: %f / Chance level: %f\" % (np.mean(scores),\n",
    "                                                          class_balance))\n",
    "\n",
    "# plot CSP patterns estimated on full data for visualization\n",
    "csp.fit_transform(X, labels)\n",
    "\n",
    "import mne\n",
    "info = epochs.average().info.copy()\n",
    "info['sfreq'] = 1\n",
    "evoked = mne.EvokedArray(csp.patterns_.T, info, tmin=0)\n",
    "\n",
    "csp.plot_patterns(epochs.info, ch_type='eeg',\n",
    "                  units='Patterns (AU)', size=1.5);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Look at performance over time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sfreq = raw.info['sfreq']\n",
    "w_length = int(sfreq * 0.5)   # running classifier: window length\n",
    "w_step = int(sfreq * 0.1)  # running classifier: window step size\n",
    "w_start = np.arange(0, X.shape[2] - w_length, w_step)\n",
    "\n",
    "scores_windows = []\n",
    "\n",
    "for train_idx, test_idx in cv.split(X, labels):\n",
    "    y_train, y_test = labels[train_idx], labels[test_idx]\n",
    "\n",
    "    XX_train = csp.fit_transform(X_train[train_idx], y_train)\n",
    "    XX_test = csp.transform(X_train[test_idx])\n",
    "\n",
    "    # fit classifier\n",
    "    lda.fit(XX_train, y_train)\n",
    "\n",
    "    # running classifier: test classifier on sliding window\n",
    "    score_this_window = []\n",
    "    for n in w_start:\n",
    "        XX_test = csp.transform(X[test_idx][:, :, n:(n + w_length)])\n",
    "        score_this_window.append(lda.score(XX_test, y_test))\n",
    "    scores_windows.append(score_this_window)\n",
    "\n",
    "# Plot scores over time\n",
    "w_times = (w_start + w_length / 2.) / sfreq + epochs.tmin\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(w_times, np.mean(scores_windows, 0), label='Score')\n",
    "plt.axvline(0, linestyle='--', color='k', label='Onset')\n",
    "plt.axhline(0.5, linestyle='-', color='k', label='Chance')\n",
    "plt.xlabel('time (s)')\n",
    "plt.ylabel('classification accuracy')\n",
    "plt.title('Classification score over time')\n",
    "plt.legend(loc='lower right')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    <b>EXERCISES:</b>\n",
    "     <ul>\n",
    "      <li>How is the performance of the CSP affected by the choice of its parameter such as n_components, reg or log? You will plot the performance of the model as a function of n_components.</li>\n",
    "      <li>How you would tune these parameters without the risk of overfitting?</li>\n",
    "    </ul>\n",
    "</div>\n",
    "\n",
    "See documentation: https://mne.tools/stable/generated/mne.decoding.CSP.html#mne.decoding.CSP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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
}