nb-7-1-moabb · Decoding and BCI (L7.1)¶
Lesson L7.1 · Level 7 · Status draft — for expert review; uncertain points carry TODO(confirm).
Three pipelines, two paradigms, three places to put the fold boundary. The question the lesson ends on is "which pipeline wins cross-subject and by how much", and the honest answer turns out to be about the and by how much.
- Motor imagery on
ds-bci-iv-2a— CSP + LDA, a Riemannian tangent-space classifier, and filter-bank CSP, evaluated within session, across sessions (different days) and across subjects. - P300 on
ds-brain-invaders— the same three families adapted to an evoked response, within subject and across subjects. - Leakage, measured rather than asserted: a random split pooled over two sessions against the train-on-day-one / test-on-day-two split the data set's own design gives, and a check of whether the EOG channels the authors forbid as features actually carry the label.
Data.
ds-bci-iv-2a— BCI Competition IV data set 2a / BNCI Horizon 001-2014, Tangermann et al. (2012), Review of the BCI Competition IV, Frontiers in Neuroscience 6, 55, DOI 10.3389/fnins.2012.00055. No dataset DOI exists (TODO(confirm)). Licence CC BY-ND 4.0,access: open. 9 subjects, 2 sessions on different days, 22 EEG + 3 EOG at 250 Hz.ds-brain-invaders— Brain Invaders bi2014a, Korczowski et al. (2019), GIPSA-lab research report HAL hal-02171575; dataset DOI 10.5281/zenodo.3266223. Licence CC BY 4.0,access: open. 16 dry electrodes at 512 Hz.
No derivatives ship from ds-bci-iv-2a. CC BY-ND is a no-derivatives licence and every snippet this site
ships is a derivative, so this notebook downloads the data, computes on it and publishes tables and figures of
its own results — a benchmark table is a report of findings, not a redistribution of the recording — and ships
nothing cut from the signal. w-csp-explorer keeps using ds-eegbci. The licence was verified from primary
material; see site/notes/data-p4-licences.md.
Two instructions from the data set's own documentation are followed and one of them is tested. The EOG channels "must not be used for classification", and the two sessions are on different days so a random pooled split leaks. Section 4 measures both rather than taking them on trust.
The compact CNN spec §6 asks for is not run. braindecode and torch are not in the pinned stack; filter-bank
CSP stands in as the third pipeline family, and section 0 prints exactly what that substitution does and does not
answer. No published EEGNet accuracy is quoted from memory.
# Setup: dependencies, the shared helpers, non-interactive plotting, a quiet downloader.
import importlib.util
import subprocess
import sys
import time
import warnings
from pathlib import Path
_needed = ("mne", "scipy", "matplotlib", "pandas", "pooch", "sklearn", "moabb", "pyriemann")
_missing = [p for p in _needed if importlib.util.find_spec(p) is None]
if _missing:
_req = next((d / "requirements.txt" for d in (Path.cwd(), *Path.cwd().parents)
if (d / "requirements.txt").exists()), None)
_cmd = [sys.executable, "-m", "pip", "install", "-q"]
_cmd += ["-r", str(_req)] if _req else ["mne==1.10.2", "moabb==1.7.2", "pyriemann", "scikit-learn"]
subprocess.check_call(_cmd)
_shared = next((d / "_shared" for d in (Path.cwd(), *Path.cwd().parents)
if (d / "_shared" / "helpers_l7.py").exists()), None)
if _shared is None:
raise FileNotFoundError("start the kernel in notebooks/L7/ (or notebooks/) so that _shared/helpers_l7.py is found")
sys.path.insert(0, str(_shared))
import helpers
import helpers_l6 as L6
import helpers_l7 as L7
import matplotlib.pyplot as plt
import numpy as np
import mne
import pooch
from scipy import stats
from sklearn.base import clone
from sklearn.metrics import confusion_matrix, get_scorer
from sklearn.model_selection import StratifiedKFold, cross_val_predict, cross_val_score
mne.set_log_level("WARNING")
# Silence the downloader at the source rather than scrubbing cache paths out of stored outputs
# afterwards: a re-execution would put them straight back (Phase 3 note in notebooks/README.md).
pooch.get_logger().setLevel("WARNING")
plt.rcParams["figure.dpi"] = 72
SEED = L7.SEED
MI_SUBJECTS = list(range(1, 10)) # every ds-bci-iv-2a subject
P300_SUBJECTS = [1, 2, 3, 4, 5, 6] # a documented subset of bi2014a's 64
FOLDS = 5
print(f"MNE {mne.__version__}; helpers_l7 imported from notebooks/_shared")
print(f"moabb download provider pinned to: {L7.use_upstream_moabb()}")
print(f"seed {SEED}; {len(MI_SUBJECTS)} motor-imagery subjects, {len(P300_SUBJECTS)} P300 subjects, "
f"{FOLDS}-fold cross-validation")
0 · What may ship, what the data set's own rules are, and what is missing¶
Read before any number: a licence decides what leaves this notebook, and two of the data set's own instructions decide how the pipelines are allowed to be built.
for ds in ("ds-bci-iv-2a", "ds-brain-invaders"):
print(L7.dataset_line(ds))
print()
print("The data set's own instructions (from its description document, not from its licence):")
for k, v in L7.BCI_IV_2A_RULES.items():
print(f" [{k}] {v}")
print("\nThe third pipeline spec section 6 asks for:")
for k in ("asked_for", "status", "why", "substitute", "what_it_does_not_answer"):
print(f" [{k}] {L7.CNN_ARM[k]}")
1 · One extraction, three pipelines¶
The three pipelines see the same trials, the same channels, the same window and the same folds. They differ in the classifier and in nothing else, which is the only way a benchmark table means anything.
Band-pass filtering does not use the labels, so it cannot leak and is done once, outside the cross-validation
loop. CSP, the filter bank's CSPs and the feature selection do use the labels and are re-fitted inside every
training fold. That line is the whole of pf-decoding-leakage and section 4 measures what crossing it costs.
Each subject is downloaded, epoched and deleted before the next one is fetched, so peak disk is about 87 MB rather than the 780 MB the whole set would take.
for k, v in L7.MI_SPEC.items():
print(f" {k:24s} {v}")
L7.disk_report("before any download",
folders={"BNCI cache": L7.moabb_root(), "NEMAR cache": L7.nemar_root()})
already_on_disk = [f for s in MI_SUBJECTS for f in L7.bnci_subject_files(s)]
print(f"{len(already_on_disk)} of this dataset's files were already cached before this run; "
f"they are kept, everything this run fetches is deleted")
t0 = time.time()
mi = {}
for s in MI_SUBJECTS:
mi[s] = L7.load_bci_iv_2a_subject(s, keep=already_on_disk, four_class=(s == MI_SUBJECTS[0]),
default_picks=(s == MI_SUBJECTS[0]))
print(f"\n{len(mi)} subjects in {time.time() - t0:.0f} s; "
f"{sum(d['downloaded_mb'] for d in mi.values()):.0f} MB fetched, "
f"{sum(d['freed_mb'] for d in mi.values()):.0f} MB deleted")
assert len(mi) == len(MI_SUBJECTS), "a subject is missing: the cohort must be complete before any mean"
L7.disk_report("after the download loop",
folders={"BNCI cache": L7.moabb_root(), "NEMAR cache": L7.nemar_root()})
The EOG channels are absent from the features, and that is checked rather than assumed¶
The data set's description document says the three EOG channels "must not be used for classification". moabb's
paradigm already drops them: asked for the default picks it returns the 22 EEG channels and nothing else. This
notebook asks for all 25 on purpose, so that section 4 can test what is in the EOG, and every pipeline that
reports a score is fitted on the 22.
d0 = mi[MI_SUBJECTS[0]]
print(f"moabb's default picks (channels=None): {len(d0['default_ch_names'])} channels")
print(f" {d0['default_ch_names']}")
print(f"EOG among them: {[c for c in d0['default_ch_names'] if c.upper().startswith('EOG')] or 'none'}")
print(f"\nwhat the pipelines are fitted on: {len(d0['ch_names'])} channels, "
f"{[c for c in d0['ch_names'] if c.upper().startswith('EOG')] or 'no EOG'}")
print(f"what section 4 tests separately: {d0['eog_names']}")
assert not any(c.upper().startswith("EOG") for c in d0["ch_names"])
print(f"\nepochs per subject: {d0['X'].shape[0]} two-class trials, "
f"{d0['X'].shape[1]} channels x {d0['n_times']} samples at {d0['sfreq']:.0f} Hz")
print(f"filter bank: {d0['X_fb'].shape[1]} sub-bands {L7.MI_SPEC['fb_bands_hz']}")
print(f"arrays held in memory: "
f"{sum(v.nbytes for d in mi.values() for v in d.values() if isinstance(v, np.ndarray)) / 1e6:.0f} MB")
2 · Three evaluations¶
Where the fold boundary falls changes the number more than the choice of classifier does.
for k, v in L7.EVALUATION_NOTES.items():
print(f"[{k}]\n {v}\n")
pipelines = L7.build_mi_pipelines()
for name, p in pipelines.items():
print(f"{name:24s} {' -> '.join(s for s, _ in p.steps)}")
rows = []
t0 = time.time()
rows += L7.within_session_scores(mi, pipelines, folds=FOLDS, seed=SEED, progress=False)
print(f"within-session {time.time() - t0:5.0f} s")
t1 = time.time()
rows += L7.cross_session_scores(mi, pipelines, progress=False)
print(f"cross-session {time.time() - t1:5.0f} s")
t1 = time.time()
rows += L7.cross_subject_scores(mi, pipelines, progress=False)
print(f"cross-subject {time.time() - t1:5.0f} s")
mi_df = L7.benchmark_frame(rows)
print(f"total {time.time() - t0:5.0f} s\n")
L7.print_benchmark(mi_df, chance=0.5,
title="ds-bci-iv-2a, left hand vs right hand, accuracy (dimensionless)")
lo, hi = L6.chance_band(int(mi[MI_SUBJECTS[0]]["X"].shape[0]), 0.5)
print(f"binomial 95 % band around chance for {mi[MI_SUBJECTS[0]]['X'].shape[0]} trials: {lo:.4f}-{hi:.4f}")
fig, _ = L7.plot_benchmark(mi_df, chance=0.5, scoring="accuracy",
title="ds-bci-iv-2a: decoding accuracy by pipeline and evaluation "
"(dimensionless, 9 subjects)")
plt.show()
Which pipeline wins, and by how much¶
The margin between pipelines is reported paired across subjects — every pipeline saw the same people and the same trials — because the unpaired spread across subjects is an order of magnitude larger than any difference between the pipelines, and quoting one without the other would say the opposite of what the data says.
mi_winners = {}
for ev in ("within-session", "cross-session", "cross-subject"):
w = L7.winner(mi_df, ev)
mi_winners[ev] = w
print(f"[{ev}]")
print(f" ranking: " + ", ".join(f"{k} {v:.4f}" for k, v in w["ranking"]))
print(f" {w['winner']} beats {w['runner_up']} by {w['margin']:+.4f} "
f"(paired mean {w['margin_mean_paired']:+.4f}, paired SD {w['margin_sd_paired']:.4f}, "
f"t({w['n_subjects'] - 1}) = {w['t']:.2f}, p = {w['p']:.4f})")
spread = mi_df[mi_df['evaluation'] == ev].groupby('subject', observed=True)['score'].mean().std(ddof=1)
print(f" between-subject SD of the mean score: {spread:.4f} "
f"({spread / max(abs(w['margin']), 1e-9):.0f}x the winning margin)\n")
3 · The four-class design, and why chance is not a constant¶
The data set is a four-class paradigm — left hand, right hand, both feet, tongue. The benchmark above uses
the two hand classes, so chance is 0.5 and the numbers are directly comparable with the ds-eegbci left/right
pipeline of L5.7 and L6.5. Running the same estimator on all four classes on one subject shows what "accuracy"
means when the design changes under it, and where the errors go.
d = mi[MI_SUBJECTS[0]]
m = d["session4"] == "0train"
X4, y4 = d["X4"][m].astype(np.float64), d["y4"][m]
cv = StratifiedKFold(n_splits=FOLDS, shuffle=True, random_state=SEED)
with warnings.catch_warnings():
warnings.simplefilter("ignore")
sc4 = cross_val_score(L7.make_csp_lda(8), X4, y4, cv=cv, scoring="accuracy")
pred4 = cross_val_predict(L7.make_csp_lda(8), X4, y4, cv=cv)
classes = d["classes4"]
print(f"subject {d['subject']}, session T, {len(y4)} trials, {len(classes)} classes {classes}")
print(f" four-class accuracy {sc4.mean():.4f} (chance 0.25)")
sc2 = mi_df[(mi_df['evaluation'] == 'within-session') & (mi_df['subject'] == d['subject'])
& (mi_df['pipeline'] == 'CSP + LDA') & (mi_df['session'] == '0train')]['score']
print(f" two-class accuracy on the same session {float(sc2.iloc[0]):.4f} (chance 0.50)")
cm = confusion_matrix(y4, pred4)
print("\nconfusion (rows = true, columns = predicted):")
print(f"{'':12s}" + "".join(f"{c:>12s}" for c in classes))
for i, c in enumerate(classes):
print(f"{c:12s}" + "".join(f"{int(v):12d}" for v in cm[i]))
4 · Leakage, measured twice¶
pf-decoding-leakage in two forms the data set warns about by itself.
4.1 · A random split pooled over two sessions¶
The two sessions were recorded on different days. Pooling them and cross-validating at random puts trials from the same day, the same cap placement and the same impedances on both sides of every fold boundary. Train on session T and test on session E and nothing crosses. The same estimator, the same trials, the same seed: only the split moves.
cv = StratifiedKFold(n_splits=FOLDS, shuffle=True, random_state=SEED)
acc = get_scorer("accuracy")
pooled, sessionwise = [], []
print(f"{'subject':>8s} {'pooled random':>14s} {'train T / test E':>17s} {'gap':>8s}")
for s in MI_SUBJECTS:
d = mi[s]
X, y, sess = d["X"].astype(np.float64), d["y"], d["session"]
with warnings.catch_warnings():
warnings.simplefilter("ignore")
p = float(cross_val_score(L7.make_csp_lda(8), X, y, cv=cv, scoring="accuracy").mean())
tr, te = sess == "0train", sess == "1test"
est = clone(L7.make_csp_lda(8)).fit(X[tr], y[tr])
c = float(acc(est, X[te], y[te]))
pooled.append(p)
sessionwise.append(c)
print(f"{s:8d} {p:14.4f} {c:17.4f} {p - c:+8.4f}")
pooled, sessionwise = np.asarray(pooled), np.asarray(sessionwise)
t, pv = stats.ttest_rel(pooled, sessionwise)
print(f"\nmean {pooled.mean():.4f} pooled against {sessionwise.mean():.4f} cross-session; "
f"inflation {pooled.mean() - sessionwise.mean():+.4f}")
print(f"paired SD {np.std(pooled - sessionwise, ddof=1):.4f}, t({len(pooled) - 1}) = {t:.2f}, p = {pv:.4f}; "
f"largest single-subject gap {np.max(pooled - sessionwise):+.4f}")
4.2 · Are the forbidden channels actually informative?¶
The data set says the EOG channels must not be used for classification. That is an instruction, and it is worth
knowing whether it is also a live risk on this data: a pipeline that leaves the EOG in its feature set is
decoding eye movement, and pf-decoding-leakage is about exactly that kind of shortcut. The test is the same
estimator on the three EOG channels alone.
eog_acc = []
for s in MI_SUBJECTS:
d = mi[s]
with warnings.catch_warnings():
warnings.simplefilter("ignore")
v = float(cross_val_score(L7.make_csp_lda(2), d["X_eog"].astype(np.float64), d["y"],
cv=cv, scoring="accuracy").mean())
eog_acc.append(v)
eog_acc = np.asarray(eog_acc)
n_trials = int(mi[MI_SUBJECTS[0]]["X"].shape[0])
lo, hi = L6.chance_band(n_trials, 0.5)
above = eog_acc > hi
print(f"{'subject':>8s} {'EOG-only accuracy':>18s}")
for s, v in zip(MI_SUBJECTS, eog_acc):
print(f"{s:8d} {v:18.4f}{' above the chance band' if v > hi else ''}")
print(f"\nmean {eog_acc.mean():.4f}, range {eog_acc.min():.4f}-{eog_acc.max():.4f}")
print(f"binomial 95 % band around chance for {n_trials} trials: {lo:.4f}-{hi:.4f}")
print(f"{int(above.sum())} of {len(eog_acc)} subjects are above it — the prohibition is not a formality, "
f"and it is also not a licence to claim the EEG result is an EOG artifact")
5 · The P300 paradigm on dry electrodes¶
The same three families, adapted to an evoked response: a supervised spatial filter with a linear classifier, the
same spatial filtering into a Riemannian tangent space, and no spatial filter at all. moabb collapses this
data set's game sessions into one, so there is a within-subject and a cross-subject evaluation and no
cross-session one. The score is ROC AUC, not accuracy, because the oddball is about one target to five
non-targets and accuracy would reward a classifier that always answered "non-target".
for k, v in L7.P300_SPEC.items():
print(f" {k:20s} {v}")
bi_root = helpers.download_root() / "MNE-braininvaders2014a-data"
bi_pre = [p for p in bi_root.rglob("*") if p.is_file()] if bi_root.exists() else []
L7.disk_report("before the P300 downloads", folders={"bi2014a cache": bi_root})
t0 = time.time()
p300 = {}
for s in P300_SUBJECTS:
p300[s] = L7.load_bi2014a_subject(s, keep=bi_pre)
print(f"\n{len(p300)} subjects in {time.time() - t0:.0f} s; "
f"{sum(d['downloaded_mb'] for d in p300.values()):.0f} MB fetched, "
f"{sum(d['freed_mb'] for d in p300.values()):.0f} MB deleted")
assert len(p300) == len(P300_SUBJECTS)
L7.disk_report("after the P300 downloads", folders={"bi2014a cache": bi_root})
p300_pipelines = L7.build_p300_pipelines()
p300_inputs = {k: "X" for k in p300_pipelines}
rows = []
t0 = time.time()
rows += L7.within_session_scores(p300, p300_pipelines, inputs=p300_inputs, folds=FOLDS, seed=SEED,
scoring="roc_auc", session_key="__single__", progress=False)
rows += L7.cross_subject_scores(p300, p300_pipelines, inputs=p300_inputs, scoring="roc_auc",
progress=False)
p300_df = L7.benchmark_frame(rows)
print(f"{time.time() - t0:.0f} s\n")
L7.print_benchmark(p300_df, scoring="roc_auc", chance=0.5,
title="ds-brain-invaders bi2014a, target vs non-target, ROC AUC (dimensionless)")
p300_winners = {}
for ev in p300_df["evaluation"].cat.categories:
w = L7.winner(p300_df, ev)
p300_winners[ev] = w
print(f"\n[{ev}] " + ", ".join(f"{k} {v:.4f}" for k, v in w["ranking"]))
print(f" {w['winner']} beats {w['runner_up']} by {w['margin']:+.4f} "
f"(paired SD {w['margin_sd_paired']:.4f}, t({w['n_subjects'] - 1}) = {w['t']:.2f}, "
f"p = {w['p']:.4f})")
fig, _ = L7.plot_benchmark(p300_df, chance=0.5, scoring="roc_auc",
title="ds-brain-invaders bi2014a: P300 decoding AUC by pipeline and "
"evaluation (dimensionless, 6 subjects)")
plt.show()
6 · What the lesson's exercise asks for¶
"Which pipeline wins cross-subject, and by how much?" — printed below for both paradigms, with the paired margin, its spread and the between-subject spread beside it, because the by how much is the part that decides whether the ranking means anything.
print("=" * 96)
print("nb-7-1-moabb — the numbers the L7.1 exercise asks for")
print("=" * 96)
print("\n[1] ds-bci-iv-2a, left vs right hand, accuracy (chance 0.5), 9 subjects, "
f"{L7.MI_SPEC['analysis_band_hz'][0]:g}-{L7.MI_SPEC['analysis_band_hz'][1]:g} Hz, "
f"{L7.MI_SPEC['analysis_window_s'][0]:g}-{L7.MI_SPEC['analysis_window_s'][1]:g} s after the cue")
for ev in ("within-session", "cross-session", "cross-subject"):
g = mi_df[mi_df["evaluation"] == ev]
print(f" {ev:15s} " + " ".join(
f"{p}: {g[g['pipeline'] == p]['score'].mean():.4f}" for p in L7.MI_INPUT))
w = mi_winners["cross-subject"]
print(f"\n[2] CROSS-SUBJECT WINNER (motor imagery): {w['winner']} at {w['winner_mean']:.4f}, "
f"ahead of {w['runner_up']} at {w['runner_up_mean']:.4f}")
print(f" margin {w['margin']:+.4f} ({100 * w['margin']:+.2f} percentage points); "
f"paired SD across the 9 subjects {w['margin_sd_paired']:.4f}; "
f"t({w['n_subjects'] - 1}) = {w['t']:.2f}, p = {w['p']:.4f}")
print(f" ranking: " + ", ".join(f"{k} {v:.4f}" for k, v in w["ranking"]))
print(f" READ THIS BEFORE QUOTING THE WINNER: the margin is "
f"{100 * w['margin']:.1f} percentage points and the paired standard deviation is "
f"{100 * w['margin_sd_paired']:.1f}. The three pipelines are not distinguishable "
f"cross-subject on 9 subjects. The finding is the DROP, not the ranking.")
wi = mi_winners["within-session"]
print(f"\n[3] The drop that is the real result: within-session "
f"{mi_df[mi_df['evaluation'] == 'within-session']['score'].mean():.4f} -> cross-session "
f"{mi_df[mi_df['evaluation'] == 'cross-session']['score'].mean():.4f} -> cross-subject "
f"{mi_df[mi_df['evaluation'] == 'cross-subject']['score'].mean():.4f} "
f"(mean over the three pipelines)")
pw = p300_winners["cross-subject"]
print(f"\n[4] CROSS-SUBJECT WINNER (P300, ROC AUC): {pw['winner']} at {pw['winner_mean']:.4f}, "
f"ahead of {pw['runner_up']} at {pw['runner_up_mean']:.4f}, margin {pw['margin']:+.4f} "
f"(paired SD {pw['margin_sd_paired']:.4f}, p = {pw['p']:.4f})")
pwi = p300_winners["within-session"]
print(f" within subject the order is different and the margin IS separable: {pwi['winner']} at "
f"{pwi['winner_mean']:.4f} over {pwi['runner_up']} at {pwi['runner_up_mean']:.4f} "
f"(margin {pwi['margin']:+.4f}, paired SD {pwi['margin_sd_paired']:.4f}, p = {pwi['p']:.4f}); "
f"{pwi['winner']} falls to {float(p300_df[(p300_df['evaluation'] == 'cross-subject') & (p300_df['pipeline'] == pwi['winner'])]['score'].mean()):.4f} "
f"cross-subject")
print(f"\n[5] Leakage — a pooled random split over two days: {pooled.mean():.4f} against "
f"{sessionwise.mean():.4f} for train-on-T / test-on-E, inflation "
f"{pooled.mean() - sessionwise.mean():+.4f} "
f"(paired SD {np.std(pooled - sessionwise, ddof=1):.4f}, p = {pv:.4f}, "
f"largest subject {np.max(pooled - sessionwise):+.4f})")
print(f"[6] Leakage — EOG-only decoding, the channels the data set forbids as features: mean "
f"{eog_acc.mean():.4f}, range {eog_acc.min():.4f}-{eog_acc.max():.4f}; "
f"{int(above.sum())} of {len(eog_acc)} subjects above the {lo:.4f}-{hi:.4f} chance band")
print(f"[7] Four classes instead of two, subject {mi[MI_SUBJECTS[0]]['subject']}, session T: "
f"{sc4.mean():.4f} against chance 0.25 (two-class {float(sc2.iloc[0]):.4f} against chance 0.50)")
print(f"\n[8] NOT ANSWERED: {L7.CNN_ARM['what_it_does_not_answer']}")
print("=" * 96)
# Nothing this notebook downloaded is left on disk: every loader deletes its subject in a finally,
# and the counts above are the evidence. This is the final check.
L7.disk_report("at the end",
folders={"BNCI cache": L7.moabb_root(), "NEMAR cache": L7.nemar_root(),
"bi2014a cache": bi_root})
left = [f.name for s in MI_SUBJECTS for f in L7.bnci_subject_files(s)]
print(f"ds-bci-iv-2a files still on disk: {sorted(left) or 'none'} "
f"(pre-existing before this run: {sorted(f.name for f in already_on_disk) or 'none'})")
assert set(left) <= {f.name for f in already_on_disk}, "this run left a download behind"
print("no download made by this run was left behind")