{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n\n# Computing a covariance matrix\n\nMany methods in MNE, including source estimation and some classification algorithms,\nrequire covariance estimations from the recordings. In this tutorial we cover the basics\nof sensor covariance computations and construct a noise covariance matrix that can be\nused when computing the minimum-norm inverse solution. For more information, see\n`minimum_norm_estimates`.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Authors: The MNE-Python contributors.\n# License: BSD-3-Clause\n# Copyright the MNE-Python contributors." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import mne\nfrom mne.datasets import sample" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Source estimation methods such as MNE require noise estimates from the\nrecordings. In this tutorial we cover the basics of noise covariance and\nconstruct a noise covariance matrix that can be used when computing the\ninverse solution. For more information, see `minimum_norm_estimates`.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "data_path = sample.data_path()\nraw_empty_room_fname = data_path / \"MEG\" / \"sample\" / \"ernoise_raw.fif\"\nraw_empty_room = mne.io.read_raw_fif(raw_empty_room_fname)\nraw_fname = data_path / \"MEG\" / \"sample\" / \"sample_audvis_raw.fif\"\nraw = mne.io.read_raw_fif(raw_fname)\nraw.set_eeg_reference(\"average\", projection=True)\nraw.info[\"bads\"] += [\"EEG 053\"] # bads + 1 more" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The definition of noise depends on the paradigm. In MEG it is quite common\nto use empty room measurements for the estimation of sensor noise. However if\nyou are dealing with evoked responses, you might want to also consider\nresting state brain activity as noise.\nFirst we compute the noise using empty room recording. Note that you can also\nuse only a part of the recording with tmin and tmax arguments. That can be\nuseful if you use resting state as a noise baseline. Here we use the whole\nempty room recording to compute the noise covariance (``tmax=None`` is the\nsame as the end of the recording, see :func:`mne.compute_raw_covariance`).\n\nKeep in mind that you want to match your empty room dataset to your\nactual MEG data, processing-wise. Ensure that filters\nare all the same and if you use ICA, apply it to your empty-room and subject\ndata equivalently. In this case we did not filter the data and\nwe don't use ICA. However, we do have bad channels and projections in\nthe MEG data, and, hence, we want to make sure they get stored in the\ncovariance object.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "raw_empty_room.info[\"bads\"] = [bb for bb in raw.info[\"bads\"] if \"EEG\" not in bb]\nraw_empty_room.add_proj(\n [pp.copy() for pp in raw.info[\"projs\"] if \"EEG\" not in pp[\"desc\"]]\n)\n\nnoise_cov = mne.compute_raw_covariance(raw_empty_room, tmin=0, tmax=None)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now that you have the covariance matrix in an MNE-Python object you can\nsave it to a file with :func:`mne.write_cov`. Later you can read it back\nusing :func:`mne.read_cov`.\n\nYou can also use the pre-stimulus baseline to estimate the noise covariance.\nFirst we have to construct the epochs. When computing the covariance, you\nshould use baseline correction when constructing the epochs. Otherwise the\ncovariance matrix will be inaccurate. In MNE this is done by default, but\njust to be sure, we define it here manually.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "events = mne.find_events(raw)\nepochs = mne.Epochs(\n raw,\n events,\n event_id=1,\n tmin=-0.2,\n tmax=0.5,\n baseline=(-0.2, 0.0),\n decim=3, # we'll decimate for speed\n verbose=\"error\",\n) # and ignore the warning about aliasing" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that this method also attenuates any activity in your\nsource estimates that resemble the baseline, if you like it or not.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "noise_cov_baseline = mne.compute_covariance(epochs, tmax=0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plot the covariance matrices\n\nTry setting proj to False to see the effect. Notice that the projectors in\nepochs are already applied, so ``proj`` parameter has no effect.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "noise_cov.plot(raw_empty_room.info, proj=True)\nnoise_cov_baseline.plot(epochs.info, proj=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n## How should I regularize the covariance matrix?\n\nThe estimated covariance can be numerically\nunstable and tends to induce correlations between estimated source amplitudes\nand the number of samples available. The MNE manual therefore suggests to\nregularize the noise covariance matrix (see\n`cov_regularization_math`), especially if only few samples are\navailable. Unfortunately it is not easy to tell the effective number of\nsamples, hence, to choose the appropriate regularization.\nIn MNE-Python, regularization is done using advanced regularization methods\ndescribed in :footcite:t:`EngemannGramfort2015`. For this the ``'auto'`` option\ncan be used. With this option, cross-validation will be used to learn the\noptimal regularization:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "noise_cov_reg = mne.compute_covariance(epochs, tmax=0.0, method=\"auto\", rank=None)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This procedure evaluates the noise covariance quantitatively by how well it\nwhitens the data using the\nnegative log-likelihood of unseen data. The final result can also be visually\ninspected.\nUnder the assumption that the baseline does not contain a systematic signal\n(time-locked to the event of interest), the whitened baseline signal should\nbe follow a multivariate Gaussian distribution, i.e.,\nwhitened baseline signals should be between -1.96 and 1.96 at a given time\nsample.\nBased on the same reasoning, the expected value for the :term:`global field\npower (GFP) ` is 1 (calculation of the GFP should take into account the\ntrue degrees of freedom, e.g. ``ddof=3`` with 2 active SSP vectors):\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "evoked = epochs.average()\nevoked.plot_white(noise_cov_reg, time_unit=\"s\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This plot displays both, the whitened evoked signals for each channels and\nthe whitened :term:`GFP`. The numbers in the GFP panel represent the\nestimated rank of the data, which amounts to the effective degrees of freedom\nby which the squared sum across sensors is divided when computing the\nwhitened :term:`GFP`. The whitened :term:`GFP` also helps detecting spurious\nlate evoked components which can be the consequence of over- or\nunder-regularization.\n\nNote that if data have been processed using signal space separation\n(SSS) :footcite:`TauluEtAl2005`,\ngradiometers and magnetometers will be displayed jointly because both are\nreconstructed from the same SSS basis vectors with the same numerical rank.\nThis also implies that both sensor types are not any longer statistically\nindependent.\nThese methods for evaluation can be used to assess model violations.\nAdditional\nintroductory materials can be found [here](https://goo.gl/ElWrxe).\n\nFor expert use cases or debugging the alternative estimators can also be\ncompared (see `ex-evoked-whitening`):\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "noise_covs = mne.compute_covariance(\n epochs, tmax=0.0, method=(\"empirical\", \"shrunk\"), return_estimators=True, rank=None\n)\nevoked.plot_white(noise_covs, time_unit=\"s\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This will plot the whitened evoked for the optimal estimator and display the\n:term:`GFP` for all estimators as separate lines in the related panel.\n\nFinally, let's have a look at the difference between empty room and\nevent related covariance, hacking the ``\"method\"`` option so that their types\nare shown in the legend of the plot.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "evoked_meg = evoked.copy().pick(\"meg\")\nnoise_cov[\"method\"] = \"empty_room\"\nnoise_cov_baseline[\"method\"] = \"baseline\"\nevoked_meg.plot_white([noise_cov_baseline, noise_cov], time_unit=\"s\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Based on the negative log-likelihood, the baseline covariance\nseems more appropriate. Improper regularization can lead to overestimation of\nsource amplitudes, see :footcite:p:`EngemannGramfort2015` for more\ninformation and examples.\n\n## References\n\n.. footbibliography::\n\n" ] } ], "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.12.2" } }, "nbformat": 4, "nbformat_minor": 0 }