{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n\n# FDR correction on T-test on sensor data\n\nOne tests if the evoked response significantly deviates from 0.\nMultiple comparison problem is addressed with\nFalse Discovery Rate (FDR) correction.\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\nimport numpy as np\nfrom scipy import stats\n\nimport mne\nfrom mne import io\nfrom mne.datasets import sample\nfrom mne.stats import bonferroni_correction, 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_filt-0-40_raw.fif\"\nevent_fname = meg_path / \"sample_audvis_filt-0-40_raw-eve.fif\"\nevent_id, tmin, tmax = 1, -0.2, 0.5\n\n#   Setup for reading the raw data\nraw = io.read_raw_fif(raw_fname)\nevents = mne.read_events(event_fname)[:30]\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)\nepochs = mne.Epochs(\n    raw, events, event_id, tmin, tmax, picks=picks, baseline=(None, 0), reject=reject\n)\nX = epochs.get_data()  # as 3D matrix\nX = X[:, 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": [
        "T, pval = stats.ttest_1samp(X, 0)\nalpha = 0.05\n\nn_samples, n_tests = X.shape\nthreshold_uncorrected = stats.t.ppf(1.0 - alpha, n_samples - 1)\n\nreject_bonferroni, pval_bonferroni = bonferroni_correction(pval, alpha=alpha)\nthreshold_bonferroni = stats.t.ppf(1.0 - alpha / n_tests, n_samples - 1)\n\nreject_fdr, pval_fdr = fdr_correction(pval, alpha=alpha, method=\"indep\")\nthreshold_fdr = np.min(np.abs(T)[reject_fdr])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plot\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "times = 1e3 * epochs.times\n\nplt.close(\"all\")\nplt.plot(times, T, \"k\", label=\"T-stat\")\nxmin, xmax = plt.xlim()\nplt.hlines(\n    threshold_uncorrected,\n    xmin,\n    xmax,\n    linestyle=\"--\",\n    colors=\"k\",\n    label=\"p=0.05 (uncorrected)\",\n    linewidth=2,\n)\nplt.hlines(\n    threshold_bonferroni,\n    xmin,\n    xmax,\n    linestyle=\"--\",\n    colors=\"r\",\n    label=\"p=0.05 (Bonferroni)\",\n    linewidth=2,\n)\nplt.hlines(\n    threshold_fdr,\n    xmin,\n    xmax,\n    linestyle=\"--\",\n    colors=\"b\",\n    label=\"p=0.05 (FDR)\",\n    linewidth=2,\n)\nplt.legend()\nplt.xlabel(\"Time (ms)\")\nplt.ylabel(\"T-stat\")\nplt.show()"
      ]
    }
  ],
  "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
}