{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n\n# Auto-generating Epochs metadata\n\nThis tutorial shows how to auto-generate metadata for `~mne.Epochs`, based on\nevents via `mne.epochs.make_metadata`.\n\nWe are going to use data from the `erp-core-dataset` (derived from\n:footcite:`Kappenman2021`). This is EEG data from a single participant\nperforming an active visual task (Eriksen flanker task).\n\n<div class=\"alert alert-info\"><h4>Note</h4><p>If you wish to skip the introductory parts of this tutorial, you may jump\n   straight to `tut-autogenerate-metadata-ern` after completing the data\n   import and event creation in the\n   `tut-autogenerate-metadata-preparation` section.</p></div>\n\nThis tutorial is loosely divided into two parts:\n\n1. We will first focus on producing ERPs time-locked to the **visual\n   stimulation**, conditional on response correctness and response time in\n   order to familiarize ourselves with the `~mne.epochs.make_metadata`\n   function.\n2. After that, we will calculate ERPs time-locked to the **responses** \u2013 again,\n   conditional on response correctness \u2013 to visualize the error-related\n   negativity (ERN), i.e. the ERP component associated with incorrect\n   behavioral responses.\n\n\n\n## Preparation\n\nLet's start by reading, filtering, and producing a simple visualization of the\nraw data. The data is pretty clean and contains very few blinks, so there's no\nneed to apply sophisticated preprocessing and data cleaning procedures.\nWe will also convert the `~mne.Annotations` contained in this dataset to events\nby calling `mne.events_from_annotations`.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Authors: The MNE-Python contributors.\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\n\ndata_dir = mne.datasets.erp_core.data_path()\ninfile = data_dir / \"ERP-CORE_Subject-001_Task-Flankers_eeg.fif\"\n\nraw = mne.io.read_raw(infile, preload=True)\nraw.filter(l_freq=0.1, h_freq=40)\nraw.plot(start=60)\n\n# extract events\nall_events, all_event_id = mne.events_from_annotations(raw)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Creating metadata from events\n\n### The basics of ``make_metadata``\n\nNow it's time to think about the time windows to use for epoching and\nmetadata generation. **It is important to understand that these time windows\nneed not be the same!** That is, the automatically generated metadata might\ninclude information about events from only a fraction of the epochs duration;\nor it might include events that occurred well outside a given epoch.\n\nLet us look at a concrete example. In the Flankers task of the ERP CORE\ndataset, participants were required to respond to visual stimuli by pressing\na button. We're interested in looking at the visual evoked responses (ERPs)\nof trials with correct responses. Assume that based on literature\nstudies, we decide that responses later than 1500 ms after stimulus onset are\nto be considered invalid, because they don't capture the neuronal processes\nof interest here. We can approach this  in the following way with the help of\n`mne.epochs.make_metadata`:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# metadata for each epoch shall include events from the range: [0.0, 1.5] s,\n# i.e. starting with stimulus onset and expanding beyond the end of the epoch\nmetadata_tmin, metadata_tmax = 0.0, 1.5\n\n# auto-create metadata:\n# this also returns a new events array and an event_id dictionary. we'll see\n# later why this is important\nmetadata, events, event_id = mne.epochs.make_metadata(\n    events=all_events,\n    event_id=all_event_id,\n    tmin=metadata_tmin,\n    tmax=metadata_tmax,\n    sfreq=raw.info[\"sfreq\"],\n)\n\n# let's look at what we got!\nmetadata"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Specifying time-locked events\n\nWe can see that the generated table has 802 rows, each one corresponding to\nan individual event in ``all_events``. The first column, ``event_name``,\ncontains the name of the respective event around which the metadata of that\nspecific column was generated \u2013 we'll call that the \"time-locked event\",\nbecause we'll assign it time point zero.\n\nThe names of the remaining columns correspond to the event names specified in\nthe ``all_event_id`` dictionary. These columns contain floats; the values\nrepresent the latency of that specific event in seconds, relative to\nthe time-locked event (the one mentioned in the ``event_name`` column).\nFor events that didn't occur within the given time window, you'll see\na value of ``NaN``, simply indicating that no event latency could be\nextracted.\n\nNow, there's a problem here. We want to investigate the visual ERPs\nconditional on responses. But the metadata that was just created contains\none row for **every** event, including responses. While we **could** create\nepochs for all events, allowing us to pass those metadata, and later subset\nthe created events, there's a more elegant way to handle this:\n`~mne.epochs.make_metadata` has a ``row_events`` parameter that\nallows us to specify for which events to create metadata **rows**, while\nstill creating **columns for all events** in the ``event_id`` dictionary.\n\nBecause the metadata, then, only pertains to a subset of our original events,\nit's important to keep the returned ``events`` and ``event_id`` around for\nlater use when we actually create our epochs, to ensure that metadata,\nevents, and event descriptions stay in sync.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "row_events = [\n    \"stimulus/compatible/target_left\",\n    \"stimulus/compatible/target_right\",\n    \"stimulus/incompatible/target_left\",\n    \"stimulus/incompatible/target_right\",\n]\n\nmetadata, events, event_id = mne.epochs.make_metadata(\n    events=all_events,\n    event_id=all_event_id,\n    tmin=metadata_tmin,\n    tmax=metadata_tmax,\n    sfreq=raw.info[\"sfreq\"],\n    row_events=row_events,\n)\n\nmetadata"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Keeping only the first events of a group\n\nThe metadata now contains 400 rows \u2013 one per stimulation \u2013 and the same\nnumber of columns as before. Great!\n\nWe have two types of responses in our data: ``response/left`` and\n``response/right``. We would like to map those to \"correct\" and \"incorrect\".\nTo make this easier, we can ask `~mne.epochs.make_metadata` to generate an\nentirely **new** column that refers to the first response observed during the\ngiven time interval. This works by passing a subset of the\n:term:`hierarchical event descriptors` (HEDs, inspired by\n:footcite:`BigdelyShamloEtAl2013`) used to name events via the ``keep_first``\nparameter. For example, in the case of the HEDs ``response/left`` and\n``response/right``, we could pass ``keep_first='response'`` to generate a new\ncolumn, ``response``, containing the latency of the respective event. This\nvalue represents the first (or, in this specific example: the only)\nresponse, regardless of side (left or right). To indicate **which** event\ntype (here: response side) was matched, a second column named\n``first_response`` is added. The values in this column are the event types\nwithout the string used for matching, as it is already encoded as the column\nname, i.e. in our example, we expect it to only contain ``'left'`` and\n``'right'``.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "keep_first = \"response\"\nmetadata, events, event_id = mne.epochs.make_metadata(\n    events=all_events,\n    event_id=all_event_id,\n    tmin=metadata_tmin,\n    tmax=metadata_tmax,\n    sfreq=raw.info[\"sfreq\"],\n    row_events=row_events,\n    keep_first=keep_first,\n)\n\n# visualize response times regardless of side\nmetadata[\"response\"].plot.hist(bins=50, title=\"Response Times\")\n\n# the \"first_response\" column contains only \"left\" and \"right\" entries, derived\n# from the initial event named \"response/left\" and \"response/right\"\nprint(metadata[\"first_response\"])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We're facing a similar issue with the stimulus events, and now there are not\nonly two, but **four** different types: ``stimulus/compatible/target_left``,\n``stimulus/compatible/target_right``, ``stimulus/incompatible/target_left``,\nand ``stimulus/incompatible/target_right``. What's more, because in the\npresent paradigm stimuli were presented in rapid succession, sometimes\nmultiple stimulus events occurred within the 1.5 second time window we used\nto generate our metadata. See for example:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "metadata.loc[\n    metadata[\"stimulus/compatible/target_left\"].notna()\n    & metadata[\"stimulus/compatible/target_right\"].notna(),\n    :,\n]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This can easily lead to confusion during later stages of processing, so let's\ncreate a column for the first stimulus \u2013 which will always be the time-locked\nstimulus, as our time interval starts at 0 seconds. We can pass a **list** of\nstrings to ``keep_first``.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "keep_first = [\"stimulus\", \"response\"]\nmetadata, events, event_id = mne.epochs.make_metadata(\n    events=all_events,\n    event_id=all_event_id,\n    tmin=metadata_tmin,\n    tmax=metadata_tmax,\n    sfreq=raw.info[\"sfreq\"],\n    row_events=row_events,\n    keep_first=keep_first,\n)\n\n# all times of the time-locked events should be zero\nassert all(metadata[\"stimulus\"] == 0)\n\n# the values in the new \"first_stimulus\" and \"first_response\" columns indicate\n# which events were selected via \"keep_first\"\nmetadata[[\"first_stimulus\", \"first_response\"]]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Adding new columns to describe stimulation side and response correctness\n\nPerfect! Now it's time to define which responses were correct and incorrect.\nWe first add a column encoding the side of stimulation, and then simply\ncheck whether the response matches the stimulation side, and add this result\nto another column.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "metadata.loc[:, \"stimulus_side\"] = \"\"  # initialize column\n\n# left-side stimulation\nmetadata.loc[\n    metadata[\"first_stimulus\"].isin(\n        [\"compatible/target_left\", \"incompatible/target_left\"]\n    ),\n    \"stimulus_side\",\n] = \"left\"\n\n# right-side stimulation\nmetadata.loc[\n    metadata[\"first_stimulus\"].isin(\n        [\"compatible/target_right\", \"incompatible/target_right\"]\n    ),\n    \"stimulus_side\",\n] = \"right\"\n\n# first assume all responses were incorrect, then mark those as correct where\n# the stimulation side matches the response side\nmetadata[\"response_correct\"] = False\nmetadata.loc[\n    metadata[\"stimulus_side\"] == metadata[\"first_response\"], \"response_correct\"\n] = True\n\n\ncorrect_response_count = metadata[\"response_correct\"].sum()\nprint(\n    f\"Correct responses: {correct_response_count}\\n\"\n    f\"Incorrect responses: {len(metadata) - correct_response_count}\"\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Creating ``Epochs`` with metadata, and visualizing ERPs\n\nIt's finally time to create our epochs! We set the metadata directly on\ninstantiation via the ``metadata`` parameter. Also, it is important to\nremember to pass ``events`` and ``event_id`` as returned from\n`~mne.epochs.make_metadata`, as we only created metadata for a subset of\nour original events by passing ``row_events``. Otherwise, the length\nof the metadata and the number of epochs would not match, which would raise\nan error.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "epochs_tmin, epochs_tmax = -0.1, 0.4  # epochs range: [-0.1, 0.4] s\nreject = {\"eeg\": 250e-6}  # exclude epochs with strong artifacts\nepochs = mne.Epochs(\n    raw=raw,\n    tmin=epochs_tmin,\n    tmax=epochs_tmax,\n    events=events,\n    event_id=event_id,\n    metadata=metadata,\n    reject=reject,\n    preload=True,\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Lastly, let's visualize the ERPs associated with the visual stimulation, once\nfor all trials with correct responses, and once for all trials with correct\nresponses and a response time greater than 0.5 seconds\n(i.e., slow responses).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "vis_erp = epochs[\"response_correct\"].average()\nvis_erp_slow = epochs[\"(not response_correct) & (response > 0.3)\"].average()\n\nfig, ax = plt.subplots(2, figsize=(6, 6), layout=\"constrained\")\nvis_erp.plot(gfp=True, spatial_colors=True, axes=ax[0])\nvis_erp_slow.plot(gfp=True, spatial_colors=True, axes=ax[1])\nax[0].set_title(\"Visual ERPs \u2013 All Correct Responses\")\nax[1].set_title(\"Visual ERPs \u2013 Slow Correct Responses\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Aside from the fact that the data for the (much fewer) slow responses looks\nnoisier \u2013\u00a0which is entirely to be expected \u2013\u00a0not much of an ERP difference\ncan be seen.\n\n\n## Applying the knowledge: visualizing the ERN component\n\nIn the following analysis, we will use the same dataset as above, but\nwe'll time-lock our epochs to the **response events,** not to the stimulus\nonset. Comparing ERPs associated with correct and incorrect behavioral\nresponses, we should be able to see the error-related negativity (ERN) in\nthe difference wave.\n\nSince we want to time-lock our analysis to responses, for the automated\nmetadata generation we'll consider events occurring up to 1500 ms before\nthe response trigger.\n\nWe only wish to consider the **last** stimulus and response in each time\nwindow: Remember that we're dealing with rapid stimulus presentations in\nthis paradigm; taking the last response (at time point zero) and the last\nstimulus (the one closest to the response) ensures that we actually create\nthe right stimulus-response pairings. We can achieve this by passing the\n``keep_last`` parameter, which works exactly like ``keep_first`` we used\npreviously, only that it keeps the **last** occurrences of the specified\nevents and stores them in columns whose names start with ``last_``.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "metadata_tmin, metadata_tmax = -1.5, 0\nrow_events = [\"response/left\", \"response/right\"]\nkeep_last = [\"stimulus\", \"response\"]\n\nmetadata, events, event_id = mne.epochs.make_metadata(\n    events=all_events,\n    event_id=all_event_id,\n    tmin=metadata_tmin,\n    tmax=metadata_tmax,\n    sfreq=raw.info[\"sfreq\"],\n    row_events=row_events,\n    keep_last=keep_last,\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Exactly like in the previous example, we create new columns ``stimulus_side``\nand ``response_correct``.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "metadata.loc[:, \"stimulus_side\"] = \"\"  # initialize column\n\n# left-side stimulation\nmetadata.loc[\n    metadata[\"last_stimulus\"].isin(\n        [\"compatible/target_left\", \"incompatible/target_left\"]\n    ),\n    \"stimulus_side\",\n] = \"left\"\n\n# right-side stimulation\nmetadata.loc[\n    metadata[\"last_stimulus\"].isin(\n        [\"compatible/target_right\", \"incompatible/target_right\"]\n    ),\n    \"stimulus_side\",\n] = \"right\"\n\n# first assume all responses were incorrect, then mark those as correct where\n# the stimulation side matches the response side\nmetadata[\"response_correct\"] = False\nmetadata.loc[\n    metadata[\"stimulus_side\"] == metadata[\"last_response\"], \"response_correct\"\n] = True\n\nmetadata"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now it's already time to epoch the data! When deciding upon the epochs\nduration for this specific analysis, we need to ensure to include quite a bit\nof signal from before and after the motor response. We also must be aware of\nthe fact that motor-/muscle-related signals will most likely be present\n**before** the response button trigger pulse appears in our data, so the time\nperiod close to the response event should not be used for baseline\ncorrection. But at the same time, we don't want to use a baseline\nperiod that extends too far away from the button event. The following values\nseem to work quite well.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "epochs_tmin, epochs_tmax = -0.6, 0.4\nbaseline = (-0.4, -0.2)\nreject = {\"eeg\": 250e-6}\nepochs = mne.Epochs(\n    raw=raw,\n    tmin=epochs_tmin,\n    tmax=epochs_tmax,\n    baseline=baseline,\n    reject=reject,\n    events=events,\n    event_id=event_id,\n    metadata=metadata,\n    preload=True,\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's do a final sanity check: we want to make sure that in every row, we\nactually have a stimulus. We use ``epochs.metadata`` (and not ``metadata``)\nbecause when creating the epochs, we passed the ``reject`` parameter, and\nMNE-Python always ensures that ``epochs.metadata`` stays in sync with the\navailable epochs.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "epochs.metadata.loc[epochs.metadata[\"last_stimulus\"].isna(), :]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Bummer! It seems the very first two responses were recorded before the\nfirst stimulus appeared: the values in the ``stimulus`` column are ``None``.\nThere is a very simple way to select only those epochs that **do** have a\nstimulus (i.e., are not ``None``):\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "epochs = epochs[\"last_stimulus.notna()\"]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Time to calculate the ERPs for correct  and incorrect responses.\nFor visualization, we'll only look at sensor ``FCz``, which is known to show\nthe ERN nicely in the given paradigm. We'll also create a topoplot to get an\nimpression of the average scalp potentials measured in the first 100 ms after\nan incorrect response.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "resp_erp_correct = epochs[\"response_correct\"].average()\nresp_erp_incorrect = epochs[\"not response_correct\"].average()\n\nmne.viz.plot_compare_evokeds(\n    {\"Correct Response\": resp_erp_correct, \"Incorrect Response\": resp_erp_incorrect},\n    picks=\"FCz\",\n    show_sensors=True,\n    title=\"ERPs at FCz, time-locked to response\",\n)\n\n# topoplot of average field from time 0.0-0.1 s\nfig = resp_erp_incorrect.plot_topomap(times=0.05, average=0.05, size=3)\nfig.suptitle(\"Avg. topography 0\u2013100 ms after incorrect responses\", fontsize=16)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can see a strong negative deflection immediately after incorrect\nresponses, compared to correct responses. The topoplot, too, leaves no doubt:\nwhat we're looking at is, in fact, the ERN.\n\nSome researchers suggest to construct the difference wave between ERPs for\ncorrect and incorrect responses, as it more clearly reveals signal\ndifferences, while ideally also improving the signal-to-noise ratio (under\nthe assumption that the noise level in \"correct\" and \"incorrect\" trials is\nsimilar). Let's do just that and put it into a publication-ready\nvisualization.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# difference wave: incorrect minus correct responses\nresp_erp_diff = mne.combine_evoked(\n    [resp_erp_incorrect, resp_erp_correct], weights=[1, -1]\n)\n\nfig, ax = plt.subplots()\nresp_erp_diff.plot(picks=\"FCz\", axes=ax, selectable=False, show=False)\n\n# make ERP trace bolder\nax.lines[0].set_linewidth(1.5)\n\n# add lines through origin\nax.axhline(0, ls=\"dotted\", lw=0.75, color=\"gray\")\nax.axvline(0, ls=(0, (10, 10)), lw=0.75, color=\"gray\", label=\"response trigger\")\n\n# mark trough\ntrough_time_idx = resp_erp_diff.copy().pick(\"FCz\").data.argmin()\ntrough_time = resp_erp_diff.times[trough_time_idx]\nax.axvline(trough_time, ls=(0, (10, 10)), lw=0.75, color=\"red\", label=\"max. negativity\")\n\n# legend, axis labels, title\nax.legend(loc=\"lower left\")\nax.set_xlabel(\"Time (s)\", fontweight=\"bold\")\nax.set_ylabel(\"Amplitude (\u00b5V)\", fontweight=\"bold\")\nax.set_title(\"Channel: FCz\")\nfig.suptitle(\"ERN (Difference Wave)\", fontweight=\"bold\")\n\nfig"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## References\n.. footbibliography::\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
}