{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n\n# Non-parametric between conditions cluster statistic on single trial power\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 of:\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" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Authors: Alexandre Gramfort \n#\n# License: BSD-3-Clause\n# Copyright the MNE-Python contributors." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import matplotlib.pyplot as plt\nimport numpy as np\n\nimport mne\nfrom mne.datasets import sample\nfrom mne.stats import permutation_cluster_test\n\nprint(__doc__)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Set parameters\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "data_path = sample.data_path()\nmeg_path = data_path / \"MEG\" / \"sample\"\nraw_fname = meg_path / \"sample_audvis_raw.fif\"\nevent_fname = meg_path / \"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(\n raw.info,\n meg=\"grad\",\n eeg=False,\n eog=True,\n stim=False,\n include=include,\n exclude=\"bads\",\n)\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(\n raw,\n events,\n event_id,\n tmin,\n tmax,\n picks=picks,\n baseline=(None, 0),\n reject=reject,\n preload=True,\n)\nepochs_condition_1.pick([ch_name])\n\n# Load condition 2\nevent_id = 2\nepochs_condition_2 = mne.Epochs(\n raw,\n events,\n event_id,\n tmin,\n tmax,\n picks=picks,\n baseline=(None, 0),\n reject=reject,\n preload=True,\n)\nepochs_condition_2.pick([ch_name])" ] }, { "cell_type": "markdown", "metadata": {}, "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": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "decim = 2\nfreqs = np.arange(7, 30, 3) # define frequencies of interest\nn_cycles = 1.5\ntfr_kwargs = dict(\n method=\"morlet\",\n freqs=freqs,\n n_cycles=n_cycles,\n decim=decim,\n return_itc=False,\n average=False,\n)\n\ntfr_epochs_1 = epochs_condition_1.compute_tfr(**tfr_kwargs)\ntfr_epochs_2 = epochs_condition_2.compute_tfr(**tfr_kwargs)\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" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Compute statistic\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "threshold = 6.0\nF_obs, clusters, cluster_p_values, H0 = permutation_cluster_test(\n [epochs_power_1, epochs_power_2],\n out_type=\"mask\",\n n_permutations=100,\n threshold=threshold,\n tail=0,\n seed=np.random.default_rng(seed=8675309),\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## View time-frequency plots\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "times = 1e3 * epochs_condition_1.times # change unit to ms\n\nfig, (ax, ax2) = plt.subplots(2, 1, figsize=(6, 4), layout=\"constrained\")\n\n# Compute the difference in evoked to determine which was greater since\n# we used a 1-way ANOVA which tested for a difference in population means\nevoked_power_1 = epochs_power_1.mean(axis=0)\nevoked_power_2 = epochs_power_2.mean(axis=0)\nevoked_power_contrast = evoked_power_1 - evoked_power_2\nsigns = np.sign(evoked_power_contrast)\n\n# Create new stats image with only significant clusters\nF_obs_plot = np.nan * np.ones_like(F_obs)\nfor c, p_val in zip(clusters, cluster_p_values):\n if p_val <= 0.05:\n F_obs_plot[c] = F_obs[c] * signs[c]\n\nax.imshow(\n F_obs,\n extent=[times[0], times[-1], freqs[0], freqs[-1]],\n aspect=\"auto\",\n origin=\"lower\",\n cmap=\"gray\",\n)\nmax_F = np.nanmax(abs(F_obs_plot))\nax.imshow(\n F_obs_plot,\n extent=[times[0], times[-1], freqs[0], freqs[-1]],\n aspect=\"auto\",\n origin=\"lower\",\n cmap=\"RdBu_r\",\n vmin=-max_F,\n vmax=max_F,\n)\n\nax.set_xlabel(\"Time (ms)\")\nax.set_ylabel(\"Frequency (Hz)\")\nax.set_title(f\"Induced power ({ch_name})\")\n\n# plot evoked\nevoked_condition_1 = epochs_condition_1.average()\nevoked_condition_2 = epochs_condition_2.average()\nevoked_contrast = mne.combine_evoked(\n [evoked_condition_1, evoked_condition_2], weights=[1, -1]\n)\nevoked_contrast.plot(axes=ax2, time_unit=\"s\")" ] } ], "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 }