nb-2-7-asr-vs-ica · ASR, SSP and alternatives (L2.7)¶
Lesson L2.7 · Level 2 · Status draft — for expert review; uncertain points carry TODO(confirm).
One recording, cleaned four ways on identical trials, so that nothing in the comparison is a rejection effect:
- uncorrected — the baseline everything is measured against;
- ICA, exactly as fitted in
nb-2-6-ica(two-pass, extended Infomax, recorded seed, components chosen by a stated policy); - ASR, artifact subspace reconstruction;
- SSP, one EOG projector — the method whose cost is countable, because each projector removes exactly one dimension.
For each: the continuous traces, the ERP, the spectra, the ERP signal-to-noise ratio and the rank of the output. The last two carry the lesson — ASR's output looks full rank and there is no list of what it removed.
One asymmetry is unavoidable and is itself a finding. ASR's calibration compares each window's principal-component variances against a clean reference. On data high-passed at 0.1 Hz, slow drift dominates that covariance, and section 3 measures what happens: ASR replaces essentially the whole recording. Every implementation therefore expects a ~1 Hz high-passed input. So the comparison below runs ICA and SSP on the 0.1–30 Hz ERP band and ASR on a 1–30 Hz band, and reports the matched uncorrected baseline for each band so that the filter's own contribution is visible and is not mistaken for the cleaning's.
ASR is not part of MNE. This notebook uses asrpy (a 20 kB pure-Python wheel whose only dependencies are MNE, NumPy and SciPy), installs it in the first cell, and records its version and cutoff. It also runs a documented NumPy implementation of the published calibration/cutoff idea as a cross-check, because the two disagree on this recording by more than they should — which is reported rather than hidden, and marked TODO(confirm).
Data. ds-erpcore P3 (CC BY-SA 4.0 per data/directory.yaml, contested at source; open; per-subject downloadable): sub-001, 30 EEG + 3 EOG, 1024 Hz, CMS reference, 60 Hz mains, no software filters (~56 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).
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
# ASR is not in MNE. asrpy is a 20 kB pure-Python wheel (dependencies: mne, numpy, scipy) --
# no torch, no compiled extension. If it cannot be installed, only the notebook's own
# documented implementation runs, and every table below says so.
ASR_BACKEND = None
if importlib.util.find_spec("asrpy") is None:
try:
subprocess.check_call([sys.executable, "-m", "pip", "install", "-q", "asrpy"])
except Exception as _e:
print(f"asrpy could not be installed ({_e}); only the documented implementation will run")
try:
import asrpy
ASR_BACKEND = f"asrpy {getattr(asrpy, '__version__', '0.0.8')}"
except Exception as _e:
asrpy = None
ASR_BACKEND = f"unavailable ({type(_e).__name__})"
# 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. Prepare the recording once, in two bands¶
The canonical order through re-referencing, split at the filter step into three branches that differ only in their pass band: the 0.1–30 Hz ERP band, a 1–30 Hz band for ASR, and a 1–100 Hz branch for the ICA fit. Bad channels and the reference are identical in all three.
In [2]:
import time
SID = "sub-001"
RESAMPLE_HZ = 256.0
ERP_BAND, ASR_BAND, ICA_FIT_HZ = (0.1, 30.0), (1.0, 30.0), (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
ASR_CUTOFF = 20.0 # standard deviations; swept in section 6
SELECT_PTP_UV = 150.0
raw0, facts = l2.load_erpcore("P3", SID, resample_hz=RESAMPLE_HZ, verbose=True)
n_eeg = len(mne.pick_types(raw0.info, eeg=True))
det = l2.detect_bad_channels(raw0, highpass_hz=1.0, ransac=True, seed=l2.SEED)
dec = l2.select_bads(det, n_eeg, min_criteria=2, max_fraction=0.10)
branches = {}
for name, band in (("erp", ERP_BAND), ("asr", ASR_BAND), ("ica_fit", ICA_FIT_HZ)):
r = raw0.copy().filter(band[0], band[1], picks=["eeg", "eog"], verbose=False)
r.info["bads"] = list(dec["bads"])
if dec["bads"]:
r.interpolate_bads(reset_bads=True, verbose=False)
r.set_eeg_reference("average", verbose=False)
branches[name] = r
RANK = l2.data_rank(n_eeg, n_interpolated=len(dec["bads"]), average_reference=True)
events, _, ev_info = l2.p3_events(branches["erp"])
print(f"{SID}: interpolated {dec['bads'] or 'none'} -> {RANK['arithmetic']}")
print(f" {ev_info['n_target']} target / {ev_info['n_standard']} standard stimuli; "
f"{branches['erp'].times[-1]:.0f} s at {branches['erp'].info['sfreq']:g} Hz")
print(f" branches: ERP {ERP_BAND[0]:g}-{ERP_BAND[1]:g} Hz, ASR {ASR_BAND[0]:g}-{ASR_BAND[1]:g} Hz, "
f"ICA fit {ICA_FIT_HZ[0]:g}-{ICA_FIT_HZ[1]:g} Hz -- same bads, same reference")
print(f" signal RMS over the 30 EEG channels: "
+ ", ".join(f"{k} {np.sqrt((v.get_data(picks='eeg') ** 2).mean()) * 1e6:.1f} uV" for k, v in branches.items()))
print(f"ASR backend: {ASR_BACKEND}")
2. ICA (the reference method)¶
Identical settings to nb-2-6-ica, so the two notebooks describe the same decomposition: fitted on the 1–100 Hz branch, applied to the 0.1–30 Hz ERP band.
In [3]:
ica, raw_for_ica, ica_log = l2.two_pass_ica(branches["erp"], n_components=RANK["rank"],
fit_on=branches["ica_fit"], fit_l_freq=ICA_FIT_HZ[0],
fit_h_freq=ICA_FIT_HZ[1], seed=l2.SEED)
cls = l2.classify_components(raw_for_ica, ica)
exclude = [i for i, (lab, p) in enumerate(zip(cls["labels"], cls["probabilities"]))
if lab in REMOVE_CLASSES and p >= MIN_PROBABILITY]
exclude = sorted(sorted(exclude, key=lambda i: -cls["probabilities"][i])[:MAX_REMOVED])
ica.exclude = exclude
raw_ica = ica.apply(branches["erp"].copy(), verbose=False)
print(f"ICA: {ica_log['n_components']} components, {ica_log['method']} (extended), seed {ica_log['seed']}, "
f"fitted in {ica_log['duration_s']:.0f} s on the {ICA_FIT_HZ[0]:g}-{ica_log['fit_lowpass_hz']:g} Hz branch")
print(f" classifier: {cls['tool']}")
print(f" removed {len(exclude)} of {ica.n_components_}: "
+ (", ".join(f"IC{i} {cls['labels'][i]} (p {cls['probabilities'][i]:.2f})" for i in exclude) or "nothing"))
print(f" rank after cleaning: {RANK['rank'] - len(exclude)} (one per removed component) -- "
"and the list of what went is the point")
3. ASR, and why the pass band decides what it does¶
Two stages. Calibration: find a stretch of relatively clean data in the same recording, compute its channel covariance and its principal components, and record how much variance each normally carries. Processing: slide a short window (about half a second) over the recording; in each window, any principal component whose variance exceeds the calibration value by more than a cutoff (in standard deviations) is declared an artifact subspace, and the window is reconstructed from the remaining components rather than having the artifact deleted.
Everything in that description is a statement about the covariance, and on 0.1 Hz-filtered EEG the covariance is dominated by drift. The first table measures the consequence: the same tool, the same cutoff, two pass bands.
In [4]:
def asr_documented(data_uv, sfreq, *, cutoff=ASR_CUTOFF, win_len=0.5, calib_seconds=60.0):
"""The published calibration/cutoff idea written out in NumPy, as a cross-check on the package.
Calibration: the `calib_seconds` of the recording with the lowest total variance are the clean
reference; its covariance gives the principal directions and the variance each normally carries.
Processing: in each `win_len` window, any direction whose variance exceeds `cutoff` standard-deviation
equivalents of its calibration value is dropped and the window is reconstructed from the rest.
This is a teaching implementation, not a re-implementation of any specific tool: it has no robust
(RANSAC-style) calibration search and no blending between windows. TODO(confirm) against a reference
implementation before quoting its numbers as ASR's.
"""
x = np.asarray(data_uv, float)
n_ch, n_t = x.shape
w = int(round(win_len * sfreq))
n_win = n_t // w
trimmed = x[:, : n_win * w].reshape(n_ch, n_win, w)
var_per_window = trimmed.var(axis=2).sum(axis=0)
n_calib = max(1, int(round(calib_seconds / win_len)))
calib = np.argsort(var_per_window)[:n_calib]
ref = trimmed[:, calib, :].reshape(n_ch, -1)
evals, evecs = np.linalg.eigh(np.cov(ref))
thresh = evals * (cutoff ** 2)
out, n_removed, n_windows_touched = trimmed.copy(), 0, 0
for i in range(n_win):
seg = trimmed[:, i, :]
proj = evecs.T @ (seg - seg.mean(axis=1, keepdims=True))
bad = proj.var(axis=1) > thresh
if bad.any():
n_removed += int(bad.sum())
n_windows_touched += 1
keep = evecs[:, ~bad]
out[:, i, :] = keep @ (keep.T @ seg)
rebuilt = np.concatenate([out.reshape(n_ch, -1), x[:, n_win * w:]], axis=1)
return rebuilt, {"tool": "helpers-local NumPy implementation of the documented idea (TODO(confirm))",
"cutoff_sd": cutoff, "win_len_s": win_len, "n_windows": n_win,
"calibration_windows": int(n_calib), "windows_touched": n_windows_touched,
"component_windows_removed": n_removed}
def run_asr(raw_in, *, cutoff=ASR_CUTOFF, backend="asrpy"):
"""ASR on one Raw; returns (cleaned Raw, log)."""
t0 = time.time()
r = raw_in.copy().pick("eeg")
if backend == "asrpy" and asrpy is not None:
a = asrpy.ASR(sfreq=r.info["sfreq"], cutoff=cutoff)
a.fit(r.copy())
out = a.transform(r.copy())
log = {"tool": ASR_BACKEND, "cutoff_sd": cutoff, "win_len_s": 0.5,
"calibration": "the package's own clean-window search on this recording"}
else:
rebuilt, log = asr_documented(r.get_data() * 1e6, r.info["sfreq"], cutoff=cutoff)
out = mne.io.RawArray(rebuilt * 1e-6, r.info.copy(), verbose=False)
log["duration_s"] = round(time.time() - t0, 1)
x0, x1 = r.get_data() * 1e6, out.get_data() * 1e6
log["rms_before_uv"] = float(np.sqrt((x0 ** 2).mean()))
log["rms_after_uv"] = float(np.sqrt((x1 ** 2).mean()))
log["changed_uv_rms"] = float(np.sqrt(((x0 - x1) ** 2).mean()))
log["percent_of_signal_changed"] = 100 * log["changed_uv_rms"] / log["rms_before_uv"]
return out, log
rows = []
for band_name, src in (("0.1-30 Hz (the ERP band)", branches["erp"]), ("1-30 Hz (high-passed first)", branches["asr"])):
for backend in (["asrpy", "documented"] if asrpy is not None else ["documented"]):
_, lg = run_asr(src, cutoff=ASR_CUTOFF, backend=backend)
rows.append({"input band": band_name, "implementation": lg["tool"][:44],
"signal RMS before (uV)": lg["rms_before_uv"], "after (uV)": lg["rms_after_uv"],
"replaced (uV RMS)": lg["changed_uv_rms"],
"% of the signal replaced": lg["percent_of_signal_changed"], "seconds": lg["duration_s"]})
print(f"ASR at cutoff {ASR_CUTOFF:g} SD, two input bands:\n")
print(l2.fmt_table(rows, list(rows[0]), floatfmt="{:.2f}"))
print()
print("Read the last column. Given 0.1 Hz-filtered data, ASR replaces almost the entire recording: the drift "
"subspace carries most of the variance, so almost every window looks like an outlier against a clean "
"reference that happens to hold less of it. High-passing first is not a convenience, it is a "
"precondition -- and this is exactly the failure mode the lesson warns about, because the output still "
"looks like continuous, full-rank EEG.")
print()
print("The two implementations also disagree on how much they take at the same nominal cutoff. That "
"disagreement is reported, not reconciled: TODO(confirm) which calibration a reference implementation "
"performs before treating either number as 'what ASR does'.")
In [5]:
raw_asr, asr_log = run_asr(branches["asr"], cutoff=ASR_CUTOFF, backend="asrpy" if asrpy is not None else "documented")
raw_asr_doc, asr_doc_log = run_asr(branches["asr"], cutoff=ASR_CUTOFF, backend="documented")
print(f"ASR used for the comparison below: {asr_log['tool']}, cutoff {ASR_CUTOFF:g} SD, "
f"window {asr_log.get('win_len_s', 0.5):g} s, input {ASR_BAND[0]:g}-{ASR_BAND[1]:g} Hz, "
f"{asr_log['duration_s']:.0f} s")
for k, v in asr_log.items():
if k not in ("tool", "duration_s"):
print(f" {k}: {v if not isinstance(v, float) else round(v, 2)}")
x_after = raw_asr.get_data() * 1e6
rank_asr = int(np.linalg.matrix_rank(x_after, tol=1e-6 * np.abs(x_after).max()))
print(f" numerical rank of the ASR output: {rank_asr}, against a carried rank of {RANK['rank']} before "
"cleaning. ASR reconstructs rather than deletes, so nothing in the output announces what it took.")
4. SSP, where the cost is countable¶
In [6]:
projs, _ = mne.preprocessing.compute_proj_eog(branches["erp"].copy(), n_eeg=1, reject=None, no_proj=True,
verbose=False)
raw_ssp = branches["erp"].copy().add_proj(projs, verbose=False).apply_proj(verbose=False)
print(f"SSP: {len(projs)} EOG projector(s) -- {[p['desc'] for p in projs]}")
print(f" each projector removes exactly one dimension: "
f"{l2.data_rank(n_eeg, n_interpolated=len(dec['bads']), average_reference=True, n_projectors=len(projs))['arithmetic']}")
x_ssp = raw_ssp.get_data(picks="eeg") * 1e6
print(f" numerical rank of the SSP output: {int(np.linalg.matrix_rank(x_ssp, tol=1e-6 * np.abs(x_ssp).max()))}")
print(" the projector's topography can be plotted and argued about, and the count is in the file. That "
"explicitness is SSP's main virtue over ASR.")
fig, ax = plt.subplots(figsize=(3.6, 3.6))
mne.viz.plot_projs_topomap(projs, info=mne.pick_info(branches["erp"].info,
mne.pick_types(branches["erp"].info, eeg=True)),
axes=ax, show=False)
ax.set_title("SSP: the one EOG projector's topography (arbitrary units)", fontsize=9)
fig.tight_layout()
plt.show() # render the static figure(s) of this cell inline
5. The comparison, on identical trials¶
The trial set is chosen once, on the uncorrected 0.1–30 Hz epochs, with a fixed 150 µV peak-to-peak criterion, and the same trials are then averaged under every method. Otherwise a method that removes an artifact would also change which trials survive, and the comparison would confound cleaning with rejection.
Each row names its input band, and each band's uncorrected row is in the table, so the reader can separate what the filter did from what the cleaning did.
The SNR is |mean amplitude in the a-priori window| / SD of the averaged pre-stimulus baseline at Pz, on the target-minus-standard difference wave — both terms from the same waveform, so the ratio is dimensionless.
In [7]:
methods = {
"uncorrected (0.1-30 Hz)": (branches["erp"], "-"),
"ICA (0.1-30 Hz)": (raw_ica, f"{len(exclude)} named components"),
"SSP (0.1-30 Hz)": (raw_ssp, f"{len(projs)} fixed projector(s)"),
"uncorrected (1-30 Hz)": (branches["asr"], "-"),
f"ASR, {asr_log['tool'].split()[0]} (1-30 Hz)": (raw_asr, "per-window subspaces, not listed"),
"ASR, documented implementation (1-30 Hz)": (raw_asr_doc, "per-window subspaces, not listed"),
}
ep_ref = l2.epochs_p3(branches["erp"], events, tmin=TMIN, tmax=TMAX, baseline=BASELINE)
KEEP = l2.epoch_ptp_uv(ep_ref) <= SELECT_PTP_UV
print(f"trial set fixed on the uncorrected 0.1-30 Hz epochs at {SELECT_PTP_UV:g} uV peak-to-peak: "
f"{int(KEEP.sum())} of {len(KEEP)} trials "
f"({int((ep_ref.events[KEEP, 2] == l2.P3_EVENT_ID['target']).sum())} target)")
comparison, waves, spectra = [], {}, {}
for name, (r, removed_note) in methods.items():
ep = l2.epochs_p3(r, events, tmin=TMIN, tmax=TMAX, baseline=BASELINE)[KEEP]
d = l2.difference_wave(ep)
snr = l2.erp_snr(d, CH, WINDOW, BASELINE)
waves[name] = d
x = r.get_data(picks="eeg") * 1e6
spec = r.compute_psd(method="welch", picks="eeg", fmin=1, fmax=min(60, r.info["sfreq"] / 2 - 1),
n_fft=int(4 * r.info["sfreq"]), verbose=False)
psd, freqs = spec.get_data(return_freqs=True)
spectra[name] = (freqs, 10 * np.log10(psd.mean(axis=0) * 1e12))
comparison.append({
"method": name,
f"{CH} mean {WINDOW[0] * 1000:.0f}-{WINDOW[1] * 1000:.0f} ms (uV)": snr["signal_uv"],
"baseline noise (uV)": snr["noise_uv"], "ERP SNR": snr["snr"], "epochs": int(KEEP.sum()),
"numerical rank": int(np.linalg.matrix_rank(x, tol=1e-6 * np.abs(x).max())),
"signal RMS (uV)": float(np.sqrt((x ** 2).mean())), "what was removed": removed_note,
})
print()
print(l2.fmt_table(comparison, list(comparison[0]), floatfmt="{:.2f}"))
In [8]:
fig, axes = plt.subplots(2, 2, figsize=(14, 8.5))
colors = {k: c for k, c in zip(methods, ("0.45", "tab:green", "tab:blue", "0.7", "tab:red", "tab:purple"))}
t0_s, sf = 200.0, branches["erp"].info["sfreq"]
sl = slice(int(t0_s * sf), int((t0_s + 10) * sf))
for name, (r, _) in methods.items():
axes[0, 0].plot(branches["erp"].times[sl], r.get_data(picks=["Fp1"])[0][sl] * 1e6, lw=0.8,
color=colors[name], label=name)
axes[0, 0].set(xlabel="Time (s)", ylabel="Amplitude (uV)",
title="Fp1, 10 s of continuous data under six cleanings (uV)")
axes[0, 0].legend(fontsize=7); axes[0, 0].grid(alpha=0.3)
for name in methods:
d = waves[name]
snr = [c["ERP SNR"] for c in comparison if c["method"] == name][0]
axes[0, 1].plot(d.times * 1000, d.data[d.ch_names.index(CH)] * 1e6, lw=1.4, color=colors[name],
label=f"{name} (SNR {snr:.2f})")
axes[0, 1].axvspan(WINDOW[0] * 1000, WINDOW[1] * 1000, color="tab:orange", alpha=0.18)
axes[0, 1].axhline(0, color="gray", lw=0.6); axes[0, 1].axvline(0, color="gray", lw=0.6)
axes[0, 1].set(xlabel="Time from stimulus (ms)", ylabel="Amplitude (uV)",
title=f"Target minus standard at {CH}, identical trials (uV, positive up)")
axes[0, 1].legend(fontsize=7); axes[0, 1].grid(alpha=0.3)
for name in methods:
f_, p_ = spectra[name]
axes[1, 0].plot(f_, p_, lw=1.1, color=colors[name], label=name)
axes[1, 0].set(xlabel="Frequency (Hz)", ylabel="PSD (dB re 1 uV^2/Hz)",
title="Spectrum averaged over the 30 EEG channels (dB re 1 uV^2/Hz)")
axes[1, 0].legend(fontsize=7); axes[1, 0].grid(alpha=0.3)
for base, group in (("uncorrected (0.1-30 Hz)", ["ICA (0.1-30 Hz)", "SSP (0.1-30 Hz)"]),
("uncorrected (1-30 Hz)", [k for k in methods if k.startswith("ASR")])):
f_ref, p_ref = spectra[base]
for name in group:
f_, p_ = spectra[name]
axes[1, 1].plot(f_, p_ - p_ref, lw=1.2, color=colors[name], label=f"{name}\n minus {base}")
axes[1, 1].axhline(0, color="gray", lw=0.6)
axes[1, 1].set(xlabel="Frequency (Hz)", ylabel="Change in PSD (dB)",
title="What each method removed, against its own matched baseline (dB)")
axes[1, 1].legend(fontsize=6); axes[1, 1].grid(alpha=0.3)
fig.tight_layout()
plt.show() # render the static figure(s) of this cell inline
6. The cutoff is the whole dial¶
ASR's aggressiveness is one number. Sweeping it makes the trade-off visible: how much of the data is replaced, and what happens to the ERP measured on identical trials. Both implementations are swept, because they disagree, and a reader who is going to quote a cutoff should see that the number does not mean the same thing in both.
In [9]:
CUTOFF_SWEEP = [5.0, 10.0, 20.0, 40.0]
sweep_rows = []
for backend in (["asrpy", "documented"] if asrpy is not None else ["documented"]):
for c in CUTOFF_SWEEP:
r, lg = run_asr(branches["asr"], cutoff=c, backend=backend)
ep = l2.epochs_p3(r, events, tmin=TMIN, tmax=TMAX, baseline=BASELINE)[KEEP]
snr = l2.erp_snr(l2.difference_wave(ep), CH, WINDOW, BASELINE)
sweep_rows.append({"implementation": lg["tool"].split()[0], "cutoff (SD)": c,
"replaced (uV RMS)": lg["changed_uv_rms"],
"% of the signal replaced": lg["percent_of_signal_changed"],
f"{CH} mean (uV)": snr["signal_uv"], "baseline noise (uV)": snr["noise_uv"],
"ERP SNR": snr["snr"], "seconds": lg["duration_s"]})
print(l2.fmt_table(sweep_rows, list(sweep_rows[0]), floatfmt="{:.2f}"))
print()
print("As the cutoff falls, more of the recording is replaced by its reconstruction -- and there is no list "
"anywhere of what was replaced. That is the lesson's point: ICA fails by removing a component you "
"wanted and leaves a record of it; ASR fails by removing high-amplitude activity and leaves none "
"(pf-overcleaning-ica in its ASR form).")
fig, ax = plt.subplots(1, 2, figsize=(12, 4))
for backend in sorted({r["implementation"] for r in sweep_rows}):
rowsb = [r for r in sweep_rows if r["implementation"] == backend]
ax[0].plot([r["cutoff (SD)"] for r in rowsb], [r["% of the signal replaced"] for r in rowsb], "o-",
label=backend)
ax[1].plot([r["cutoff (SD)"] for r in rowsb], [r["ERP SNR"] for r in rowsb], "o-", label=backend)
ax[0].set(xlabel="ASR cutoff (standard deviations)", ylabel="Signal replaced (%)",
title="How much ASR replaces, against its cutoff (%)")
ax[0].legend(fontsize=8); ax[0].grid(alpha=0.3)
snr_of = {c["method"]: c["ERP SNR"] for c in comparison}
ax[1].axhline(snr_of["ICA (0.1-30 Hz)"], color="tab:green", ls="--", lw=1.2,
label=f"ICA, 0.1-30 Hz ({snr_of['ICA (0.1-30 Hz)']:.2f})")
ax[1].axhline(snr_of["uncorrected (1-30 Hz)"], color="0.5", ls=":", lw=1.2,
label=f"uncorrected, 1-30 Hz ({snr_of['uncorrected (1-30 Hz)']:.2f})")
ax[1].set(xlabel="ASR cutoff (standard deviations)", ylabel="ERP SNR (dimensionless)",
title=f"ERP SNR at {CH} against ASR cutoff (dimensionless, identical trials)")
ax[1].legend(fontsize=7); ax[1].grid(alpha=0.3)
fig.tight_layout()
plt.show() # render the static figure(s) of this cell inline
7. The numbers¶
In [10]:
ASR_ROW = next(k for k in methods if k.startswith("ASR") and "documented" not in k)
ASR_DOC_ROW = "ASR, documented implementation (1-30 Hz)"
get = lambda name, key: [c for c in comparison if c["method"] == name][0][key]
AMP = f"{CH} mean {WINDOW[0] * 1000:.0f}-{WINDOW[1] * 1000:.0f} ms (uV)"
print("nb-2-7-asr-vs-ica -- L2.7 answer key (draft; TODO(confirm) at author review)")
print(f"Data: ds-erpcore P3, {SID} (CC BY-SA 4.0 per data/directory.yaml, contested at source, open, per-subject downloadable). "
f"30 EEG + 3 EOG, 1024 Hz resampled to {RESAMPLE_HZ:g} Hz, CMS online reference re-referenced to the "
f"average of the 30 EEG channels.")
print(f"Shared pipeline: pyprep bad channels (>=2 criteria, cap 10 %) -> filter -> interpolate -> average "
f"reference. {RANK['arithmetic']}.")
print(f"Bands: ICA and SSP on {ERP_BAND[0]:g}-{ERP_BAND[1]:g} Hz; ASR on {ASR_BAND[0]:g}-{ASR_BAND[1]:g} Hz, "
"because ASR's calibration is a statement about the covariance and drift dominates it otherwise "
"(section 3 measures that). Each band's uncorrected baseline is reported alongside.")
print(f"Trial set: fixed once on the uncorrected {ERP_BAND[0]:g}-{ERP_BAND[1]:g} Hz epochs at "
f"{SELECT_PTP_UV:g} uV peak-to-peak ({int(KEEP.sum())} of {len(KEEP)} trials), identical under every method.")
print(f"SNR definition: |mean amplitude of the target-minus-standard difference wave at {CH} over "
f"{WINDOW[0] * 1000:.0f}-{WINDOW[1] * 1000:.0f} ms| / SD of that waveform's averaged pre-stimulus "
f"baseline ({BASELINE[0]:g} to {BASELINE[1]:g} s). Dimensionless.")
print()
print(f"ICA: two-pass, extended Infomax on the {ICA_FIT_HZ[0]:g}-{ica_log['fit_lowpass_hz']:g} Hz branch, "
f"n_components = {RANK['rank']}, random_state {l2.SEED}; classifier {cls['tool']}; removed "
f"{len(exclude)} components "
+ (", ".join(f"IC{i} ({cls['labels'][i]})" for i in exclude) or "none") + ".")
print(f"ASR: {asr_log['tool']}, cutoff {ASR_CUTOFF:g} SD, window {asr_log.get('win_len_s', 0.5):g} s, "
f"calibration {asr_log.get('calibration', 'lowest-variance windows of this recording')}; it replaced "
f"{asr_log['percent_of_signal_changed']:.0f} % of the signal RMS.")
print()
print("ex-2-7-snr-ica -- ERP signal-to-noise ratio after ICA cleaning:")
print(f" ANSWER: {snr_of['ICA (0.1-30 Hz)']:.2f} (signal {get('ICA (0.1-30 Hz)', AMP):+.2f} uV over "
f"{get('ICA (0.1-30 Hz)', 'baseline noise (uV)'):.2f} uV of baseline noise, {int(KEEP.sum())} trials, "
f"{ERP_BAND[0]:g}-{ERP_BAND[1]:g} Hz)")
print(f" matched baseline: uncorrected {ERP_BAND[0]:g}-{ERP_BAND[1]:g} Hz = "
f"{snr_of['uncorrected (0.1-30 Hz)']:.2f}")
print("ex-2-7-snr-asr -- ERP signal-to-noise ratio after ASR cleaning of the same file, same trials:")
print(f" ANSWER: {snr_of[ASR_ROW]:.2f} (signal {get(ASR_ROW, AMP):+.2f} uV over "
f"{get(ASR_ROW, 'baseline noise (uV)'):.2f} uV of baseline noise, cutoff {ASR_CUTOFF:g} SD, "
f"{ASR_BAND[0]:g}-{ASR_BAND[1]:g} Hz)")
print(f" matched baseline: uncorrected {ASR_BAND[0]:g}-{ASR_BAND[1]:g} Hz = "
f"{snr_of['uncorrected (1-30 Hz)']:.2f}")
print(f" cross-check, the documented implementation at the same cutoff: {snr_of[ASR_DOC_ROW]:.2f} "
f"(it replaced {[r for r in sweep_rows if r['implementation'].startswith('helpers') and r['cutoff (SD)'] == ASR_CUTOFF][0]['% of the signal replaced']:.0f} % "
f"of the signal against {asr_log['percent_of_signal_changed']:.0f} % for the package). "
"TODO(confirm): the two implementations do not agree about what a cutoff of 20 SD means; a reviewer "
"should settle which behaviour the key quotes before the numbers are published.")
print(f" for reference: SSP {snr_of['SSP (0.1-30 Hz)']:.2f}")
print()
print("Full comparison:")
print(l2.fmt_table(comparison, list(comparison[0]), floatfmt="{:.2f}"))
print()
print("ASR cutoff sweep (same trials throughout):")
print(l2.fmt_table(sweep_rows, list(sweep_rows[0]), floatfmt="{:.2f}"))
print()
print("What each method leaves behind -- the ex-2-7-mobile-choice argument:")
print(f" ICA -- a list of {len(exclude)} named components with classes and probabilities; the rank falls by "
f"{len(exclude)} to {RANK['rank'] - len(exclude)} and says so. The fit needs a stationary, "
"high-channel-count recording and a matrix estimated offline, so it cannot run inside a 200 ms "
"feedback loop unless the matrix was frozen in advance.")
print(f" ASR -- no list; the output's numerical rank is {rank_asr}, indistinguishable from uncleaned data, "
"and the aggressiveness is one number whose meaning is implementation-dependent. It runs causally on "
"short windows, tolerates non-stationary movement artifact and degrades gracefully at low channel "
"counts, which is why it is the answer for a mobile, real-time, 16-channel recording -- accepting that "
"what it removed cannot be audited afterwards.")
print(f" SSP -- exactly {len(projs)} dimension(s), stated in the file and countable in the rank.")
print()
print(f"Package versions: {l2.package_versions(('asrpy',))}")