{ "nbformat_minor": 0, "nbformat": 4, "cells": [ { "execution_count": null, "cell_type": "code", "source": [ "%matplotlib inline" ], "outputs": [], "metadata": { "collapsed": false } }, { "source": [ "\n==============================================\nCompute effect-matched-spatial filtering (EMS)\n==============================================\n\nThis example computes the EMS to reconstruct the time course of the\nexperimental effect as described in [1]_.\n\nThis technique is used to create spatial filters based on the difference\nbetween two conditions. By projecting the trial onto the corresponding spatial\nfilters, surrogate single trials are created in which multi-sensor activity is\nreduced to one time series which exposes experimental effects, if present.\n\nWe will first plot a trials x times image of the single trials and order the\ntrials by condition. A second plot shows the average time series for each\ncondition. Finally a topographic plot is created which exhibits the temporal\nevolution of the spatial filters.\n\nReferences\n----------\n\n.. [1] Aaron Schurger, Sebastien Marti, and Stanislas Dehaene, \"Reducing\n multi-sensor data to a single time course that reveals experimental\n effects\", BMC Neuroscience 2013, 14:122.\n\n" ], "cell_type": "markdown", "metadata": {} }, { "execution_count": null, "cell_type": "code", "source": [ "# Author: Denis Engemann \n# Jean-Remi King \n#\n# License: BSD (3-clause)\n\nimport numpy as np\nimport matplotlib.pyplot as plt\n\nimport mne\nfrom mne import io, EvokedArray\nfrom mne.datasets import sample\nfrom mne.decoding import EMS, compute_ems\nfrom sklearn.cross_validation import StratifiedKFold\n\nprint(__doc__)\n\ndata_path = sample.data_path()\n\n# Preprocess the data\nraw_fname = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw.fif'\nevent_fname = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw-eve.fif'\nevent_ids = {'AudL': 1, 'VisL': 3}\n\n# Read data and create epochs\nraw = io.read_raw_fif(raw_fname, preload=True)\nraw.filter(0.5, 45, l_trans_bandwidth='auto', h_trans_bandwidth='auto',\n filter_length='auto', phase='zero')\nevents = mne.read_events(event_fname)\n\npicks = mne.pick_types(raw.info, meg='grad', eeg=False, stim=False, eog=True,\n exclude='bads')\n\nepochs = mne.Epochs(raw, events, event_ids, tmin=-0.2, tmax=0.5, picks=picks,\n baseline=None, reject=dict(grad=4000e-13, eog=150e-6),\n preload=True)\nepochs.drop_bad()\nepochs.pick_types(meg='grad')\n\n# Setup the data to use it a scikit-learn way:\nX = epochs.get_data() # The MEG data\ny = epochs.events[:, 2] # The conditions indices\nn_epochs, n_channels, n_times = X.shape" ], "outputs": [], "metadata": { "collapsed": false } }, { "execution_count": null, "cell_type": "code", "source": [ "# Initialize EMS transformer\nems = EMS()\n\n# Initialize the variables of interest\nX_transform = np.zeros((n_epochs, n_times)) # Data after EMS transformation\nfilters = list() # Spatial filters at each time point\n\n# In the original paper, the cross-validation is a leave-one-out. However,\n# we recommend using a Stratified KFold, because leave-one-out tends\n# to overfit and cannot be used to estimate the variance of the\n# prediction within a given fold.\n\nfor train, test in StratifiedKFold(y):\n # In the original paper, the z-scoring is applied outside the CV.\n # However, we recommend to apply this preprocessing inside the CV.\n # Note that such scaling should be done separately for each channels if the\n # data contains multiple channel types.\n X_scaled = X / np.std(X[train])\n\n # Fit and store the spatial filters\n ems.fit(X_scaled[train], y[train])\n\n # Store filters for future plotting\n filters.append(ems.filters_)\n\n # Generate the transformed data\n X_transform[test] = ems.transform(X_scaled[test])\n\n# Average the spatial filters across folds\nfilters = np.mean(filters, axis=0)\n\n# Plot individual trials\nplt.figure()\nplt.title('single trial surrogates')\nplt.imshow(X_transform[y.argsort()], origin='lower', aspect='auto',\n extent=[epochs.times[0], epochs.times[-1], 1, len(X_transform)],\n cmap='RdBu_r')\nplt.xlabel('Time (ms)')\nplt.ylabel('Trials (reordered by condition)')\n\n# Plot average response\nplt.figure()\nplt.title('Average EMS signal')\nmappings = [(key, value) for key, value in event_ids.items()]\nfor key, value in mappings:\n ems_ave = X_transform[y == value]\n plt.plot(epochs.times, ems_ave.mean(0), label=key)\nplt.xlabel('Time (ms)')\nplt.ylabel('a.u.')\nplt.legend(loc='best')\nplt.show()\n\n# Visualize spatial filters across time\nevoked = EvokedArray(filters, epochs.info, tmin=epochs.tmin)\nevoked.plot_topomap()" ], "outputs": [], "metadata": { "collapsed": false } }, { "source": [ "Note that a similar transformation can be applied with `compute_ems`\nHowever, this function replicates Schurger et al's original paper, and thus\napplies the normalization outside a leave-one-out cross-validation, which we\nrecommend not to do.\n\n" ], "cell_type": "markdown", "metadata": {} }, { "execution_count": null, "cell_type": "code", "source": [ "epochs.equalize_event_counts(event_ids)\nX_transform, filters, classes = compute_ems(epochs)" ], "outputs": [], "metadata": { "collapsed": false } } ], "metadata": { "kernelspec": { "display_name": "Python 2", "name": "python2", "language": "python" }, "language_info": { "mimetype": "text/x-python", "nbconvert_exporter": "python", "name": "python", "file_extension": ".py", "version": "2.7.13", "pygments_lexer": "ipython2", "codemirror_mode": { "version": 2, "name": "ipython" } } } }