EOG regression and ICA

Why this matters

Blinks are larger than any ERP you will ever measure, and they are not random with respect to your conditions. Removing them is not optional — but every method that removes them also removes some brain activity, and the methods that remove the most artifact remove the most brain. This is the longest lesson in the level because component classification is a skill, not a setting: the drill is the point, and the reading afterwards is what turns a guess into evidence.

Concepts

Regression: the simple, honest method

If you have recorded EOG channels, the oldest correction is also the most transparent. Model each EEG channel as the sum of brain activity and a scaled copy of the EOG signal, estimate the scaling by regression, and subtract:

EEG_corrected(i, t) = EEG(i, t) − b_i · EOG(t)

with one propagation coefficient b_i per channel, estimated over the whole recording or over blink-centred segments. This is the Gratton-style family of corrections, and in MNE it is mne.preprocessing.EOGRegression.

Its virtues: one number per channel, estimable from a few blinks, no assumption about independence, and completely inspectable — you can plot the coefficients as a topography and see whether they look like a blink field.

Its limitation is the reason ICA exists: the EOG channels also record brain activity. An electrode below the eye picks up frontal cortex, so subtracting a scaled copy of it subtracts some frontal EEG too. Regression also assumes one fixed propagation pattern, whereas a blink and a horizontal saccade have different fields, and it needs EOG channels that many datasets do not have. ERP CORE ships three EOG channels, so regression is available there; many resting datasets do not.

A practical refinement: estimate the coefficients on data from which the evoked response has been removed (or on blink segments outside the task), so that condition-locked brain activity cannot leak into the estimate.

ICA at intuition level

Independent component analysis starts from a physical claim that is very nearly true for EEG: the scalp signals are an instantaneous linear mixture of a smaller number of sources, because volume conduction is linear and effectively instantaneous at EEG frequencies (L0.1, L5.1). Write it as

X = A · S

where X is channels × time (what you recorded), S is components × time (the sources’ time courses), and A is the mixing matrix whose columns are each source’s scalp projection — its topography. ICA estimates an unmixing matrix W ≈ A⁻¹ so that S = W · X.

What makes the problem solvable is the assumption that the sources are statistically independent and non-Gaussian. Independence is stronger than uncorrelatedness: PCA finds uncorrelated components ordered by variance and stops there, and its components are constrained to be orthogonal, which the real scalp projections are not. ICA keeps looking for a rotation that makes the components as independent — as non-Gaussian — as possible. A blink is an excellent candidate: it is a large, spatially fixed, temporally sparse and decidedly non-Gaussian source.

Three consequences that matter in practice:

• Order and sign are arbitrary. ICA does not rank components by importance, and the sign of a component (and of its topography) can flip between runs. Never interpret component order as importance, and never compare component numbers across subjects.
• It is not deterministic unless you make it so. Most implementations start from a random initialisation. Set random_state (and record it, L2.8) or two runs on the same data will give you two different decompositions.
• The “components” are not neurons or even cortical patches. They are directions in channel space that happen to be maximally independent. A component that looks like a dipolar field over occipital cortex is a good candidate for a real generator; a component with a speckled, non-dipolar map is a good candidate for noise. That is a heuristic, not a proof.

Algorithms. infomax (and its extended variant, the EEGLAB default) maximises an information-theoretic objective; fastica uses a fixed-point iteration on negentropy and is fast but more sensitive to initialisation; picard is a modern preconditioned solver that converges to the same kind of solution faster and is MNE’s recommended default in most cases. For component classification the differences matter far less than data quantity and rank do.

Rank and data quantity: the two things that actually break ICA

Rank. 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 estimator responds by manufacturing components out of numerical noise: near-duplicate pairs, components with implausible maps, or outright failure to converge. Your rank is the channel count minus one per interpolated channel (L2.2), minus one for an average reference (L2.3), minus anything else linear you did. Compute it, pass it, and record it:

rank = mne.compute_rank(raw, rank='info')          # or carry your own arithmetic
ica = mne.preprocessing.ICA(n_components=rank['eeg'], method='picard', random_state=97)
ica.fit(raw_for_ica)                                # the 1 Hz copy, L2.4

That is pf-interpolation-rank, below, and it is the single most common way to get a bad decomposition from good data.

Data quantity. Estimating an N×N unmixing matrix needs enough samples to estimate N² parameters. The rule of thumb in the literature is that the number of samples should be some multiple of the squared channel count — (Chaumon, 2015) discusses the requirement and its consequences. TODO(confirm) the multiplier to quote: values in the range of tens are commonly cited, and the practical test is whether the decomposition is stable when you re-run it on half the data. Short recordings and high channel counts are the bad combination; reducing rank deliberately (a PCA step before ICA) is a legitimate response, and it must be reported because it changes what can be separated.

