Stage 6: ICA artefact decomposition

Author: Hadi Zaatiti <hadi.zaatiti@nyu.edu>

What & why

Stage 6 of the KIT general pipeline decomposes every stage-2 recording with independent component analysis (ICA) so that stereotyped artefacts (heartbeat, eye blinks and movements, muscle) can be removed before the data is filtered, epoched and averaged. It follows the recipe of the COGNESTIC course (08_EEG_MEG/A_1_Artifact_Correction.ipynb, itself based on the MNE-Python ICA tutorial):

  1. fit the ICA on a copy of the data high-pass filtered at 1 Hz (slow drifts are detrimental to the decomposition; this stage also low-passes that copy at 40 Hz, see step 1 below), with ICA(n_components=0.90, method="fastica", max_iter="auto", random_state=97);

  2. look at the solution: plot_sources, plot_components, plot_properties, plot_overlay, the explained variance;

  3. let the data propose artefact components: find_bads_ecg on a synthetic ECG (correlation, threshold auto), find_bads_eog when an EOG or proxy channel exists, with plot_scores and the sources of the heartbeat-locked average as evidence;

  4. set ica.exclude and apply the solution to the unfiltered data.

As in stage 5, no component is removed without a human: the automatic matches are proposals, pre-marked in MNE’s source browser next to the figures above; the reviewer confirms, rejects or extends them and closes the window, and that is the decision. Stage 7 then applies the reviewed solution to the unfiltered recording right before filtering, so no second copy of the raw data is written.

The decomposition uses the 207 KIT head gradiometers (MNE type mag) minus the bad channels decided by stage 5. The reference magnetometers (ref_meg, including MEG 092) never take part.

What the ICA computes

