nb-2-6-ica · EOG regression and ICA (L2.6)¶
Lesson L2.6 · Level 2 · Status draft — for expert review; uncertain points carry TODO(confirm).
- EOG regression first, because it is the transparent method and its limitation is the reason ICA exists.
- Rank, demonstrated rather than asserted: the same data decomposed with the correct number of components and with too many, and what "too many" actually produces.
- The two-pass strategy: fit on a 1 Hz-filtered copy with a recorded seed, apply the unmixing to the 0.1 Hz analysis data.
- Classification from the four views of the lesson — topography, time course, spectrum, ERP image — with
mne-icalabelas the second opinion where it installs (it needsonnxruntime, nottorch), and documented feature heuristics where it does not. - Before and after, so that the cost of each removed component is visible; then over-cleaning, by removing one more component that a classifier calls brain and measuring what that costs.
Data. ds-erpcore P3 (CC BY-SA 4.0 per data/directory.yaml, contested at source; open; per-subject downloadable): sub-001, sub-002, sub-003 — the same three subjects the site's w-ica-component-gallery decomposes, with the same fit settings (1–100 Hz, average reference, 256 Hz, extended Infomax, seed 20260917), so a component index here means the same component there. ~170 MB on an empty cache. TODO(confirm): the author mirrors the ERP CORE entry into the catalog registry and signs off the dataset page (§10.11 item 8).
Labels produced here are label_source: algorithmic (§4.5). No person has reviewed them.
In [1]:
# Setup: dependencies, the shared helpers, non-interactive plotting.
import importlib.util
import subprocess
import sys
import warnings
from pathlib import Path
# 1. Dependencies are pinned in notebooks/requirements.txt. Nothing is installed
# when the pinned stack is already present (local runs, CI); a fresh Colab or
# Binder kernel installs it once. On Colab, run from a clone of the repository
# so that notebooks/_shared/ is available (repository URL: TODO(confirm), spec
# section 13 item 3).
_needed = ('mne', 'scipy', 'matplotlib', 'pooch', 'pyprep')
_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", "pooch>=1.8"]
if "pyprep" in _missing:
_cmd += ["pyprep>=0.9"]
subprocess.check_call(_cmd)
# 2. Shared helpers (notebooks/_shared/helpers.py and helpers_l2.py), located
# relative to the working directory -- notebooks/<level>/ or notebooks/ --
# never through an absolute path.
_shared = next((d / "_shared" for d in (Path.cwd(), *Path.cwd().parents)
if (d / "_shared" / "helpers_l2.py").exists()), None)
if _shared is None:
raise FileNotFoundError("start the kernel in notebooks/L2/ (or notebooks/) so that _shared/helpers_l2.py is found")
sys.path.insert(0, str(_shared))
import helpers
import helpers_l2 as l2
# 3. Plotting: Jupyter's default inline backend renders static PNGs through Agg
# (no windows, nothing blocks); outside Jupyter the helpers select Agg. Every
# MNE figure is requested with show=False, and plt.show() renders each cell's
# figures in place.
import matplotlib.pyplot as plt
import numpy as np
import mne
# mne-icalabel is optional: it needs onnxruntime (a ~15 MB wheel), never torch. When it
# cannot be installed the helpers fall back to documented feature heuristics and say so.
if importlib.util.find_spec("mne_icalabel") is None:
try:
subprocess.check_call([sys.executable, "-m", "pip", "install", "-q", "mne-icalabel>=0.7", "onnxruntime"])
except Exception as _e:
print(f"mne-icalabel not installed ({_e}); documented feature heuristics will be used instead")
# Warnings are worth reading, so they are not silenced -- but their default format prints the
# absolute path of the file that raised them, which is nobody else's business and would put this
# machine's directory layout into the saved outputs. Only the class and the message are shown.
warnings.formatwarning = lambda message, category, *a, **k: f"{category.__name__}: {message}\n"
mne.viz.set_browser_backend("matplotlib", verbose=False)
mne.set_log_level("WARNING")
plt.rcParams["figure.dpi"] = 72
print(f"MNE {mne.__version__}; helpers imported from notebooks/_shared")
print("ERP CORE downloads (~56 MB per subject) are cached under data/downloads/erpcore/ in a repository "
"clone, otherwise under MNE's data directory; nothing is re-fetched.")
1. Load and prepare, in the canonical order¶
Everything up to (but not including) ICA: resample, detect bad channels, split into two branches at the filter step, interpolate both, average-reference both. The rank is carried forward from here.
The split is the part that is easy to get wrong. The analysis branch is filtered 0.1–30 Hz; the ICA branch is filtered 1–100 Hz from the same unfiltered recording, not from the analysis branch. Building the fit copy by high-passing data that has already been low-passed at 30 Hz leaves it with no content above 30 Hz — muscle and line-noise components then cannot be separated at all, and ICLabel is being asked to classify data far outside the 1–100 Hz regime it documents. Both branches then get the identical interpolation and the identical reference, so they differ only in their pass band.
In [2]:
import time
SUBJECTS = ["sub-001", "sub-002", "sub-003"]
RESAMPLE_HZ, L_FREQ, H_FREQ = 256.0, 0.1, 30.0
ICA_FIT_HZ = (1.0, 100.0)
TMIN, TMAX, BASELINE = -0.2, 0.8, (-0.2, 0.0)
CH, WINDOW = l2.P3_CHANNEL, l2.P3_WINDOW
REMOVE_CLASSES = ("eye", "muscle", "heart", "line noise", "channel noise")
MIN_PROBABILITY, MAX_REMOVED = 0.80, 6
prep = {}
t0 = time.time()
for sid in SUBJECTS:
raw, facts = l2.load_erpcore("P3", sid, resample_hz=RESAMPLE_HZ)
n_eeg = len(mne.pick_types(raw.info, eeg=True))
det = l2.detect_bad_channels(raw, highpass_hz=1.0, ransac=True, seed=l2.SEED)
dec = l2.select_bads(det, n_eeg, min_criteria=2, max_fraction=0.10)
ica_src = raw.copy().filter(ICA_FIT_HZ[0], ICA_FIT_HZ[1], picks=["eeg", "eog"], verbose=False) # ICA branch
raw.filter(L_FREQ, H_FREQ, picks=["eeg", "eog"], verbose=False) # analysis branch
for r in (raw, ica_src):
r.info["bads"] = list(dec["bads"])
if dec["bads"]:
r.interpolate_bads(reset_bads=True, verbose=False)
r.set_eeg_reference("average", verbose=False)
rank = l2.data_rank(n_eeg, n_interpolated=len(dec["bads"]), average_reference=True)
events, _, ev_info = l2.p3_events(raw)
prep[sid] = {"raw": raw, "ica_src": ica_src, "rank": rank, "bads": dec["bads"], "events": events,
"n_eeg": n_eeg}
print(f" {sid}: interpolated {dec['bads'] or 'none'} -> {rank['arithmetic']}; "
f"{ev_info['n_target']} target / {ev_info['n_standard']} standard stimuli; "
f"{raw.n_times} samples at {raw.info['sfreq']:g} Hz; analysis branch "
f"{L_FREQ:g}-{H_FREQ:g} Hz, ICA branch {ICA_FIT_HZ[0]:g}-{ICA_FIT_HZ[1]:g} Hz")
print(f"\nprepared {len(prep)} subjects in {time.time() - t0:.0f} s")
print("Removal policy, fixed before any component was looked at: remove a component whose class is one of "
f"{REMOVE_CLASSES} with probability >= {MIN_PROBABILITY:.2f}, at most {MAX_REMOVED} components.")
2. EOG regression: the transparent method, and its limitation¶
Model each EEG channel as brain activity plus a scaled copy of the EOG signal, estimate one propagation coefficient per channel, subtract. ERP CORE ships three EOG channels, so this is available here; many resting datasets have none.
Its virtue is that the whole model is one number per channel and can be plotted as a topography. Its limitation is that the EOG electrodes also record brain activity, so subtracting a scaled copy of them subtracts some frontal EEG too — which is measured below as the change the regression makes during blink-free stretches.
In [3]:
SID = SUBJECTS[0]
raw_r = prep[SID]["raw"]
eog_names = [raw_r.ch_names[i] for i in mne.pick_types(raw_r.info, eog=True)]
print(f"{SID}: EOG channels {eog_names}")
blinks = mne.preprocessing.find_eog_events(raw_r, ch_name=eog_names[-1], verbose=False)
print(f" {len(blinks)} blink events detected on {eog_names[-1]} "
f"({60 * len(blinks) / raw_r.times[-1]:.1f} per minute)")
model = mne.preprocessing.EOGRegression(picks="eeg", picks_artifact="eog", proj=False).fit(raw_r)
raw_reg = model.apply(raw_r.copy(), copy=False)
coefs = model.coef_ # n_eeg x n_eog
eeg_names = [raw_r.ch_names[i] for i in mne.pick_types(raw_r.info, eeg=True)]
print(f" regression coefficients: {coefs.shape[0]} EEG channels x {coefs.shape[1]} EOG channels, "
f"range {coefs.min():+.3f} to {coefs.max():+.3f} (dimensionless)")
for j, name in enumerate(eog_names):
big = np.argsort(-np.abs(coefs[:, j]))[:3]
print(f" {name:12s}: largest |coefficient| at " +
", ".join(f"{eeg_names[i]} {coefs[i, j]:+.2f}" for i in big))
# What does the regression change where there is no blink?
veog = raw_r.get_data(picks=[eog_names[-1]])[0]
sf = raw_r.info["sfreq"]
quiet = np.abs(veog - np.median(veog)) < 2 * np.std(veog)
for probe in ("Fp1", CH):
x0 = raw_r.get_data(picks=[probe])[0] * 1e6
x1 = raw_reg.get_data(picks=[probe])[0] * 1e6
print(f" {probe}: regression changes the trace by {np.sqrt(((x0 - x1) ** 2).mean()):5.2f} uV RMS overall, "
f"and by {np.sqrt(((x0 - x1)[quiet] ** 2).mean()):5.2f} uV RMS during blink-free samples "
f"({100 * quiet.mean():.0f} % of the recording) -- the second number is brain activity removed with "
"the artifact")
fig, axes = plt.subplots(1, 3, figsize=(14, 3.9))
for ax, j in zip(axes[:2], (0, len(eog_names) - 1)):
helpers.plot_topomap_values(raw_r, coefs[:, j], unit="coefficient (dimensionless)", ax=ax,
title=f"EOG regression weights for {eog_names[j]}")
t = raw_r.times
sl = (t >= blinks[2, 0] / sf - 1.5) & (t <= blinks[2, 0] / sf + 1.5) if len(blinks) > 2 else (t < 3)
axes[2].plot(t[sl], raw_r.get_data(picks=["Fp1"])[0][sl] * 1e6, color="k", lw=0.9, label="Fp1 before")
axes[2].plot(t[sl], raw_reg.get_data(picks=["Fp1"])[0][sl] * 1e6, color="tab:green", lw=0.9,
label="Fp1 after regression")
axes[2].plot(t[sl], raw_r.get_data(picks=[eog_names[-1]])[0][sl] * 1e6, color="tab:red", lw=0.7, alpha=0.6,
label=eog_names[-1])
axes[2].set(xlabel="Time (s)", ylabel="Amplitude (uV)", title="One blink, before and after regression (uV)")
axes[2].legend(fontsize=7); axes[2].grid(alpha=0.3)
fig.tight_layout()
plt.show() # render the static figure(s) of this cell inline
3. Rank: what "too many components" actually produces¶
ICA is a change of basis, so asking for more components than the data has independent dimensions is asking for something that does not exist. The claim is easy to test: fit the same data twice, once with the rank carried down from L2.2/L2.3 and once with the full channel count, and look at what the extra components are.
In [4]:
raw_rank, ica_src_rank = prep[SID]["raw"], prep[SID]["ica_src"]
rank_ok = prep[SID]["rank"]["rank"]
rank_too_many = prep[SID]["n_eeg"]
print(f"{SID}: {prep[SID]['rank']['arithmetic']}; the array still has {rank_too_many} rows.")
fits = {}
for label, n in ((f"correct rank ({rank_ok})", rank_ok), (f"channel count ({rank_too_many})", rank_too_many)):
t0 = time.time()
with warnings.catch_warnings():
warnings.simplefilter("ignore") # MNE warns when n_components exceeds the estimated rank
ica_x, raw_fit_x, log_x = l2.two_pass_ica(raw_rank, n_components=n, fit_on=ica_src_rank,
fit_l_freq=ICA_FIT_HZ[0], fit_h_freq=ICA_FIT_HZ[1],
seed=l2.SEED)
mix = ica_x.get_components()
corr = np.corrcoef(mix.T)
np.fill_diagonal(corr, 0.0)
evr = ica_x.pca_explained_variance_[:ica_x.n_components_]
fits[label] = {"ica": ica_x, "raw_fit": raw_fit_x, "n": ica_x.n_components_,
"max |r| between component maps": float(np.abs(corr).max()),
"pairs with |r| > 0.9": int((np.abs(np.triu(corr, 1)) > 0.9).sum()),
"smallest / largest PCA variance": float(evr.min() / evr.max()),
"fit seconds": round(time.time() - t0, 1)}
print(f" {label}: requested {n}, obtained {ica_x.n_components_} components in {fits[label]['fit seconds']:.0f} s")
print()
print(l2.fmt_table([dict(decomposition=k, **{kk: vv for kk, vv in v.items()
if kk not in ("ica", "raw_fit")}) for k, v in fits.items()],
["decomposition", "n", "max |r| between component maps", "pairs with |r| > 0.9",
"smallest / largest PCA variance", "fit seconds"], floatfmt="{:.4g}"))
print()
print("The ratio of the smallest to the largest PCA variance is the number to watch: when it approaches zero "
"the last components are being estimated from numerical noise. Asking for the channel count on "
"average-referenced, interpolated data is asking for exactly that (pf-interpolation-rank).")
ICA, RAW_FIT = fits[f"correct rank ({rank_ok})"]["ica"], fits[f"correct rank ({rank_ok})"]["raw_fit"]
4. The two-pass fit, on all three subjects¶
ica.fit sees a 1–100 Hz, average-referenced copy; the unmixing is applied to the 0.1–30 Hz analysis data. The fit filter is also what ICLabel documents as its training regime, so the classifier below is used inside the conditions it was trained for.
The data-quantity rule of thumb from the lesson is that the number of samples should be some multiple of the squared component count — TODO(confirm) which multiplier to quote; the ratio is printed so that the reader can apply whichever value their reference gives.
In [5]:
for sid in SUBJECTS:
p = prep[sid]
ica, raw_fit, log = l2.two_pass_ica(p["raw"], n_components=p["rank"]["rank"], fit_on=p["ica_src"],
fit_l_freq=ICA_FIT_HZ[0], fit_h_freq=ICA_FIT_HZ[1], seed=l2.SEED)
cls = l2.classify_components(raw_fit, ica)
p.update(ica=ica, raw_fit=raw_fit, ica_log=log, cls=cls)
print(f" {sid}: {log['n_components']} components, {log['method']} (extended), seed {log['seed']}, "
f"{log['duration_s']:.0f} s; {log['n_samples_fit']} samples = "
f"{log['samples_per_squared_channel']:.0f} x n_components^2; fit copy "
f"{log['fit_filter_hz'][0]:g}-{log['fit_lowpass_hz']:g} Hz ({log['fit_copy']})")
counts = {c: cls["labels"].count(c) for c in l2.ICLABEL_CLASSES if cls["labels"].count(c)}
print(f" classifier: {cls['tool']}")
print(f" label counts: {counts}")
5. The four views¶
One component per class that the classifier is most confident about, shown the way the lesson says a component must be read: topography, time course, spectrum and ERP image together. A label that only one view supports is a label to distrust.
In [6]:
from scipy import signal as sp_signal
p = prep[SID]
ica, raw_fit, cls = p["ica"], p["raw_fit"], p["cls"]
sources = ica.get_sources(raw_fit).get_data()
mixing = ica.get_components()
sf = raw_fit.info["sfreq"]
ep_src = mne.Epochs(ica.get_sources(raw_fit), p["events"], dict(l2.P3_EVENT_ID), tmin=TMIN, tmax=TMAX,
baseline=BASELINE, preload=True, verbose=False)
# One exemplar per class present, the most confident of each.
show = []
for klass in l2.ICLABEL_CLASSES:
idx = [i for i, lab in enumerate(cls["labels"]) if lab == klass]
if idx:
show.append(max(idx, key=lambda i: cls["probabilities"][i]))
show = show[:4]
info_eeg = mne.pick_info(raw_fit.info, mne.pick_types(raw_fit.info, eeg=True))
nper = int(min(4 * sf, sources.shape[1]))
freqs, psd = sp_signal.welch(sources, fs=sf, nperseg=nper, noverlap=nper // 2, axis=-1)
fig, axes = plt.subplots(len(show), 4, figsize=(15, 3.1 * len(show)))
axes = np.atleast_2d(axes)
for r, comp in enumerate(show):
im, _ = mne.viz.plot_topomap(mixing[:, comp], info_eeg, axes=axes[r, 0], show=False, contours=4)
axes[r, 0].set_title(f"IC{comp}: {cls['labels'][comp]} (p {cls['probabilities'][comp]:.2f})", fontsize=9)
t0_s = 200.0
sl = slice(int(t0_s * sf), int((t0_s + 10) * sf))
axes[r, 1].plot(raw_fit.times[sl], sources[comp][sl], "k", lw=0.6)
axes[r, 1].set(xlabel="Time (s)", ylabel="Component amplitude (a.u.)",
title=f"IC{comp} time course, 10 s (arbitrary units)")
axes[r, 1].grid(alpha=0.3)
axes[r, 2].semilogy(freqs, psd[comp], "k", lw=0.8)
axes[r, 2].set(xlim=(0, 60), xlabel="Frequency (Hz)", ylabel="Power (a.u.^2/Hz)",
title=f"IC{comp} spectrum (arbitrary units)")
axes[r, 2].grid(alpha=0.3, which="both")
img = ep_src.get_data(picks=[comp])[:, 0, :]
vmax = np.percentile(np.abs(img), 98)
axes[r, 3].imshow(img, aspect="auto", cmap="RdBu_r", vmin=-vmax, vmax=vmax, origin="lower",
extent=[ep_src.times[0] * 1000, ep_src.times[-1] * 1000, 0, len(img)])
axes[r, 3].axvline(0, color="k", lw=0.6)
axes[r, 3].set(xlabel="Time from stimulus (ms)", ylabel="Trial",
title=f"IC{comp} ERP image (arbitrary units)")
fig.suptitle(f"{SID}: four views of one component per class (ICA on a "
f"{ICA_FIT_HZ[0]:g}-{ICA_FIT_HZ[1]:g} Hz average-referenced copy, seed {l2.SEED})", y=1.005)
fig.tight_layout()
plt.show() # render the static figure(s) of this cell inline
print("Evidence recorded for each shown component:")
for comp in show:
w = np.abs(mixing[:, comp])
band = (freqs >= 20) & (freqs <= 45)
slope = float(np.polyfit(np.log10(freqs[band]), np.log10(psd[comp][band] + 1e-30), 1)[0])
eog = raw_fit.get_data(picks="eog")
r_eog = max(abs(float(np.corrcoef(sources[comp], e)[0, 1])) for e in eog)
print(f" IC{comp:2d} {cls['labels'][comp]:14s} p {cls['probabilities'][comp]:.2f} | "
f"map: {100 * w.max() / w.sum():.0f} % of the weight on {raw_fit.ch_names[int(np.argmax(w))]} | "
f"spectrum 20-45 Hz slope {slope:+.2f} | |r| with EOG {r_eog:.2f} | {cls['evidence'][comp]}")
6. Before and after¶
The policy was fixed in section 1. Applying it: the components it selects, the ERP before and after, and the change in the measured amplitude. The trial set is identical in both, so nothing here is a rejection effect.
In [7]:
summary = []
for sid in SUBJECTS:
p = prep[sid]
cls, ica = p["cls"], p["ica"]
exclude = [i for i, (lab, prob) in enumerate(zip(cls["labels"], cls["probabilities"]))
if lab in REMOVE_CLASSES and prob >= MIN_PROBABILITY]
exclude = sorted(sorted(exclude, key=lambda i: -cls["probabilities"][i])[:MAX_REMOVED])
ica.exclude = exclude
p["exclude"] = exclude
before = l2.epochs_p3(p["raw"], p["events"], tmin=TMIN, tmax=TMAX, baseline=BASELINE)
after = l2.epochs_p3(ica.apply(p["raw"].copy(), verbose=False), p["events"], tmin=TMIN, tmax=TMAX,
baseline=BASELINE)
keep = l2.epoch_ptp_uv(after) <= 150.0 # one trial set, decided on the cleaned data
p["before"], p["after"], p["keep"] = before, after, keep
d_b, d_a = l2.difference_wave(before[keep]), l2.difference_wave(after[keep])
p["d_before"], p["d_after"] = d_b, d_a
snr_b, snr_a = l2.erp_snr(d_b, CH, WINDOW, BASELINE), l2.erp_snr(d_a, CH, WINDOW, BASELINE)
summary.append({"subject": sid, "rank": p["rank"]["rank"], "removed": len(exclude),
"components": exclude or ["-"],
"classes": [f"{cls['labels'][i]}" for i in exclude] or ["-"],
"trials kept": int(keep.sum()),
f"{CH} mean before (uV)": snr_b["signal_uv"], f"{CH} mean after (uV)": snr_a["signal_uv"],
"baseline noise before (uV)": snr_b["noise_uv"], "baseline noise after (uV)": snr_a["noise_uv"],
"SNR before": snr_b["snr"], "SNR after": snr_a["snr"],
"rank after cleaning": p["rank"]["rank"] - len(exclude)})
print(l2.fmt_table(summary, ["subject", "rank", "removed", "components", "classes", "trials kept",
f"{CH} mean before (uV)", f"{CH} mean after (uV)", "baseline noise before (uV)",
"baseline noise after (uV)", "SNR before", "SNR after", "rank after cleaning"],
floatfmt="{:+.2f}"))
In [8]:
fig, axes = plt.subplots(1, 3, figsize=(15, 4.2), sharey=True)
for ax, sid in zip(axes, SUBJECTS):
p = prep[sid]
for d, color, label in ((p["d_before"], "0.45", "before ICA"), (p["d_after"], "tab:green", "after ICA")):
i = d.ch_names.index(CH)
ax.plot(d.times * 1000, d.data[i] * 1e6, lw=1.5, color=color, label=f"{label} (n = {d.nave})")
ax.axvspan(WINDOW[0] * 1000, WINDOW[1] * 1000, color="tab:orange", alpha=0.18)
ax.axhline(0, color="gray", lw=0.6); ax.axvline(0, color="gray", lw=0.6)
ax.set(xlabel="Time from stimulus (ms)",
title=f"{sid}: {len(p['exclude'])} removed {p['exclude'] or ''}")
ax.grid(alpha=0.3); ax.legend(fontsize=8)
axes[0].set_ylabel("Amplitude (uV)")
fig.suptitle(f"Target minus standard at {CH}, before and after ICA cleaning (uV, positive up, identical trials)",
y=1.02)
fig.tight_layout()
plt.show() # render the static figure(s) of this cell inline
7. Over-cleaning¶
The error in the other direction. The policy above removes components a classifier calls eye, muscle, heart, line noise or channel noise. Suppose the policy were sloppier — "remove anything that does not obviously look like a signal" — and one more component went with them: the highest-confidence brain component. The cost is measured the same way as the benefit.
In [9]:
p = prep[SID]
cls, ica = p["cls"], p["ica"]
brain = [i for i, lab in enumerate(cls["labels"]) if lab == "brain"]
# Which brain component would cost the most if it went? Removing each one in turn answers that exactly,
# and the answer is the component a sloppier policy is most dangerous to.
base_amp = None
cost = {}
for i in brain:
ica.exclude = sorted(p["exclude"] + [i])
ep_i = l2.epochs_p3(ica.apply(p["raw"].copy(), verbose=False), p["events"], tmin=TMIN, tmax=TMAX,
baseline=BASELINE)[p["keep"]]
cost[i] = l2.mean_amplitude(l2.difference_wave(ep_i), CH, WINDOW)
ica.exclude = list(p["exclude"])
base_amp = l2.mean_amplitude(p["d_after"], CH, WINDOW)
extra = min(cost, key=lambda i: cost[i]) if cost else None
if cost:
print(f"cost of removing each of the {len(brain)} components the classifier calls brain, on top of the "
f"stated policy ({CH} mean amplitude, {base_amp:+.2f} uV before any of them goes):")
for i in sorted(cost, key=lambda i: cost[i])[:5]:
print(f" IC{i:<3d} p {cls['probabilities'][i]:.2f} -> {cost[i]:+.2f} uV "
f"({cost[i] - base_amp:+.2f} uV)")
print(f" the most expensive is IC{extra} at {cost[extra] - base_amp:+.2f} uV; that is the component the "
"sloppier policy below removes.")
mixing = ica.get_components()
sources = ica.get_sources(p["raw_fit"]).get_data()
nper = int(min(4 * p["raw_fit"].info["sfreq"], sources.shape[1]))
freqs2, psd2 = sp_signal.welch(sources, fs=p["raw_fit"].info["sfreq"], nperseg=nper, noverlap=nper // 2, axis=-1)
rows = []
policies = {
"policy as stated": p["exclude"],
f"policy + one brain component (IC{extra})": sorted(p["exclude"] + ([extra] if extra is not None else [])),
"nothing removed": [],
}
for name, ex in policies.items():
ica.exclude = list(ex)
cleaned = ica.apply(p["raw"].copy(), verbose=False)
ep = l2.epochs_p3(cleaned, p["events"], tmin=TMIN, tmax=TMAX, baseline=BASELINE)[p["keep"]]
d = l2.difference_wave(ep)
snr = l2.erp_snr(d, CH, WINDOW, BASELINE)
rows.append({"removed": name, "components": ex or ["-"], f"{CH} mean (uV)": snr["signal_uv"],
"baseline noise (uV)": snr["noise_uv"], "SNR": snr["snr"],
"rank after cleaning": p["rank"]["rank"] - len(ex)})
ica.exclude = list(p["exclude"]) # restore the stated policy
print(l2.fmt_table(rows, list(rows[0]), floatfmt="{:+.2f}"))
if extra is not None:
band = (freqs2 >= 8) & (freqs2 <= 12)
alpha_share = float(psd2[extra][band].sum() / psd2[extra].sum())
w = np.abs(mixing[:, extra])
print(f"\nIC{extra} is the component the sloppier policy removes: classifier label "
f"{cls['labels'][extra]!r} with probability {cls['probabilities'][extra]:.2f}, "
f"{100 * alpha_share:.0f} % of its power between 8 and 12 Hz, map weight concentrated at "
f"{p['raw_fit'].ch_names[int(np.argmax(w))]} ({100 * w.max() / w.sum():.0f} % of the total). "
"It is the brain component whose removal costs the most, found by removing each one in turn "
"rather than by guessing.")
lost = rows[1][f"{CH} mean (uV)"] - rows[0][f"{CH} mean (uV)"]
print(f"Removing it changes the measured P3 by {lost:+.2f} uV "
f"({100 * lost / abs(rows[0][f'{CH} mean (uV)']):+.0f} %) and drops the rank by one more. "
"The symptom list of pf-overcleaning-ica starts with exactly this: the effect shrinks after cleaning.")
fig, axes = plt.subplots(1, 2, figsize=(12, 4.2))
ica.exclude = list(p["exclude"])
for name, ex, color in (("nothing removed", [], "0.45"),
("policy as stated", p["exclude"], "tab:green"),
(f"policy + IC{extra}", policies[f"policy + one brain component (IC{extra})"], "tab:red")):
ica.exclude = list(ex)
ep = l2.epochs_p3(ica.apply(p["raw"].copy(), verbose=False), p["events"], tmin=TMIN, tmax=TMAX,
baseline=BASELINE)[p["keep"]]
d = l2.difference_wave(ep)
axes[0].plot(d.times * 1000, d.data[d.ch_names.index(CH)] * 1e6, lw=1.5, color=color, label=name)
ica.exclude = list(p["exclude"])
axes[0].axvspan(WINDOW[0] * 1000, WINDOW[1] * 1000, color="tab:orange", alpha=0.18)
axes[0].axhline(0, color="gray", lw=0.6); axes[0].axvline(0, color="gray", lw=0.6)
axes[0].set(xlabel="Time from stimulus (ms)", ylabel="Amplitude (uV)",
title=f"{SID}: what one extra removed component costs ({CH}, uV, positive up)")
axes[0].grid(alpha=0.3); axes[0].legend(fontsize=8)
if extra is not None:
im, _ = mne.viz.plot_topomap(mixing[:, extra], info_eeg, axes=axes[1], show=False, contours=4)
axes[1].set_title(f"IC{extra}: {cls['labels'][extra]} (p {cls['probabilities'][extra]:.2f}), "
f"{100 * alpha_share:.0f} % of power 8-12 Hz", fontsize=9)
fig.tight_layout()
plt.show() # render the static figure(s) of this cell inline
8. The numbers¶
In [10]:
print("nb-2-6-ica -- L2.6 numbers (draft; TODO(confirm) at author review)")
print(f"Data: ds-erpcore P3, {', '.join(SUBJECTS)} (CC BY-SA 4.0 per data/directory.yaml, contested at source, open, per-subject "
f"downloadable). 30 EEG + 3 EOG, 1024 Hz resampled to {RESAMPLE_HZ:g} Hz.")
print(f"Pipeline: pyprep bad channels (>=2 criteria, cap 10 %) -> FIR zero-phase {L_FREQ:g}-{H_FREQ:g} Hz on "
f"the continuous data -> interpolate -> average reference -> ICA -> epoch {TMIN:g}..{TMAX:g} s, baseline "
f"{BASELINE[0]:g}..{BASELINE[1]:g} s.")
print(f"ICA: two-pass -- fitted on a {ICA_FIT_HZ[0]:g}-{ICA_FIT_HZ[1]:g} Hz copy of the same data, extended "
f"Infomax, n_components = the carried rank, random_state {l2.SEED}; the unmixing is applied to the "
f"{L_FREQ:g}-{H_FREQ:g} Hz analysis data.")
print(f"Classifier: {prep[SID]['cls']['tool']}. Labels are label_source: algorithmic (spec 4.5) -- no person "
"has reviewed them.")
print(f"Policy: remove classes {REMOVE_CLASSES} with probability >= {MIN_PROBABILITY:.2f}, at most "
f"{MAX_REMOVED} components.")
print()
print("Components removed, with their classes and probabilities:")
for sid in SUBJECTS:
p = prep[sid]
cls = p["cls"]
print(f" {sid} (rank {p['rank']['rank']}, {len(cls['labels'])} components):")
if not p["exclude"]:
print(" nothing met the policy")
for i in p["exclude"]:
print(f" IC{i:<3d} {cls['labels'][i]:14s} p {cls['probabilities'][i]:.3f} {cls['evidence'][i]}")
counts = {c: cls["labels"].count(c) for c in l2.ICLABEL_CLASSES if cls["labels"].count(c)}
print(f" all labels: {counts}; rank after cleaning {p['rank']['rank'] - len(p['exclude'])}")
print()
print(f"Effect of the cleaning on the measured P3 ({CH}, {WINDOW[0] * 1000:.0f}-"
f"{WINDOW[1] * 1000:.0f} ms, identical trials before and after):")
print(l2.fmt_table(summary, ["subject", "removed", "trials kept", f"{CH} mean before (uV)",
f"{CH} mean after (uV)", "baseline noise before (uV)",
"baseline noise after (uV)", "SNR before", "SNR after"], floatfmt="{:+.2f}"))
print()
print(f"Over-cleaning ({SID}); every component the classifier calls brain was removed in turn to find the "
f"most expensive one (IC{extra}):")
print(l2.fmt_table(rows, ["removed", "components", f"{CH} mean (uV)", "baseline noise (uV)", "SNR",
"rank after cleaning"], floatfmt="{:+.2f}"))
print()
_ok_fit, _bad_fit = fits[f"correct rank ({rank_ok})"], fits[f"channel count ({rank_too_many})"]
print(f"Rank check ({SID}): asking for {rank_too_many} components instead of {rank_ok} drives the ratio of the "
f"smallest to the largest PCA variance from {_ok_fit['smallest / largest PCA variance']:.3g} to "
f"{_bad_fit['smallest / largest PCA variance']:.3g} -- the last components are being estimated from "
f"numerical noise -- and raises the largest correlation between two component topographies from "
f"{_ok_fit['max |r| between component maps']:.2f} to {_bad_fit['max |r| between component maps']:.2f}. "
f"The fit also takes {_bad_fit['fit seconds'] / max(_ok_fit['fit seconds'], 0.1):.1f} times as long. "
f"Outright near-duplicate pairs (|r| > 0.9) did not appear here, which is worth saying: the symptom "
"list of pf-interpolation-rank is a list of things that can happen, not a checklist that always fires.")
print()
print("L2.6's exercises are a drill in w-ica-component-gallery (>= 85 % agreement on 20 components) and a "
"free response on one ambiguous component; this notebook supplies the decompositions and the evidence "
"rather than a numeric key.")