{ "nbformat_minor": 0, "nbformat": 4, "cells": [ { "execution_count": null, "cell_type": "code", "source": [ "%matplotlib inline" ], "outputs": [], "metadata": { "collapsed": false } }, { "source": [ "\n# Demonstrate impact of whitening on source estimates\n\n\nThis example demonstrates the relationship between the noise covariance\nestimate and the MNE / dSPM source amplitudes. It computes source estimates for\nthe SPM faces data and compares proper regularization with insufficient\nregularization based on the methods described in [1]_. The example demonstrates\nthat improper regularization can lead to overestimation of source amplitudes.\nThis example makes use of the previous, non-optimized code path that was used\nbefore implementing the suggestions presented in [1]_. Please do not copy the\npatterns presented here for your own analysis, this is example is purely\nillustrative.\n\n

Note

This example does quite a bit of processing, so even on a\n fast machine it can take a couple of minutes to complete.

\n\nReferences\n----------\n.. [1] Engemann D. and Gramfort A. (2015) Automated model selection in\n covariance estimation and spatial whitening of MEG and EEG signals,\n vol. 108, 328-342, NeuroImage.\n\n" ], "cell_type": "markdown", "metadata": {} }, { "execution_count": null, "cell_type": "code", "source": [ "# Author: Denis A. Engemann \n#\n# License: BSD (3-clause)\n\nimport os\nimport os.path as op\n\nimport numpy as np\nfrom scipy.misc import imread\nimport matplotlib.pyplot as plt\n\nimport mne\nfrom mne import io\nfrom mne.datasets import spm_face\nfrom mne.minimum_norm import apply_inverse, make_inverse_operator\nfrom mne.cov import compute_covariance\n\nprint(__doc__)" ], "outputs": [], "metadata": { "collapsed": false } }, { "source": [ "Get data\n\n" ], "cell_type": "markdown", "metadata": {} }, { "execution_count": null, "cell_type": "code", "source": [ "data_path = spm_face.data_path()\nsubjects_dir = data_path + '/subjects'\n\nraw_fname = data_path + '/MEG/spm/SPM_CTF_MEG_example_faces%d_3D.ds'\n\nraw = io.read_raw_ctf(raw_fname % 1) # Take first run\n# To save time and memory for this demo, we'll just use the first\n# 2.5 minutes (all we need to get 30 total events) and heavily\n# resample 480->60 Hz (usually you wouldn't do either of these!)\nraw = raw.crop(0, 150.).load_data().resample(60, npad='auto')\n\npicks = mne.pick_types(raw.info, meg=True, exclude='bads')\nraw.filter(1, None, n_jobs=1)\n\nevents = mne.find_events(raw, stim_channel='UPPT001')\n\nevent_ids = {\"faces\": 1, \"scrambled\": 2}\ntmin, tmax = -0.2, 0.5\nbaseline = None # no baseline as high-pass is applied\nreject = dict(mag=3e-12)\n\n# Make source space\n\ntrans = data_path + '/MEG/spm/SPM_CTF_MEG_example_faces1_3D_raw-trans.fif'\nsrc = mne.setup_source_space('spm', fname=None, spacing='oct6',\n subjects_dir=subjects_dir, add_dist=False)\nbem = data_path + '/subjects/spm/bem/spm-5120-5120-5120-bem-sol.fif'\nforward = mne.make_forward_solution(raw.info, trans, src, bem)\nforward = mne.convert_forward_solution(forward, surf_ori=True)\ndel src\n\n# inverse parameters\nconditions = 'faces', 'scrambled'\nsnr = 3.0\nlambda2 = 1.0 / snr ** 2\nmethod = 'dSPM'\nclim = dict(kind='value', lims=[0, 2.5, 5])" ], "outputs": [], "metadata": { "collapsed": false } }, { "source": [ "Estimate covariances\n\n" ], "cell_type": "markdown", "metadata": {} }, { "execution_count": null, "cell_type": "code", "source": [ "samples_epochs = 5, 15,\nmethod = 'empirical', 'shrunk'\ncolors = 'steelblue', 'red'\n\nevokeds = list()\nstcs = list()\nmethods_ordered = list()\nfor n_train in samples_epochs:\n # estimate covs based on a subset of samples\n # make sure we have the same number of conditions.\n events_ = np.concatenate([events[events[:, 2] == id_][:n_train]\n for id_ in [event_ids[k] for k in conditions]])\n epochs_train = mne.Epochs(raw, events_, event_ids, tmin, tmax, picks=picks,\n baseline=baseline, preload=True, reject=reject)\n epochs_train.equalize_event_counts(event_ids)\n assert len(epochs_train) == 2 * n_train\n\n noise_covs = compute_covariance(\n epochs_train, method=method, tmin=None, tmax=0, # baseline only\n return_estimators=True) # returns list\n # prepare contrast\n evokeds = [epochs_train[k].average() for k in conditions]\n del epochs_train, events_\n # do contrast\n\n # We skip empirical rank estimation that we introduced in response to\n # the findings in reference [1] to use the naive code path that\n # triggered the behavior described in [1]. The expected true rank is\n # 274 for this dataset. Please do not do this with your data but\n # rely on the default rank estimator that helps regularizing the\n # covariance.\n stcs.append(list())\n methods_ordered.append(list())\n for cov in noise_covs:\n inverse_operator = make_inverse_operator(evokeds[0].info, forward,\n cov, loose=0.2, depth=0.8,\n rank=274)\n stc_a, stc_b = (apply_inverse(e, inverse_operator, lambda2, \"dSPM\",\n pick_ori=None) for e in evokeds)\n stc = stc_a - stc_b\n methods_ordered[-1].append(cov['method'])\n stcs[-1].append(stc)\n del inverse_operator, evokeds, cov, noise_covs, stc, stc_a, stc_b\ndel raw, forward # save some memory" ], "outputs": [], "metadata": { "collapsed": false } }, { "source": [ "Show the resulting source estimates\n\n" ], "cell_type": "markdown", "metadata": {} }, { "execution_count": null, "cell_type": "code", "source": [ "fig, (axes1, axes2) = plt.subplots(2, 3, figsize=(9.5, 6))\n\n\ndef brain_to_mpl(brain):\n \"\"\"convert image to be usable with matplotlib\"\"\"\n tmp_path = op.abspath(op.join(op.curdir, 'my_tmp'))\n brain.save_imageset(tmp_path, views=['ven'])\n im = imread(tmp_path + '_ven.png')\n os.remove(tmp_path + '_ven.png')\n return im\n\n\nfor ni, (n_train, axes) in enumerate(zip(samples_epochs, (axes1, axes2))):\n # compute stc based on worst and best\n ax_dynamics = axes[1]\n for stc, ax, method, kind, color in zip(stcs[ni],\n axes[::2],\n methods_ordered[ni],\n ['best', 'worst'],\n colors):\n brain = stc.plot(subjects_dir=subjects_dir, hemi='both', clim=clim,\n initial_time=0.175)\n\n im = brain_to_mpl(brain)\n brain.close()\n del brain\n ax.axis('off')\n ax.get_xaxis().set_visible(False)\n ax.get_yaxis().set_visible(False)\n ax.imshow(im)\n ax.set_title('{0} ({1} epochs)'.format(kind, n_train * 2))\n\n # plot spatial mean\n stc_mean = stc.data.mean(0)\n ax_dynamics.plot(stc.times * 1e3, stc_mean,\n label='{0} ({1})'.format(method, kind),\n color=color)\n # plot spatial std\n stc_var = stc.data.std(0)\n ax_dynamics.fill_between(stc.times * 1e3, stc_mean - stc_var,\n stc_mean + stc_var, alpha=0.2, color=color)\n\n # signal dynamics worst and best\n ax_dynamics.set_title('{0} epochs'.format(n_train * 2))\n ax_dynamics.set_xlabel('Time (ms)')\n ax_dynamics.set_ylabel('Source Activation (dSPM)')\n ax_dynamics.set_xlim(tmin * 1e3, tmax * 1e3)\n ax_dynamics.set_ylim(-3, 3)\n ax_dynamics.legend(loc='upper left', fontsize=10)\n\nfig.subplots_adjust(hspace=0.4, left=0.03, right=0.98, wspace=0.07)\nfig.canvas.draw()\nfig.show()" ], "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" } } } }