Notation: the fit copy is \(X \in \mathbb{R}^{n \times T}\), \(n\) channels (207 gradiometers minus the bad ones) by \(T\) samples (the band-passed recording, every decim-th sample: about 510 000 samples for a 42-minute KidLang split). MNE’s steps, in order:

  1. Scaling and centering. \(\sigma\) = standard deviation of all entries of \(X\) (one scalar, since there is one channel type); \(X_1 = X/\sigma\). \(m \in \mathbb{R}^n\) = mean of each row of \(X_1\); \(X_2 = X_1 - m\,\mathbf{1}^{\top}\).

  2. PCA. Thin SVD \(X_2^{\top} = U\,\Sigma\,V^{\top}\); \(P := V^{\top} \in \mathbb{R}^{n \times n}\) has orthonormal rows (\(P P^{\top} = P^{\top} P = I_n\)), \(\lambda_j = \Sigma_{jj}^2/(T-1)\) decreasing. \(k\) = smallest integer with \(\sum_{j \le k} \lambda_j / \sum_j \lambda_j > 0.9\) (n_components: 0.9; 23 for Y0312 split 1). \(P_k\) = first \(k\) rows, \(P_{\perp}\) = the other \(n-k\). Whitened scores \(Z = \Lambda_k^{-1/2} P_k X_2 \in \mathbb{R}^{k \times T}\), with \(Z Z^{\top}/(T-1) = I_k\).

  3. FastICA. Find an orthogonal \(R \in \mathbb{R}^{k \times k}\) such that the rows of \(S = R\,Z\) are as non-Gaussian, hence as statistically independent, as possible: maximise \(\sum_i \big|\,\mathbb{E}[G(r_i^{\top} z)] - \mathbb{E}[G(\nu)]\,\big|\) with \(G(u) = \log\cosh u\), \(\nu \sim \mathcal{N}(0,1)\), subject to \(R R^{\top} = I\). Justification: a mixture of independent non-Gaussian sources is more Gaussian than the sources, so the most non-Gaussian orthogonal directions of the whitened data are the unmixed ones. Algorithm (scikit-learn, algorithm="parallel"): \(R\) starts random (random_state: 97); fixed-point update \(R \leftarrow \mathbb{E}[g(RZ)\,Z^{\top}] - \mathrm{diag}(\mathbb{E}[g'(RZ)])\,R\) with \(g = \tanh\), then symmetric orthogonalisation \(R \leftarrow (R R^{\top})^{-1/2} R\); stop when every row moved by less than \(10^{-4}\) or after 1000 iterations (max_iter: auto). Only expectations over the \(T\) samples enter, so the order of samples is irrelevant and decimation is harmless. There is no reconstruction error to minimise: every invertible \(R\) reconstructs \(Z\) exactly.

  4. Stored solution. \(W = R\,\Lambda_k^{-1/2}\) (unmixing_matrix_, \(k \times k\)) and \(M = W^{-1} = \Lambda_k^{1/2} R^{\top}\) (mixing_matrix_). Component time courses: \(S = W P_k X_2\) (row \(i\) = component \(i\), the trace in the source browser). Topography of component \(i\): \(a_i = \sigma\,P_k^{\top} M e_i \in \mathbb{R}^n\) (the map in plot_components). The fit data decompose exactly, with no error term:

    \[X = \sum_{i=1}^{k} a_i\, s_i^{\top} + \sigma\left(P_{\perp}^{\top} P_{\perp} X_2 + m\,\mathbf{1}^{\top}\right),\]

    the second term being the 10 % of the variance left out of the ICA, orthogonal to every component. Identifiability: the components are defined up to permutation and sign (their scale is fixed by the whitening), which is what the seed pins down. Assumptions: mutually independent sources, at most one Gaussian, a mixing that does not change during the recording (sensor and source geometry), \(n \ge k\).

  5. Cleaning (stage 7, ica.apply, all \(n\) PCA components kept). For the unfiltered recording \(X' \in \mathbb{R}^{n \times T'}\) (same channels), \(S' = W P_k (X'/\sigma - m\,\mathbf{1}^{\top})\) with the same \(\sigma, m, P_k, W\); nothing is refitted. With \(E\) the excluded set,

    \[\tilde X' = X' - \sum_{i \in E} a_i\, s_i'^{\top} = X' - \sigma\,P_k^{\top} M_{:,E}\, W_{E,:}\, P_k \left(X'/\sigma - m\,\mathbf{1}^{\top}\right),\]

    which is MNE’s \(\sigma[(I - P_k^{\top} M_{:,E} W_{E,:} P_k)(X'/\sigma - m\mathbf{1}^{\top}) + m\mathbf{1}^{\top}]\) simplified with \(P_k^{\top} M W P_k + P_{\perp}^{\top} P_{\perp} = I_n\). Everything outside \(\mathrm{span}\{a_i : i \in E\}\) is untouched, including the variance the ICA never saw. The operation is linear and sample-wise, which is why a \(W\) estimated on the filtered copy (where drifts and high-frequency noise do not dominate the statistics) applies unchanged to the unfiltered data.

What you need to provide

  • The stage-2 outputs (kit2fiff_summary.csv and the FIFs).

  • The stage 5 records: a recording without one is decomposed with the FIF’s own info['bads']; a recording whose candidates are still pending is decomposed with only its human-provided bad channels (the stage warns in both cases).

  • A desktop session with a Qt binding and mne-qt-browser for the review; scikit-learn for FastICA (both in the neurowaves recipe). Fitting needs memory: the 207 channels of a 42-minute KidLang recording take about 8 GB while they are loaded and filtered.

Configuration file

The ica block of the YAML config (defaults shown; the KidLang settings are the same):

ica:
  enabled: true
  fit:
    l_freq: 1.0            # high-pass of the copy the ICA is fitted on (COGNESTIC)
    h_freq: 40             # low-pass of that copy (null = none); see step 1
    n_components: 0.9      # fraction of variance to keep (< 1) or a number of components
    method: fastica
    max_iter: auto
    random_state: 97       # fixed seed: reruns give identical components
    decim: 5               # use every 5th sample when fitting (speed)
    max_duration_s: 0      # fit on the first N seconds only (0 = whole recording)
    reject_t: null         # peak-to-peak rejection while fitting, in T (e.g. 5.0e-12)
  ecg:
    enabled: true
    method: correlation    # as in COGNESTIC ('ctps' is MNE's other option)
    measure: zscore        # 'zscore' = outliers among the components; 'correlation' = absolute value
    threshold: auto        # 3 standard deviations with zscore; e.g. 0.9 with correlation
    ch_name: null          # a real ECG channel name, if one was recorded
  eog:
    enabled: false         # KIT records no EOG; set the channels to use (e.g. MISC lines)
    channels: []
    threshold: 3.0
  muscle:
    enabled: false         # find_bads_muscle (MNE >= 1.1), opt-in
    threshold: 0.5
  manual_exclude: {}       # components always removed, per subject: {Y0312: [3, 7]}
  review:
    n_cuts: 8              # cuts of the filtered fit copy shown in the browser
    cut_duration_s: 30
    skip_edges_s: 5
    resample_hz: 150
    window_s: 10           # seconds per browser page
    antialiasing: false
    use_opengl: false
    whole_recording: false # true = browse the entire filtered recording (slow to open)
    max_properties: 8      # property figures opened for the proposed components
  save_plots: true
  save_cleaned: false      # also write <entities>_desc-icaclean_raw.fif (large)

manual_exclude is keyed by subject label with or without sub- and applies to every recording of that subject. reports.scripts.6-kit_ica enables or disables the report refresh after the stage.

How to run

From the repository root with the neurowaves environment activated:

python -m pipeline.mne_pipelines.kit_general_pipelines.6-kit_ica \
    -c pipeline/mne_pipelines/kit_general_pipelines/pipeline_config_files/kidlang_config.yml

Options:

Flag

Effect

-c, --config PATH

YAML config (always pass the dataset’s file).

(no flag)

Interactive run: fit (or reuse the saved solution), propose, open the figures and the source browser for every recording without a saved review. While you review one recording, the next one is fitted in a background process, so after the first there is no wait.

--no-gui

Fit, propose and save the figures only; exclude nothing (review_status: pending). Needs no display; use it to precompute the solutions, the GUI run then opens straight into the review.

--review-all

Redo the review of reviewed recordings, starting from the previous decision.

--reset-review

Discard the saved reviews; recordings become pending again.

--overwrite

Refit even when a solution with the same input, settings and bad channels exists (reviews are kept but worth redoing).

--no-plots

Skip the saved figures (the live figures of the review still open).

--save-cleaned

Also write <entities>_desc-icaclean_raw.fif, the unfiltered recording with the reviewed components removed (as large as the input). Stage 7 does not need it.

A solution is reused when the FIF (size and content digest; the modification time is only a fast path, so a synced copy still matches), the fit/ecg/eog/muscle settings and the bad-channel list are unchanged. A changed ecg/eog/muscle setting recomputes the proposals on the same solution and keeps the reviewer’s decision (new proposals are put to the reviewer); only a refit (new input, fit settings or bad channels) sets the review aside. The FIF’s own info['bads'] always stay out of the decomposition, with or without a stage-5 record. The VS Code entries (Kidlang) ICA: fit + review (GUI), (Kidlang) ICA: fit only (no GUI) and (Kidlang) ICA: redo all reviews (GUI) run the three variants.

Step by step

Step 1: fit on a filtered copy

The 207 gradiometers minus the bad channels (plus the configured EOG or proxy channels, if any) are loaded, band-pass filtered between fit.l_freq and fit.h_freq (1 to 40 Hz; one zero-phase FIR pass, raw.filter(l_freq, h_freq), which is the same as a 1 Hz high-pass followed by a 40 Hz low-pass: both are zero-phase, and the two kernels combine to the same frequency response), and the ICA is fitted on them with decim (every n-th sample), optionally on the first max_duration_s seconds and with peak-to-peak rejection reject_t. COGNESTIC only high-passes at 1 Hz, but its sample data was already low-passed at 40 Hz; KIT raw data is broadband at 1000 Hz, and without the low-pass the sensor noise above 40 Hz costs dozens of extra components (measured on a 42-minute KidLang recording: 90 components for 10 minutes of broadband data and a fit that did not finish in 17 minutes on the whole recording, against 23 components and 43 s with the 1 to 40 Hz band). decim subsamples without an anti-alias filter, which is another reason to low-pass first. n_components: 0.9 keeps as many components as needed to explain 90 % of the variance; the fixed random_state makes the solution reproducible. The explained-variance ratio of the solution and of each component is computed on a 60 s crop and stored. Loading and filtering one recording takes about a minute and 8 GB of memory.

Step 2: the figures

Saved under figures/ and, during the review, opened live:

  • plot_components: the topographies of all components (several pages of 20); live, clicking a topography opens its properties;

  • plot_scores: the match score of every component for each proposer, proposed components in red;

  • plot_properties of each proposed component (up to review.max_properties): topography, source segments, spectrum, variance;

  • plot_sources of the ECG-evoked response (the sources averaged around the detected heartbeats): heartbeat components show a clear QRS shape;

  • plot_overlay: the data before (red) and after (black) removing the proposals, and after the review the same with the final exclusions (overlay-final).

Step 3: proposals

  • ECG: find_bads_ecg correlates every component with the ECG signal. KIT records no ECG channel, so MNE builds a synthetic one from the average of the MEG channels (ch_name: null), which contains the cardiac field clearly enough. With measure: zscore (default) the correlations are z-scored and threshold: auto proposes the components more than 3 standard deviations from the others; with measure: correlation the threshold is an absolute correlation (for example 0.9). The heartbeat often appears as two components (the cardiac field is a travelling dipole). On KIT data the synthetic ECG is a weak reference (the cardiac field is a small part of the MEG average), so a heartbeat component may not stand out statistically: the stage therefore also prints and stores the three highest ECG correlations as hints, and the ECG-evoked sources figure shows which components carry the QRS shape. measure: correlation with threshold: 0.5 proposes more liberally.

  • EOG: find_bads_eog against the configured channels, only when eog.enabled and eog.channels are set. Without an EOG channel the blink components are not proposed, but they are easy to spot in the review: a frontal, dipolar topography and a blink-shaped time course.

  • Muscle: find_bads_muscle when enabled (rarely needed for MEG).

  • manual_exclude components are always removed and pre-marked.

The union is the list of candidates; they are proposals and nothing more.

Step 7: human review and decision

What you decide. Which components are artefacts (heartbeat, eye blinks and movements, muscle, a jumping sensor) rather than brain activity. Every excluded component is subtracted from the recording for all later stages, so exclude a component only when at least two kinds of evidence agree (its map and its time course, or its map and its spectrum or heartbeat average); when in doubt, keep it. A clean KIT recording typically yields 1 to 4 artefact components.

The windows, in the order they open (all titled with the recording):

  1. Terminal: the proposals of step 3 and, on KIT, the three components with the highest correlation to the synthetic ECG, as hints to check first.

  2. components (click one for its properties): a grid of all component maps (ICA000, ICA001, …), each the view of the sensor array from above. Heartbeat: a smooth gradient across the whole array, usually strongest at the lower edge nearest the heart, often two components. Blinks: a strong frontal pattern, symmetric left/right. Horizontal eye movements: frontal, antisymmetric left/right. Muscle: patchy, confined to edge sensors (temples, neck). A single-sensor problem: one isolated hot spot. Brain components: dipolar patterns over the scalp, often with left/right counterparts. Click a map to open its properties window.

  3. Properties window (one per clicked or proposed component): the map, an image of its time course cut into consecutive 2 s segments and stacked (one row per segment, time left to right), the mean segment, the power spectrum and the variance per segment. Heartbeat: regular vertical stripes in the image, a spectrum with a peak near 1 Hz and its harmonics. Blinks: a few large slow deflections, power concentrated below 4 Hz. Muscle: power rising towards 40 Hz. Brain: a 1/f spectrum, usually with an alpha bump near 10 Hz, no stripes.

  4. ECG match scores (when the proposer ran): one bar per component, its correlation with the synthetic ECG; red bars were proposed. On KIT the synthetic ECG is weak, so read the bars as a ranking, not a verdict.

  5. sources of the ECG-locked average: every component averaged around the detected heartbeats; a heartbeat component shows a sharp deflection at time 0, the others stay flat.

  6. Overlay (when there are proposals): the sensor data before (red) and after (black) removing the proposals; the difference should look like the artefact alone.

  7. Review ICA: ..., MNE’s source browser (plot_sources, Qt backend from mne-qt-browser), opened last; the stage waits for it. It shows the component time courses on review.n_cuts cuts of review.cut_duration_s seconds taken evenly across the recording (an annotation names each cut and where it comes from; 4 minutes in total by default, resampled to review.resample_hz; whole_recording: true shows everything, slowly), review.window_s seconds per page, with the proposed and manual components already greyed out. Keys: left/right arrows scroll in time, up/down arrows (or Page Up/Page Down) change the component page, +/- change the scale. Heartbeat: sharp spikes every 0.7 to 1.2 s, perfectly regular. Blinks: large slow bumps of about 300 ms every few seconds. Eye movements: step-like changes. Muscle: bursts of fast activity. Sensor jump: a square step. Brain: irregular waves with alpha bursts, no fixed rhythm.

How to mark a component. In the browser, click the component name in the left margin: grey = excluded, black = kept; clicking again reverses it. Right-click the name to open that component’s properties window. Compare what you see with the hints and the proposals, and check the components the browser pre-marked: they are proposals, unmark any you disagree with.

How to finish. Close the browser window. The components greyed at that moment, plus manual_exclude, are the final exclusions (review_status: reviewed); the stage prints Decision: excluded [...], saves the solution with them (ica.save), updates the scores table, draws the final overlay, closes the figures and opens the next recording.

With --no-gui the review does not run: review_status stays pending, exclude holds only the manual components, and stage 7 warns that proposals are waiting. Proposals that appear after a review (changed proposer settings) are put to the reviewer on the next interactive run; a headless run records them as unseen_candidates.

Outputs

<bids_root>/derivatives/6-kit_ica/
  config_<timestamp>.yml
  ica_summary.csv                       # one row per recording, whole cohort, rewritten after every recording
  sub-<id>/
    ica_log.txt
    [ses-<id>/]
      sub-<id>..._desc-ica_ica.fif        # the ICA solution incl. exclude (mne.preprocessing.read_ica)
      sub-<id>..._desc-ica.json           # proposals, scores, review, settings, provenance
      sub-<id>..._desc-ica_scores.tsv     # one row per component
      sub-<id>..._desc-icareview_raw.fif  # review data, kept only until the review is done
      sub-<id>..._desc-icaclean_raw.fif   # only with --save-cleaned
      figures/
        sub-<id>..._desc-ica_components-<n>.png
        sub-<id>..._desc-ica_scores-ecg.png        (and -eog, -muscle when enabled)
        sub-<id>..._desc-ica_properties-ICA<nnn>.png
        sub-<id>..._desc-ica_ecg_sources.png
        sub-<id>..._desc-ica_overlay.png
        sub-<id>..._desc-ica_overlay-final.png

The JSON record holds exclude and review_status (what stage 7 uses), candidates, proposals (per proposer: indices, scores, method, threshold), manual_exclude, unseen_candidates, review (pre-marked, decision, timestamp, browser backend, data shown, started_from), n_components, explained_variance, variance_per_component, channels_fit, bad_channels (names, status, record used), fit_settings, detection_settings, input (FIF identity), ica_fif, scores_tsv, figures, fit_duration_s and mne_version.

The scores TSV has, per component: variance_ratio, ecg_score, eog_score, muscle_score, the *_candidate flags, candidate, manual and excluded.

ica_summary.csv has one row per recording: the BIDS entities, in_fif, n_channels_fit, bad_channels, bad_channel_status, fit_duration_s, n_components, explained_variance, the *_candidates lists, candidates, manual_exclude, review_status, excluded, n_excluded, ica_fif, record_json, scores_tsv, cleaned_fif, status, error and mne_version.

How later stages use the result

utilities.final_ica_for_recording(bids_root, entities, source_fif) returns (ica_fif, exclude, status, record_path) with status ok, pending_review (only the manual components in exclude), stale_input (record made from another version of the FIF), missing_solution or no_record. Stage 7 calls it after resolving the bad channels, warns on pending_review and stale_input, applies ica.apply with the exclusions before the notch and band-pass filters (COGNESTIC applies the solution fitted on filtered data to the original data; filtering afterwards is a linear operation and commutes with it), stores ica_excluded, ica_source, ica_record and ica_fif in the provenance sidecar and the ica_excluded / ica_source columns of filter_summary.csv, and reports an existing filtered FIF whose exclusions differ as existing_output_unverified (rerun stage 7 with --overwrite). Nothing is applied when the record excludes no component.

Report

The section 6. ICA artefact decomposition of the Quarto report explains the four steps with the parameter values of the run, shows the cohort table (reviewed and pending counts, components removed per recording) and, per recording, a collapsed block with the proposals, the decision, the scores of the proposed and excluded components and the figures. The report is refreshed in the background after every recording.

Edge cases & safeguards

  • Short recordings or no channels: a recording shorter than 10 s, or one with fewer than two gradiometers left after the bad channels, is refused with a clear error (status: error); the review cuts of a short recording shrink to what exists.

  • FIF changed since the fit (stage 2 rerun): the solution is refitted on the next run; stage 7 reports stale_input until then. A file whose content is the same (a sync client rewrote its modification time) is not refitted.

  • Bad-channel list changed (stage 5 rerun): the solution is refitted, since the channel set changed; the previous review is set aside with it (its component indices refer to other components) and the GUI starts from the new proposals.

  • FastICA does not converge: MNE warns; raise fit.max_iter or lower n_components.

  • Interrupted run: finished recordings keep their solution, record and summary row; a background fit of the next recording is terminated and redone on the next run.

  • No Qt binding / no display: the interactive run stops at start; use --no-gui.

Troubleshooting

  • Fitting is slow or yields very many components: make sure fit.h_freq is set (40 Hz); broadband data needs several times more components and the fit can run for tens of minutes. fit.max_duration_s (for example 600 s) fits on the first 10 minutes only. The console prints the fit time and the number of components.

  • No ECG proposal although a heartbeat component is visible: with measure: zscore a component is proposed only when it stands out from the others; try measure: correlation with threshold: 0.9, or just mark the component in the review (the ECG-evoked sources figure shows which ones carry the QRS shape).

  • Too many proposals: raise the threshold, or ignore them in the review; only what is marked at close is removed.

  • The browser is slow: same remedies as for stage 5 (review.antialiasing: false, review.use_opengl: true on a machine with a GPU, mne-qt-browser installed, no debugger).