{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n\n# Mass-univariate twoway repeated measures ANOVA on single trial power\n\nThis script shows how to conduct a mass-univariate repeated measures\nANOVA. As the model to be fitted assumes two fully crossed factors,\nwe will study the interplay between perceptual modality\n(auditory VS visual) and the location of stimulus presentation\n(left VS right). Here we use single trials as replications\n(subjects) while iterating over time slices plus frequency bands\nfor to fit our mass-univariate model. For the sake of simplicity we\nwill confine this analysis to one single channel of which we know\nthat it exposes a strong induced response. We will then visualize\neach effect by creating a corresponding mass-univariate effect\nimage. We conclude with accounting for multiple comparisons by\nperforming a permutation clustering test using the ANOVA as\nclustering function. The results final will be compared to multiple\ncomparisons using False Discovery Rate correction.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Authors: Denis Engemann \n# Eric Larson \n# Alexandre Gramfort \n# Alex Rockhill \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 f_mway_rm, f_threshold_mway_rm, fdr_correction\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\"] # bads\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\"\n\n# Load conditions\nreject = dict(grad=4000e-13, eog=150e-6)\nevent_id = dict(aud_l=1, aud_r=2, vis_l=3, vis_r=4)\nepochs = mne.Epochs(\n raw,\n events,\n event_id,\n tmin,\n tmax,\n picks=picks,\n baseline=(None, 0),\n preload=True,\n reject=reject,\n)\nepochs.pick([ch_name]) # restrict example to one channel" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We have to make sure all conditions have the same counts, as the ANOVA\nexpects a fully balanced data matrix and does not forgive imbalances that\ngenerously (risk of type-I error).\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "epochs.equalize_event_counts(event_id)\n\n# Factor to down-sample the temporal dimension of the TFR computed by\n# tfr_morlet.\ndecim = 2\nfreqs = np.arange(7, 30, 3) # define frequencies of interest\nn_cycles = freqs / freqs[0]\nzero_mean = False # don't correct morlet wavelet to be of mean zero\n# To have a true wavelet zero_mean should be True but here for illustration\n# purposes it helps to spot the evoked response." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Create TFR representations for all conditions\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "epochs_power = list()\nfor condition in [epochs[k] for k in event_id]:\n this_tfr = condition.compute_tfr(\n \"morlet\",\n freqs,\n n_cycles=n_cycles,\n decim=decim,\n average=False,\n zero_mean=zero_mean,\n return_itc=False,\n )\n this_tfr.apply_baseline(mode=\"ratio\", baseline=(None, 0))\n this_power = this_tfr.data[:, 0, :, :] # we only have one channel.\n epochs_power.append(this_power)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Setup repeated measures ANOVA\n\nWe will tell the ANOVA how to interpret the data matrix in terms of factors.\nThis is done via the factor levels argument which is a list of the number\nfactor levels for each factor.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "n_conditions = len(epochs.event_id)\nn_replications = epochs.events.shape[0] // n_conditions\n\nfactor_levels = [2, 2] # number of levels in each factor\neffects = \"A*B\" # this is the default signature for computing all effects\n# Other possible options are 'A' or 'B' for the corresponding main effects\n# or 'A:B' for the interaction effect only (this notation is borrowed from the\n# R formula language)\nn_freqs = len(freqs)\ntimes = 1e3 * epochs.times[::decim]\nn_times = len(times)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we'll assemble the data matrix and swap axes so the trial replications\nare the first dimension and the conditions are the second dimension.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "data = np.swapaxes(np.asarray(epochs_power), 1, 0)\n\n# so we have replications \u00d7 conditions \u00d7 observations\n# where the time-frequency observations are freqs \u00d7 times:\nprint(data.shape)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "While the iteration scheme used above for assembling the data matrix\nmakes sure the first two dimensions are organized as expected (with A =\nmodality and B = location):\n\n.. table:: Sample data layout\n\n ===== ==== ==== ==== ====\n trial A1B1 A1B2 A2B1 B2B2\n ===== ==== ==== ==== ====\n 1 1.34 2.53 0.97 1.74\n ... ... ... ... ...\n 56 2.45 7.90 3.09 4.76\n ===== ==== ==== ==== ====\n\nNow we're ready to run our repeated measures ANOVA.\n\nNote. As we treat trials as subjects, the test only accounts for\ntime locked responses despite the 'induced' approach.\nFor analysis for induced power at the group level averaged TRFs\nare required.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fvals, pvals = f_mway_rm(data, factor_levels, effects=effects)\n\neffect_labels = [\"modality\", \"location\", \"modality by location\"]\n\nfig, axes = plt.subplots(3, 1, figsize=(6, 6), layout=\"constrained\")\n\n# let's visualize our effects by computing f-images\nfor effect, sig, effect_label, ax in zip(fvals, pvals, effect_labels, axes):\n # show naive F-values in gray\n ax.imshow(\n effect,\n cmap=\"gray\",\n aspect=\"auto\",\n origin=\"lower\",\n extent=[times[0], times[-1], freqs[0], freqs[-1]],\n )\n # create mask for significant time-frequency locations\n effect[sig >= 0.05] = np.nan\n c = ax.imshow(\n effect,\n cmap=\"autumn\",\n aspect=\"auto\",\n origin=\"lower\",\n extent=[times[0], times[-1], freqs[0], freqs[-1]],\n )\n fig.colorbar(c, ax=ax)\n ax.set_xlabel(\"Time (ms)\")\n ax.set_ylabel(\"Frequency (Hz)\")\n ax.set_title(f'Time-locked response for \"{effect_label}\" ({ch_name})')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Account for multiple comparisons using FDR versus permutation clustering test\n\nFirst we need to slightly modify the ANOVA function to be suitable for\nthe clustering procedure. Also want to set some defaults.\nLet's first override effects to confine the analysis to the interaction\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "effects = \"A:B\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A stat_fun must deal with a variable number of input arguments.\nInside the clustering function each condition will be passed as flattened\narray, necessitated by the clustering procedure. The ANOVA however expects an\ninput array of dimensions: subjects \u00d7 conditions \u00d7 observations (optional).\nThe following function catches the list input and swaps the first and\nthe second dimension and finally calls the ANOVA function.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def stat_fun(*args):\n return f_mway_rm(\n np.swapaxes(args, 1, 0),\n factor_levels=factor_levels,\n effects=effects,\n return_pvals=False,\n )[0]\n\n\n# The ANOVA returns a tuple f-values and p-values, we will pick the former.\npthresh = 0.001 # set threshold rather high to save some time\nf_thresh = f_threshold_mway_rm(n_replications, factor_levels, effects, pthresh)\ntail = 1 # f-test, so tail > 0\nn_permutations = 256 # Save some time (the test won't be too sensitive ...)\nF_obs, clusters, cluster_p_values, h0 = mne.stats.permutation_cluster_test(\n epochs_power,\n stat_fun=stat_fun,\n threshold=f_thresh,\n tail=tail,\n n_jobs=None,\n n_permutations=n_permutations,\n buffer_size=None,\n out_type=\"mask\",\n seed=0,\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create new stats image with only significant clusters:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "good_clusters = np.where(cluster_p_values < 0.05)[0]\nF_obs_plot = np.full_like(F_obs, np.nan)\nfor ii in good_clusters:\n F_obs_plot[clusters[ii]] = F_obs[clusters[ii]]\n\nfig, ax = plt.subplots(figsize=(6, 4), layout=\"constrained\")\nfor f_image, cmap in zip([F_obs, F_obs_plot], [\"gray\", \"autumn\"]):\n c = ax.imshow(\n f_image,\n cmap=cmap,\n aspect=\"auto\",\n origin=\"lower\",\n extent=[times[0], times[-1], freqs[0], freqs[-1]],\n )\n\nfig.colorbar(c, ax=ax)\nax.set_xlabel(\"Time (ms)\")\nax.set_ylabel(\"Frequency (Hz)\")\nax.set_title(\n f'Time-locked response for \"modality by location\" ({ch_name})\\n'\n \"cluster-level corrected (p <= 0.05)\"\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now using FDR:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "mask, _ = fdr_correction(pvals[2])\nF_obs_plot2 = F_obs.copy()\nF_obs_plot2[~mask.reshape(F_obs_plot.shape)] = np.nan\n\nfig, ax = plt.subplots(figsize=(6, 4), layout=\"constrained\")\nfor f_image, cmap in zip([F_obs, F_obs_plot2], [\"gray\", \"autumn\"]):\n c = ax.imshow(\n f_image,\n cmap=cmap,\n aspect=\"auto\",\n origin=\"lower\",\n extent=[times[0], times[-1], freqs[0], freqs[-1]],\n )\n\nfig.colorbar(c, ax=ax)\nax.set_xlabel(\"Time (ms)\")\nax.set_ylabel(\"Frequency (Hz)\")\nax.set_title(\n f'Time-locked response for \"modality by location\" ({ch_name})\\n'\n \"FDR corrected (p <= 0.05)\"\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Both cluster-level and FDR correction help get rid of potential\nfalse-positives that we saw in the naive f-images. The cluster permutation\ncorrection is biased toward time-frequencies with contiguous areas of high\nor low power, which is likely appropriate given the highly correlated nature\nof this data. This is the most likely explanation for why one cluster was\npreserved by the cluster permutation correction, but no time-frequencies\nwere significant using the FDR correction.\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 }