Filtering. Fit on the 1 Hz-filtered copy and apply to the 0.1 Hz analysis data, exactly as L2.4 sets out. Also: fit on continuous data rather than epochs where you can, and exclude the grossly artifactual segments first — a single enormous pop can dominate the objective and produce a component that models nothing but that pop.

Classifying a component: four views and a rationale

A component is identified from the convergence of several views, never from one. The gallery gives you four:

Topography. Where the component projects on the scalp. A vertical blink is a broad frontal-polar field with the same polarity at Fp1 and Fp2; a horizontal saccade is a lateral frontal field with opposite polarity on the two sides; muscle is focal at the rim (temporal, neck) with a sharp, non-dipolar edge; a cardiac component is a broad, shallow gradient across the whole head; a channel-noise component is a single hot electrode with nothing around it; a brain component is smooth and dipolar.

Time course. Blinks are isolated, stereotyped, large deflections lasting a few hundred milliseconds. Saccades are step-like. Muscle is dense bursts tied to jaw or neck. Cardiac is periodic at roughly 1 Hz with a recognisable sharp beat. A brain component is continuous and modulated, without stereotyped repeats.

Spectrum. Brain components usually show a 1/f-like background with a peak — alpha near 10 Hz, mu with its arch harmonic. Muscle power rises with frequency across the band. Line noise is a single spike at 50 or 60 Hz and its harmonics. A blink component’s power is concentrated below a few hertz.

ERP image. Single trials stacked as rows, sorted by trial or by reaction time. A brain component that carries the effect shows a consistent deflection aligned to the event across trials; a blink component shows bright bands at whatever latency the participant blinked, often after the response; a muscle component shows scattered noise.

The classes worth distinguishing are the ICLabel classes: brain, muscle, eye, heart, line noise, channel noise, other. They are the labels the drill uses.

ICLabel as a second opinion

(Pion-Tonachini et al., 2019) trained a classifier on a large crowd-labelled set of components and released it as ICLabel; mne-icalabel is the Python port. It returns a probability per class per component, which is genuinely useful and genuinely limited:

• Use it as a second opinion and a triage, not an oracle. Components it labels with high confidence and you agree with are settled; components where you and it disagree are the ones to look at properly.
• Its training data came from particular montages, references and filter settings; it is least reliable on data far from those conditions (TODO(confirm) which conditions the released model documents as its training regime).
• A probability is not a decision. What you remove is a decision, and the reason you removed it belongs in the QC report (L2.8).

If mne-icalabel cannot be installed in your environment, the site’s component gallery labels its components from documented feature heuristics instead, and says so in its data file and its README.

Over-cleaning

The error in the other direction is quieter and more damaging. Every component you remove is subtracted from your data, and if it carried brain activity, so is your effect. The classic case: a posterior component with a clear alpha peak and a dipolar occipital map, removed because “alpha is noise” — and with it, the oscillatory effect the study was about. A second case: removing many “other” components on the grounds that they are not obviously brain, when “other” is precisely the class that means “the classifier cannot tell”.

Symptoms of over-cleaning: the effect shrinks after cleaning; the number of removed components varies wildly between subjects; the removed set includes components with dipolar maps and physiological spectra; the cleaned data has visibly lower rank than it should. That is pf-overcleaning-ica, below.

The discipline that prevents it: decide a policy before you look (which classes you remove, and on what evidence), keep the number removed small and comparable across subjects, record for each removed component the class, the evidence and the ICLabel probability, and — the test that settles it — compute your effect with and without the cleaning step and report both if they differ.

The data behind this lesson

• The component gallery ships at least 40 real components from one ds-erpcore subject, each with its topography weights, a 10-second time course, a spectrum, an ERP-image thumbnail, ICLabel probabilities where the classifier was available, and a rationale listing the evidence. ERP CORE is CC BY 4.0, open access and per-subject downloadable; TODO(confirm): the author mirrors the entry into the catalog registry and signs off the dataset page (§10.11 item 8).
• The gallery’s labels are label_source: algorithmic (§4.5) until the author reviews them. The drill grades against those labels; a disagreement is a disagreement with a documented procedure, and the rationale text tells you which evidence the procedure used.

Explore

Before you start scoring, browse. Find one unambiguous example of each class and look at all four views together, so you know what “obvious” looks like. Then start the drill, and for every component force yourself to name the evidence before you press a button.

Practice

The notebook fits ICA on a 1 Hz-filtered copy with the correct rank and a fixed seed, applies the unmixing to the 0.1 Hz analysis data, runs mne-icalabel where available, and plots the ERPs before and after cleaning so that the cost of each removed component is visible. Its final cell prints the components removed, their classes and probabilities, and the effect size before and after.

Exercises