{ "nbformat_minor": 0, "nbformat": 4, "cells": [ { "execution_count": null, "cell_type": "code", "source": [ "%matplotlib inline" ], "outputs": [], "metadata": { "collapsed": false } }, { "source": [ "\n\n# Non-parametric between conditions cluster statistic on single trial power\n\n\nThis script shows how to compare clusters in time-frequency\npower estimates between conditions. It uses a non-parametric\nstatistical procedure based on permutations and cluster\nlevel statistics.\n\nThe procedure consists in:\n\n - extracting epochs for 2 conditions\n - compute single trial power estimates\n - baseline line correct the power estimates (power ratios)\n - compute stats to see if the power estimates are significantly different\n between conditions.\n\n\n" ], "cell_type": "markdown", "metadata": {} }, { "execution_count": null, "cell_type": "code", "source": [ "# Authors: Alexandre Gramfort \n#\n# License: BSD (3-clause)\n\nimport numpy as np\nimport matplotlib.pyplot as plt\n\nimport mne\nfrom mne.time_frequency import tfr_morlet\nfrom mne.stats import permutation_cluster_test\nfrom mne.datasets import sample\n\nprint(__doc__)" ], "outputs": [], "metadata": { "collapsed": false } }, { "source": [ "Set parameters\n\n" ], "cell_type": "markdown", "metadata": {} }, { "execution_count": null, "cell_type": "code", "source": [ "data_path = sample.data_path()\nraw_fname = data_path + '/MEG/sample/sample_audvis_raw.fif'\nevent_fname = data_path + '/MEG/sample/sample_audvis_raw-eve.fif'\ntmin, tmax = -0.2, 0.5\n\n# Setup for reading the raw data\nraw = mne.io.read_raw_fif(raw_fname)\nevents = mne.read_events(event_fname)\n\ninclude = []\nraw.info['bads'] += ['MEG 2443', 'EEG 053'] # bads + 2 more\n\n# picks MEG gradiometers\npicks = mne.pick_types(raw.info, meg='grad', eeg=False, eog=True,\n stim=False, include=include, exclude='bads')\n\nch_name = 'MEG 1332' # restrict example to one channel\n\n# Load condition 1\nreject = dict(grad=4000e-13, eog=150e-6)\nevent_id = 1\nepochs_condition_1 = mne.Epochs(raw, events, event_id, tmin, tmax,\n picks=picks, baseline=(None, 0),\n reject=reject, preload=True)\nepochs_condition_1.pick_channels([ch_name])\n\n# Load condition 2\nevent_id = 2\nepochs_condition_2 = mne.Epochs(raw, events, event_id, tmin, tmax,\n picks=picks, baseline=(None, 0),\n reject=reject, preload=True)\nepochs_condition_2.pick_channels([ch_name])" ], "outputs": [], "metadata": { "collapsed": false } }, { "source": [ "Factor to downsample the temporal dimension of the TFR computed by\ntfr_morlet. Decimation occurs after frequency decomposition and can\nbe used to reduce memory usage (and possibly comptuational time of downstream\noperations such as nonparametric statistics) if you don't need high\nspectrotemporal resolution.\n\n" ], "cell_type": "markdown", "metadata": {} }, { "execution_count": null, "cell_type": "code", "source": [ "decim = 2\nfrequencies = np.arange(7, 30, 3) # define frequencies of interest\nn_cycles = 1.5\n\ntfr_epochs_1 = tfr_morlet(epochs_condition_1, frequencies,\n n_cycles=n_cycles, decim=decim,\n return_itc=False, average=False)\n\ntfr_epochs_2 = tfr_morlet(epochs_condition_2, frequencies,\n n_cycles=n_cycles, decim=decim,\n return_itc=False, average=False)\n\ntfr_epochs_1.apply_baseline(mode='ratio', baseline=(None, 0))\ntfr_epochs_2.apply_baseline(mode='ratio', baseline=(None, 0))\n\nepochs_power_1 = tfr_epochs_1.data[:, 0, :, :] # only 1 channel as 3D matrix\nepochs_power_2 = tfr_epochs_2.data[:, 0, :, :] # only 1 channel as 3D matrix" ], "outputs": [], "metadata": { "collapsed": false } }, { "source": [ "Compute statistic\n-----------------\n\n" ], "cell_type": "markdown", "metadata": {} }, { "execution_count": null, "cell_type": "code", "source": [ "threshold = 6.0\nT_obs, clusters, cluster_p_values, H0 = \\\n permutation_cluster_test([epochs_power_1, epochs_power_2],\n n_permutations=100, threshold=threshold, tail=0)" ], "outputs": [], "metadata": { "collapsed": false } }, { "source": [ "View time-frequency plots\n-------------------------\n\n" ], "cell_type": "markdown", "metadata": {} }, { "execution_count": null, "cell_type": "code", "source": [ "times = 1e3 * epochs_condition_1.times # change unit to ms\nevoked_condition_1 = epochs_condition_1.average()\nevoked_condition_2 = epochs_condition_2.average()\n\nplt.figure()\nplt.subplots_adjust(0.12, 0.08, 0.96, 0.94, 0.2, 0.43)\n\nplt.subplot(2, 1, 1)\n# Create new stats image with only significant clusters\nT_obs_plot = np.nan * np.ones_like(T_obs)\nfor c, p_val in zip(clusters, cluster_p_values):\n if p_val <= 0.05:\n T_obs_plot[c] = T_obs[c]\n\nplt.imshow(T_obs,\n extent=[times[0], times[-1], frequencies[0], frequencies[-1]],\n aspect='auto', origin='lower', cmap='gray')\nplt.imshow(T_obs_plot,\n extent=[times[0], times[-1], frequencies[0], frequencies[-1]],\n aspect='auto', origin='lower', cmap='RdBu_r')\n\nplt.xlabel('Time (ms)')\nplt.ylabel('Frequency (Hz)')\nplt.title('Induced power (%s)' % ch_name)\n\nax2 = plt.subplot(2, 1, 2)\nevoked_contrast = mne.combine_evoked([evoked_condition_1, evoked_condition_2],\n weights=[1, -1])\nevoked_contrast.plot(axes=ax2)\n\nplt.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" } } } }