{
  "nbformat_minor": 0, 
  "nbformat": 4, 
  "cells": [
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "%matplotlib inline"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "\n# Permutation F-test on sensor data with 1D cluster level\n\n\nOne tests if the evoked response is significantly different\nbetween conditions. Multiple comparison problem is addressed\nwith cluster level permutation test.\n\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "# Authors: Alexandre Gramfort <alexandre.gramfort@telecom-paristech.fr>\n#\n# License: BSD (3-clause)\n\nimport matplotlib.pyplot as plt\n\nimport mne\nfrom mne import io\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_filt-0-40_raw.fif'\nevent_fname = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw-eve.fif'\nevent_id = 1\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]"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "Read epochs for the channel of interest\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "picks = mne.pick_types(raw.info, meg=False, eog=True, include=include,\n                       exclude='bads')\nevent_id = 1\nreject = dict(grad=4000e-13, eog=150e-6)\nepochs1 = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks,\n                     baseline=(None, 0), reject=reject)\ncondition1 = epochs1.get_data()  # as 3D matrix\n\nevent_id = 2\nepochs2 = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks,\n                     baseline=(None, 0), reject=reject)\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"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "Compute statistic\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([condition1, condition2], n_permutations=1000,\n                             threshold=threshold, tail=1, n_jobs=1)"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "Plot\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "times = epochs1.times\nplt.close('all')\nplt.subplot(211)\nplt.title('Channel : ' + channel)\nplt.plot(times, condition1.mean(axis=0) - condition2.mean(axis=0),\n         label=\"ERF Contrast (Event 1 - Event 2)\")\nplt.ylabel(\"MEG (T / m)\")\nplt.legend()\nplt.subplot(212)\nfor i_c, c in enumerate(clusters):\n    c = c[0]\n    if cluster_p_values[i_c] <= 0.05:\n        h = plt.axvspan(times[c.start], times[c.stop - 1],\n                        color='r', alpha=0.3)\n    else:\n        plt.axvspan(times[c.start], times[c.stop - 1], color=(0.3, 0.3, 0.3),\n                    alpha=0.3)\nhf = plt.plot(times, T_obs, 'g')\nplt.legend((h, ), ('cluster p-value < 0.05', ))\nplt.xlabel(\"time (ms)\")\nplt.ylabel(\"f-values\")\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"
      }
    }
  }
}