P3 Auditory oddball¶
This builds on the P50 example.
Experimental design
Stimuli are 2 different auditory tones presented in random order, with higher probability (“standard, .66) and lower probability (”target”, .33).
One tone is higher frequency, the other is lower frequency. There are two blocks of trials, each tone serves as the standard one block, the target in the other.
The task is to respond to the target tones by pressing a button.
Event code tagging
The codemap patterns tag the stimulus and response event sequences for stimulus type, tone type, and response accuracy.
EEG quality control flags (log_evcodes > 0) are used to exclude EEG trials prior to analysis
[1]:
import os
import sys
from pathlib import Path
import re
import numpy as np
import pandas as pd
import mkpy
import spudtr
from matplotlib import pyplot as plt
from mkpy import mkh5
from spudtr import epf
# path wrangling for nbsphinx
if "MDE_HOME" in os.environ.keys():
MDE_HOME = Path(os.environ["MDE_HOME"])
else:
from conf import MDE_HOME
DOCS_DATA = MDE_HOME / "docs/_data"
print(os.environ["CONDA_DEFAULT_ENV"])
print(sys.version)
for pkg in [np, pd, mkpy, spudtr]:
print(pkg.__name__, pkg.__version__, pkg.__file__)
mkconda_dev_py39_053022
3.9.13 | packaged by conda-forge | (main, May 27 2022, 16:56:21)
[GCC 10.3.0]
numpy 1.21.6 /home/turbach/miniconda39/envs/mkconda_dev_py39_053022/lib/python3.9/site-packages/numpy/__init__.py
pandas 1.1.5 /home/turbach/miniconda39/envs/mkconda_dev_py39_053022/lib/python3.9/site-packages/pandas/__init__.py
mkpy 0.2.7 /mnt/cube/home/turbach/TPU_Projects/mkpy/mkpy/__init__.py
spudtr 0.1.0 /home/turbach/miniconda39/envs/mkconda_dev_py39_053022/lib/python3.9/site-packages/spudtr/__init__.py
mkh5 EEG data, event code log, header information
[2]:
# set filenames
crw = MDE_HOME / "mkdig/sub000p3.crw" # EEG recording
log = MDE_HOME / "mkdig/sub000p3.x.log" # events
yhdr = MDE_HOME / "mkpy/sub000p3.yhdr" # extra header info
# set calibration data filenames
cals_crw = MDE_HOME / "mkdig/sub000c.crw"
cals_log = MDE_HOME / "mkdig/sub000c.x.log"
cals_yhdr = MDE_HOME / "mkpy/sub000c.yhdr"
# HDF5 file with EEG recording, events, and header
p3_h5_f = DOCS_DATA / "sub000p3.h5"
print(p3_h5_f)
# convert to HDF5
p3_h5 = mkh5.mkh5(p3_h5_f)
p3_h5.reset_all()
p3_h5.create_mkdata("sub000", crw, log, yhdr)
# add calibration data
p3_h5.append_mkdata("sub000", cals_crw, cals_log, cals_yhdr)
# calibrate
pts, pulse, lo, hi, ccode = 5, 10, -40, 40, 0
p3_h5.calibrate_mkdata(
"sub000", # data group to calibrate with these cal pulses
n_points=pts, # pts to average
cal_size=pulse, # uV
lo_cursor=lo, # lo_cursor ms
hi_cursor=hi, # hi_cursor ms
cal_ccode=ccode, # condition code
)
/mnt/cube/home/turbach/TPU_Projects/mkpy_data_examples/docs/_data/sub000p3.h5
/mnt/cube/home/turbach/TPU_Projects/mkpy/mkpy/mkh5.py:3666: UserWarning: negative event code(s) found for cal condition code 0 -16384
warnings.warn(msg)
Found cals in /sub000/dblock_4
Calibrating block /sub000/dblock_0 of 5: (31232,)
Calibrating block /sub000/dblock_1 of 5: (32768,)
Calibrating block /sub000/dblock_2 of 5: (31744,)
Calibrating block /sub000/dblock_3 of 5: (32512,)
Calibrating block /sub000/dblock_4 of 5: (28416,)
codemap
[3]:
%%bash
head -26 ${MDE_HOME}/mkpy/p3_codemap.ytbl
# mkpy YAML codemap
# rows and columns must align, regexp is a mandatory column, all else is whatever.
#
# Note the regular expression negation *NON CAPTURNG* group (?!1040)
# to pattern-match "anything except 1040", i.e., an omitted response,
# as in (#10) (?!1040). Depending on the stimulus type, an ommited
# response is "correct" or "incorrect"
---
name: p3
columns:
[regexp, ccode, instrument, bin, tone, stim, accuracy, acc_type ]
rows:
- ['(#\d+)', 0, cal, 0, cal, cal, cal, cal ]
- ['(#11) 1040', 1, eeg, 3, hi, target, correct, hit ]
- ['(#10) (?!1040)', 1, eeg, 4, lo, standard, correct, cr ]
- ['(#11) (?!1040)', 1, eeg, 5, hi, target, incorrect, miss ]
- ['(#10) 1040', 1, eeg, 6, lo, standard, incorrect, fa ]
- ['(#21) 1040', 1, eeg, 9, lo, target, correct, hit ]
- ['(#20) (?!1040)', 1, eeg, 10, hi, standard, correct, cr ]
- ['(#21) (?!1040)', 1, eeg, 11, lo, target, incorrect, miss ]
- ['(#20) 1040', 1, eeg, 12, hi, standard, incorrect, fa ]
1. get_event_table(codemap)
[4]:
p3_event_table = p3_h5.get_event_table(MDE_HOME / "mkpy/p3_codemap.ytbl")
searching codes in: sub000/dblock_0
searching codes in: sub000/dblock_1
searching codes in: sub000/dblock_2
searching codes in: sub000/dblock_3
searching codes in: sub000/dblock_4
/mnt/cube/home/turbach/TPU_Projects/mkpy/mkpy/mkh5.py:1059: UserWarning:
As of mkpy 0.2.0 to match events with a codemap regexp pattern, the
ccode column in p3_codemap.ytbl must also match the log_ccode
in the datablock. If this behavior is not desired, delete or rename
the ccode column in the codemap.
warnings.warn(msg)
inspect the event table
[5]:
print(p3_event_table.shape)
print(p3_event_table.columns)
# select some columns to show
example_columns = [
"dblock_path", "dblock_ticks", "log_evcodes", "log_ccodes", "log_flags",
"regexp", "match_code", "instrument", "bin", "tone", "stim", "acc_type",
]
# first few stimulus events
display(p3_event_table[example_columns].head())
# last few calibration pulse events
display(p3_event_table[example_columns].tail())
pd.crosstab(
[
p3_event_table.instrument,
p3_event_table.ccode,
p3_event_table.stim,
p3_event_table.accuracy,
p3_event_table.acc_type,
p3_event_table.tone,
p3_event_table.bin,
],
[
p3_event_table.log_flags
],
margins=True
)
(601, 31)
Index(['data_group', 'dblock_path', 'dblock_tick_idx', 'dblock_ticks',
'crw_ticks', 'raw_evcodes', 'log_evcodes', 'log_ccodes', 'log_flags',
'epoch_match_tick_delta', 'epoch_ticks', 'dblock_srate', 'match_group',
'idx', 'dlim', 'anchor_str', 'match_str', 'anchor_code', 'match_code',
'anchor_tick', 'match_tick', 'anchor_tick_delta', 'is_anchor', 'regexp',
'ccode', 'instrument', 'bin', 'tone', 'stim', 'accuracy', 'acc_type'],
dtype='object')
| dblock_path | dblock_ticks | log_evcodes | log_ccodes | log_flags | regexp | match_code | instrument | bin | tone | stim | acc_type | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 125 | sub000/dblock_0 | 1468 | 11 | 1 | 0 | (#11) 1040 | 11 | eeg | 3 | hi | target | hit |
| 126 | sub000/dblock_0 | 2396 | 11 | 1 | 0 | (#11) 1040 | 11 | eeg | 3 | hi | target | hit |
| 127 | sub000/dblock_0 | 2943 | 11 | 1 | 0 | (#11) 1040 | 11 | eeg | 3 | hi | target | hit |
| 128 | sub000/dblock_0 | 3308 | 11 | 1 | 0 | (#11) 1040 | 11 | eeg | 3 | hi | target | hit |
| 129 | sub000/dblock_0 | 4574 | 11 | 1 | 0 | (#11) 1040 | 11 | eeg | 3 | hi | target | hit |
| dblock_path | dblock_ticks | log_evcodes | log_ccodes | log_flags | regexp | match_code | instrument | bin | tone | stim | acc_type | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1088 | sub000/dblock_4 | 27315 | 4 | 0 | 0 | (#\d+) | 4 | cal | 0 | cal | cal | cal |
| 1089 | sub000/dblock_4 | 27444 | 2 | 0 | 0 | (#\d+) | 2 | cal | 0 | cal | cal | cal |
| 1090 | sub000/dblock_4 | 27573 | 3 | 0 | 0 | (#\d+) | 3 | cal | 0 | cal | cal | cal |
| 1091 | sub000/dblock_4 | 27703 | 4 | 0 | 0 | (#\d+) | 4 | cal | 0 | cal | cal | cal |
| 1092 | sub000/dblock_4 | 27832 | 2 | 0 | 0 | (#\d+) | 2 | cal | 0 | cal | cal | cal |
[5]:
| log_flags | 0 | 32 | 48 | 64 | All | ||||||
|---|---|---|---|---|---|---|---|---|---|---|---|
| instrument | ccode | stim | accuracy | acc_type | tone | bin | |||||
| cal | 0 | cal | cal | cal | cal | 0 | 208 | 0 | 0 | 1 | 209 |
| eeg | 1 | standard | correct | cr | hi | 10 | 142 | 3 | 1 | 0 | 146 |
| lo | 4 | 138 | 7 | 1 | 0 | 146 | |||||
| target | correct | hit | hi | 3 | 23 | 12 | 15 | 0 | 50 | ||
| lo | 9 | 31 | 11 | 8 | 0 | 50 | |||||
| All | 542 | 33 | 25 | 1 | 601 |
2. set_epochs(name, pre, post)
Stash the the event table events and tags with a name and epoch interval boundaries.
[6]:
p3_h5.set_epochs("ms1500", p3_event_table, -750, 750)
Sanitizing event table data types for mkh5 epochs table ...
/mnt/cube/home/turbach/TPU_Projects/mkpy/mkpy/mkh5.py:1515: UserWarning: data error: pre-stimulus interval is out of bounds left ... skipping epoch (392, b'sub000', b'sub000/dblock_4', 0, 26, 26, 2, 2, 0, 64, -187, 375, 250., 1, 0, 1, b'2', b'2', 2, 2, 26, 26, 0, True, b'(#\\d+)', 0, b'cal', 0, b'cal', b'cal', b'cal', b'cal', 0, 0, 0, 26, -187, 375)
warnings.warn(
3. export_epochs(name, file_format=…)
Export the time-stamped, tagged, fixed-length segments of EEG data for analysis as HDF5, pandas HDF5 or feather
[7]:
p3_epochs_f = DOCS_DATA / "sub000p3.ms1500.epochs.feather"
p3_h5.export_epochs("ms1500", p3_epochs_f, file_format="feather")
analyze the epochs
Is there a difference between standard and targets?
If so, is the effect the same for high and low tone standards and targers?
[8]:
p3_epochs = pd.read_feather(p3_epochs_f)
# sanitize Epoch_idx for pandas index
p3_epochs['epoch_id'] = p3_epochs['epoch_id'].astype('int')
# 26 cap + EOG and A2 channel labels
p3_chans = [col for col in p3_epochs.columns if re.match(r"^[MLRAHlr]\w{1,3}$", col)]
print(f"n channels: {len(p3_chans)}\nchannels:{p3_chans}")
# lookup the epochs with stim events at time == 0) and clean artifact flags (= 0)
good_epoch_ids = p3_epochs.query("match_time == 0 and log_flags == 0").epoch_id
# select just the good epochs
p3_epochs = p3_epochs.query("epoch_id in @good_epoch_ids").copy()
# check the good event counts after dropping artifacts
p3_good_events = p3_epochs.query("match_time == 0")
print("After excluding EEG artifacts")
pd.crosstab(
[
p3_good_events.instrument,
p3_good_events.ccode,
p3_good_events.stim,
p3_good_events.accuracy,
p3_good_events.acc_type,
p3_good_events.tone,
p3_good_events.bin,
p3_good_events.log_evcodes,
],
[
p3_good_events.log_flags
],
margins=True
)
n channels: 32
channels:['lle', 'lhz', 'MiPf', 'LLPf', 'RLPf', 'LMPf', 'RMPf', 'LDFr', 'RDFr', 'LLFr', 'RLFr', 'LMFr', 'RMFr', 'LMCe', 'RMCe', 'MiCe', 'MiPa', 'LDCe', 'RDCe', 'LDPa', 'RDPa', 'LMOc', 'RMOc', 'LLTe', 'RLTe', 'LLOc', 'RLOc', 'MiOc', 'A2', 'HEOG', 'rle', 'rhz']
After excluding EEG artifacts
[8]:
| log_flags | 0 | All | |||||||
|---|---|---|---|---|---|---|---|---|---|
| instrument | ccode | stim | accuracy | acc_type | tone | bin | log_evcodes | ||
| cal | 0 | cal | cal | cal | cal | 0 | 1 | 49 | 49 |
| 2 | 47 | 47 | |||||||
| 3 | 56 | 56 | |||||||
| 4 | 56 | 56 | |||||||
| eeg | 1 | standard | correct | cr | hi | 10 | 20 | 142 | 142 |
| lo | 4 | 10 | 138 | 138 | |||||
| target | correct | hit | hi | 3 | 11 | 23 | 23 | ||
| lo | 9 | 21 | 31 | 31 | |||||
| All | 542 | 542 |
[9]:
# for illustration ...
midline = ["MiPf", "MiCe", "MiPa", "MiOc"]
# select COLUMNS: epoch index, timestamps, event tags, and midline EEG columns
midline_epochs = p3_epochs[["epoch_id", "match_time", "stim", "tone"] + midline]
# select ROWS: use tag values to pick and choose
midline_epochs = midline_epochs.query("stim in ['standard', 'target']")
# center each channel
midline_epochs = epf.center_eeg(
midline_epochs,
midline,
-750, 0,
epoch_id="epoch_id",
time="match_time"
)
time-domain average midline ERPs
[10]:
# compute domain average by stim type
midline_erps = midline_epochs.groupby(
["stim", "match_time"]
).mean().reset_index()
# plot
f, axs = plt.subplots(4, 1, figsize=(8,12), sharey=True)
for (stim), erp in midline_erps.groupby(['stim']):
for axi, chan in enumerate(midline):
# mark onset
axs[axi].axvline(0, color='gray')
axs[axi].plot(
erp["match_time"],
erp[chan],
label=f"{stim}",
lw=2,
)
# channel
axs[axi].set_title(chan, loc='left')
axs[0].legend()
axs[0].set_title("Auditory oddball P300")
f.tight_layout()
[11]:
# average by stim type
midline_erps = midline_epochs.groupby(["stim", "tone", "match_time"]).mean().reset_index()
# plot
f, axs = plt.subplots(4, 1, figsize=(8,12), sharey=True)
for (stim, tone), erp in midline_erps.groupby(['stim', "tone"]):
for axi, chan in enumerate(midline):
# mark onset
axs[axi].axvline(0, color='gray')
axs[axi].plot(
erp["match_time"],
erp[chan],
label=f"{stim} {tone}",
lw=2,
)
# channel
axs[axi].set_title(chan, loc='left')
axs[0].legend()
axs[0].set_title("Auditory oddball P300")
f.tight_layout()
