Learning curves from real training subjects

Estimated reading time:7 minutes

How does adding one training participant change validation performance? Use nested training-subject subsets and keep subject 3 fixed for validation. Two training sizes illustrate the procedure, not a scaling law.

Data: Nakanishi2015, NEMAR nm000118, subjects 1–3, session 0, run 0: approximately 21.1 MB on first download. Set EEGDASH_CACHE_DIR to reuse the cache. This processed release already includes filtering, downsampling and latency handling; do not add another latency correction. See the source study and NEMAR release.

Prerequisites: tutorial 51’s subject-disjoint fit and evaluation, plus the spectral baseline from tutorial 12. No saved model or earlier output file is needed. The output is a two-point validation curve, with both participant count and trial count reported so the horizontal axis is unambiguous.

1. Select a small, explicit cohort

Filtering subjects, session and run bounds the download. Cropping after opening a recording would reduce computation but not its download size.

import os
from functools import partial
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from braindecode.preprocessing import create_windows_from_events
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import balanced_accuracy_score
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

from eegdash import EEGDashDataset
from eegdash.features import (
    FeatureExtractor,
    extract_features,
    spectral_bands_power,
    spectral_preprocessor,
)

cache_dir = Path(os.environ.get("EEGDASH_CACHE_DIR", ".eegdash_cache"))
subjects = ["1", "2", "3"]
dataset = EEGDashDataset(
    cache_dir=cache_dir,
    dataset="nm000118",
    subject=subjects,
    session="0",
    run="0",
    task="ssvep",
    n_jobs=1,
)
assert len(dataset.datasets) == len(subjects), "Expected one recording per subject"
print(dataset.description[["subject", "session", "run"]])

  subject session run
0       1       0   0
1       2       0   0
2       3       0   0

2. Inspect real annotations and verify the signal contract

Accessing raw downloads that recording. Annotation names identify the attended stimulus frequency in Hz; they supply every classification label. All participants must have the same channel order and sampling frequency.

raw = dataset.datasets[0].raw
sfreq = raw.info["sfreq"]
channel_names = raw.ch_names
class_names = sorted(set(raw.annotations.description), key=float)
mapping = {name: index for index, name in enumerate(class_names)}
assert len(mapping) == 12, "Expected the twelve SSVEP stimulus frequencies"
for recording in dataset.datasets:
    recording_raw = recording.raw
    assert recording_raw.ch_names == channel_names
    assert recording_raw.info["sfreq"] == sfreq
    assert set(recording_raw.annotations.description) == set(mapping)
print(f"Channels: {channel_names}; sampling frequency: {sfreq} Hz")
print("Stimulus frequencies (Hz):", class_names)

Channels: ['PO7', 'PO3', 'POz', 'PO4', 'PO8', 'O1', 'Oz', 'O2']; sampling frequency: 256.0 Hz
Stimulus frequencies (Hz): ['9.25', '9.75', '10.25', '10.75', '11.25', '11.75', '12.25', '12.75', '13.25', '13.75', '14.25', '14.75']

3. Make one four-second window per annotated trial

Each annotated interval is 4.15 seconds. Keep its first four seconds and discard the remainder. Explicit size and stride avoid overlapping windows or extending the epoch beyond the recorded event duration. At 256 Hz, four seconds contain 1,024 samples. The resulting array has axes (540 trials, 8 EEG channels, 1,024 samples), with volt-valued data. The metadata has one row per array row. target is a class index, not a frequency in Hz; mapping is the explicit conversion between them. The source contains 15 trials of each of the 12 frequencies per person. A missing class is a data-contract failure, not a reason to relabel trials.

window_size = int(4 * sfreq)
windows = create_windows_from_events(
    dataset,
    mapping=mapping,
    trial_start_offset_samples=0,
    trial_stop_offset_samples=0,
    window_size_samples=window_size,
    window_stride_samples=window_size,
    on_last_window="drop",
    preload=True,
)
metadata = windows.get_metadata()
assert (metadata.i_window_in_trial == 0).all(), "Expected one window per trial"
assert not metadata.duplicated(["subject", "session", "run", "i_start_in_trial"]).any()
y = metadata["target"].to_numpy(dtype=int)
groups = metadata["subject"].astype(str).to_numpy()
X = np.stack([window[0] for window in windows])
assert X.shape == (len(metadata), len(channel_names), window_size)
assert set(groups) == set(subjects)
assert np.isfinite(X).all()
print(pd.crosstab(groups, y, rownames=["subject"], colnames=["class"]))

