{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n\n# Permutation F-test on sensor data with 1D cluster level\n\nOne tests if the evoked response is significantly different\nbetween conditions. Multiple comparison problem is addressed\nwith cluster level permutation test.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Authors: Alexandre Gramfort <alexandre.gramfort@inria.fr>\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\n\nimport mne\nfrom mne import io\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_filt-0-40_raw.fif\"\nevent_fname = meg_path / \"sample_audvis_filt-0-40_raw-eve.fif\"\ntmin = -0.2\ntmax = 0.5\n\n#   Setup for reading the raw data\nraw = io.read_raw_fif(raw_fname)\nevents = mne.read_events(event_fname)\n\nchannel = \"MEG 1332\"  # include only this channel in analysis\ninclude = [channel]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Read epochs for the channel of interest\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "picks = mne.pick_types(raw.info, meg=False, eog=True, include=include, exclude=\"bads\")\nevent_id = 1\nreject = dict(grad=4000e-13, eog=150e-6)\nepochs1 = mne.Epochs(\n    raw, events, event_id, tmin, tmax, picks=picks, baseline=(None, 0), reject=reject\n)\ncondition1 = epochs1.get_data()  # as 3D matrix\n\nevent_id = 2\nepochs2 = mne.Epochs(\n    raw, events, event_id, tmin, tmax, picks=picks, baseline=(None, 0), reject=reject\n)\ncondition2 = epochs2.get_data()  # as 3D matrix\n\ncondition1 = condition1[:, 0, :]  # take only one channel to get a 2D array\ncondition2 = condition2[:, 0, :]  # take only one channel to get a 2D array"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Compute statistic\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "threshold = 6.0\nT_obs, clusters, cluster_p_values, H0 = permutation_cluster_test(\n    [condition1, condition2],\n    n_permutations=1000,\n    threshold=threshold,\n    tail=1,\n    n_jobs=None,\n    out_type=\"mask\",\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plot\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "times = epochs1.times\nfig, (ax, ax2) = plt.subplots(2, 1, figsize=(8, 4))\nax.set_title(\"Channel : \" + channel)\nax.plot(\n    times,\n    condition1.mean(axis=0) - condition2.mean(axis=0),\n    label=\"ERF Contrast (Event 1 - Event 2)\",\n)\nax.set_ylabel(\"MEG (T / m)\")\nax.legend()\n\nfor i_c, c in enumerate(clusters):\n    c = c[0]\n    if cluster_p_values[i_c] <= 0.05:\n        h = ax2.axvspan(times[c.start], times[c.stop - 1], color=\"r\", alpha=0.3)\n    else:\n        ax2.axvspan(times[c.start], times[c.stop - 1], color=(0.3, 0.3, 0.3), alpha=0.3)\n\nhf = plt.plot(times, T_obs, \"g\")\nax2.legend((h,), (\"cluster p-value < 0.05\",))\nax2.set_xlabel(\"time (ms)\")\nax2.set_ylabel(\"f-values\")"
      ]
    }
  ],
  "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
}