class    0   1   2   3   4   5   6   7   8   9   10  11
subject
1        15  15  15  15  15  15  15  15  15  15  15  15
2        15  15  15  15  15  15  15  15  15  15  15  15
3        15  15  15  15  15  15  15  15  15  15  15  15

4. Extract spectral features from each window

SSVEP responses contain energy at the stimulus frequency. Use log spectral power around each stimulus frequency, retaining all eight posterior channels. This per-window transform learns nothing from other trials or subjects. The scaler below, in contrast, must be fitted only on training subjects. EEGDash’s shared spectral preprocessor computes a Welch PSD with one four-second Hann segment and 0.25 Hz bins. Each narrow band is centered on a documented stimulus frequency; these centers define the task, not trial-specific predictors. Retaining eight channels gives 12 × 8 = 96 features. The already processed SSVEP release does not need another EEGPrep cleaning pass or visual latency correction.

spectral_bands_power sums selected PSD bins. Multiplying by their 0.25 Hz spacing converts V²/Hz to approximate band power in V². The log compresses that scale; the StandardScaler still fits only on training participants.

bands = {
    f"hz_{name}": (float(name) - 0.125, float(name) + 0.125) for name in class_names
}
spectral = FeatureExtractor(
    {"power": partial(spectral_bands_power, bands=bands)},
    preprocessor=partial(
        spectral_preprocessor,
        fs=sfreq,
        nperseg=window_size,
        noverlap=0,
        f_min=8,
        f_max=16,
    ),
)
feature_table = extract_features(
    windows, {"spectral": spectral}, batch_size=64, n_jobs=1
).to_dataframe()
assert feature_table.shape == (len(y), len(class_names) * len(channel_names))
features = np.log(np.maximum(feature_table.to_numpy() * sfreq / window_size, 1e-30))
assert np.isfinite(features).all()

Extracting features:   0%|          | 0/3 [00:00<?, ?it/s]
Extracting features: 100%|██████████| 3/3 [00:00<00:00, 22.76it/s]
Extracting features: 100%|██████████| 3/3 [00:00<00:00, 22.71it/s]

5. Grow a nested training cohort; keep validation fixed

The seeded permutation picks subject order, never fabricates observations. Repeated inspection makes subject 3 a validation set, not a final test set. A larger study needs additional untouched test subjects and repeated orders. Both training subsets contain every frequency class. The one-participant subset contributes 180 trials; the two-participant subset contains those same trials plus the other person’s 180 trials. This nesting prevents a change in training membership from being mistaken for a change in size. Subject 3 supplies the same 180 validation trials at both points.

The feature extraction is fixed, but the scaler and classifier are refitted for each size. Reusing the larger fit would let information from the added participant enter the smaller-size result.

rng = np.random.default_rng(42)
training_order = rng.permutation(["1", "2"])
validation = np.flatnonzero(groups == "3")
rows = []
previous = set()
for n_subjects in (1, 2):
    train = np.flatnonzero(np.isin(groups, training_order[:n_subjects]))
    assert previous.issubset(set(train))
    assert set(groups[train]).isdisjoint(groups[validation])
    assert set(y[train]) == set(y[validation]) == set(mapping.values())
    previous = set(train)
    model = make_pipeline(StandardScaler(), LogisticRegression(max_iter=1000))
    model.fit(features[train], y[train])
    prediction = model.predict(features[validation])
    rows.append(
        dict(
            n_subjects=n_subjects,
            n_trials=len(train),
            balanced_accuracy=balanced_accuracy_score(y[validation], prediction),
        )
    )
results = pd.DataFrame(rows)
print("Training order:", training_order)
print(results.to_string(index=False))

Training order: ['2' '1']
 n_subjects  n_trials  balanced_accuracy
          1       180           0.144444
          2       360           0.255556

6. Display measured validation scores without extrapolation

plt.plot(results.n_subjects, results.balanced_accuracy, "o-")
plt.axhline(1 / len(mapping), color="black", linestyle="--", label="Chance (1/12)")
plt.xticks([1, 2])
plt.xlabel("Training participants")
plt.ylabel("Subject 3 validation balanced accuracy")
plt.ylim(0, 1)
plt.legend()
plt.show()

7. Read a small learning curve without over-interpreting it

A rising segment indicates improvement on this one validation participant for this particular order of training subjects. A flat or falling segment is equally valid: adding a participant changes both size and population composition. There are no error bars because only one order is evaluated.

With a larger cohort, repeat several nested training-subject orders while keeping validation identities fixed, and report the spread at each size. Reserve additional untouched people for the final chosen model. Do not extrapolate a sample-efficiency law from these two points.