Independent Component Analysis (ICA)

Initialize NeuroAnalyzer
using NeuroAnalyzer
using CairoMakie
CairoMakie.activate!(type = "png")
using Statistics
NeuroAnalyzer.progress_bar = false

# Makie figure and axis settings matching NeuroAnalyzer plots
fig_size = (900, 450)
ax_style = (titlesize = 18,
            xlabelsize = 18,
            ylabelsize = 18,
            xticklabelsize = 12,
            yticklabelsize = 12,
            xticks = LinearTicks(10),
            xminorticksvisible = true,
            xminorticks = IntervalsBetween(10),
            xautolimitmargin = (0, 0),
            yautolimitmargin = (0.1, 0.1))

# EEG with EOG and ECG channels; high-pass 1 Hz (recommended for ICA)
eeg = import_edf("files/eeg/eeg.edf")
NeuroAnalyzer.filter!(eeg;
                      ch = "all",
                      fprototype = :butterworth,
                      ftype = :hp,
                      cutoff = 1,
                      order = 4)
NeuroAnalyzer.filter!(eeg;
                      ch = "all",
                      fprototype = :fir,
                      ftype = :lp,
                      cutoff = 35,
                      bw = 10)
load_locs!(eeg; file_name = "files/locs/standard-10-20-cap19-elmiko.locs")
# 2-minute segment
eeg = crop(eeg; seg = (100, 220))

Independent Component Analysis (ICA) separates a set of mixed signals into statistically independent components. In EEG, each electrode records a mixture of the activity of many sources (brain areas, eyes, muscles, heart) conducted through the head; ICA estimates these sources without knowing them or the way they are mixed (blind source separation). Its main use in EEG is artifact removal: components representing eye blinks and movements, muscle activity or ECG are identified and removed, and the signal is reconstructed from the remaining components.

The ICA Model

The recorded signal \(X\) (channels × samples) is modeled as a linear mixture of \(n\) independent components \(S\) (\(n\) × samples):

\[ X = A \, S \]

where \(A\) (channels × \(n\)) is the mixing matrix: column \(i\) of \(A\) is the topography of component \(i\) (its weight at each channel). ICA finds the unmixing matrix \(W = A^{-1}\) such that the components \(S = W X\) are as independent as possible.

Assumptions:

  • the sources are statistically independent (not only uncorrelated, as in PCA),
  • the sources are non-Gaussian (at most one Gaussian source can be separated); ICA finds them by maximizing non-Gaussianity,
  • the mixing is linear and instantaneous, and the sources do not move during the recording (fixed topographies),
  • the number of components is at most the rank of the data (the number of independent channels).

NeuroAnalyzer uses the FastICA algorithm: the data are centered and whitened (decorrelated, with PCA), and the unmixing matrix is found iteratively by maximizing an approximation of negentropy (a measure of non-Gaussianity) with the nonlinearity f (:tanh, default, or :gaus).

Preparing the Data

  • High-pass filter at 1 Hz: slow drifts are large, non-stationary and dominate the decomposition.
  • Remove bad channels and large non-stereotyped artifacts (movements, electrode pops): ICA can model only a limited number of sources, and rare large artifacts use up components.
  • Use enough data: the number of samples should be much larger than the number of components squared (a common rule: at least 20 × \(n^2\), here 20 × 19² ≈ 7200 samples, 28 s at 256 Hz). The recording must be continuous.
  • Check the rank: after average referencing or interpolating channels, the data have fewer independent channels than channels; use at most as many components as the rank (see PCA).

Decomposition

ic, ic_mw, ic_var = ica_decompose(eeg; ch = "eeg")
  • ic: the independent components (components × samples),
  • ic_mw: the mixing matrix \(A\) (channels × components),
  • ic_var: the variance of the data accounted for by each component (%); the components are sorted by it.

By default, as many components as channels are computed (n). FastICA starts from a random unmixing matrix: the seed rng makes the result reproducible; with another seed, the components may come out in a different order and with a different sign (the sign and scale of a component are arbitrary: they are shared between the component and its topography).

for idx in 1:5
    println("IC ", idx, ": ", round(ic_var[idx]; digits = 1),
            "% of the variance")
end
IC 1: 24.9% of the variance
IC 2: 19.9% of the variance
IC 3: 12.2% of the variance
IC 4: 7.0% of the variance
IC 5: 6.5% of the variance

Identifying Artifact Components

Artifact components are identified by their time course, topography and spectrum, and, if available, by their correlation with reference channels (EOG, ECG):

Artifact Time course Topography Spectrum
eye blinks large, slow peaks (200-400 ms) frontal (Fp1, Fp2), symmetric low frequencies
horizontal eye moves steps frontal, opposite signs left/right (F7, F8) low frequencies
muscle bursts of fast activity temporal, often one electrode broad, > 20 Hz
ECG sharp peaks at the heart rate broad, often left-right gradient harmonics of the heart rate
bad electrode noise or steps in one channel one electrode broad

This recording has two EOG channels and an ECG channel. The correlation of each component with them:

eog = extract_channel(eeg; ch = "EOG1")[1, :, 1]
ecg = extract_channel(eeg; ch = "ECG")[1, :, 1]
r_eog = [abs(cor(ic[idx, :], eog)) for idx in axes(ic, 1)]
r_ecg = [abs(cor(ic[idx, :], ecg)) for idx in axes(ic, 1)]
k_eog = argmax(r_eog)
println("largest |r| with EOG: IC ", k_eog, " (",
        round(r_eog[k_eog]; digits = 2), ")")
println("largest |r| with ECG: IC ", argmax(r_ecg), " (",
        round(maximum(r_ecg); digits = 2), ")")
largest |r| with EOG: IC 2 (0.77)
largest |r| with ECG: IC 1 (0.14)

One component follows the EOG closely; no component follows the ECG (the ECG artifact is weak in the EEG channels of this recording).

Time course of the EOG component and the EOG channel:

t = eeg.time_pts
idx = round(Int, 20 * sr(eeg)):round(Int, 40 * sr(eeg))
fig = Figure(size = fig_size)
ax1 = Axis(fig[1, 1];
           ax_style...,
           title = "IC $k_eog",
           ylabel = "[a.u.]")
lines!(ax1, t[idx], ic[k_eog, idx]; color = :black)
ax2 = Axis(fig[2, 1];
           ax_style...,
           title = "EOG1",
           xlabel = "Time [s]",
           ylabel = "[μV]")
lines!(ax2, t[idx], eog[idx]; color = :black)
linkxaxes!(ax1, ax2)
fig

Its topography (the column of the mixing matrix):

plot_topo(eeg;
          data = ic_mw[:, k_eog],
          ch = "eeg",
          title = "IC $k_eog",
          ps = :m)

The left frontal topography (Fp1, F7, near the left eye recorded by EOG1) and the time course correlated with the EOG (both show the eye movements at 29-36 s) identify it as an eye component.

Removing Components

ica_remove() subtracts the back-projection of the selected components (\(A_{:,i} S_{i,:}\)) from the channels; the rest of the signal (including the channel means) is kept. With keep = true, only the selected components are kept and all others are removed.

eeg_clean = ica_remove(eeg;
                       ch = "eeg",
                       ic = ic,
                       ic_mw = ic_mw,
                       ic_idx = k_eog)
for ch in ["Fp1", "Fp2", "O1"]
    x0 = extract_channel(eeg; ch = ch)[1, :, 1]
    x1 = extract_channel(eeg_clean; ch = ch)[1, :, 1]
    println(rpad(ch, 4), "|r| with EOG: ", round(abs(cor(x0, eog)); digits = 2),
            " → ", round(abs(cor(x1, eog)); digits = 2))
end
Fp1 |r| with EOG: 0.71 → 0.09
Fp2 |r| with EOG: 0.22 → 0.1
O1  |r| with EOG: 0.18 → 0.03
fp1 = extract_channel(eeg; ch = "Fp1")[1, :, 1]
fp1_c = extract_channel(eeg_clean; ch = "Fp1")[1, :, 1]
fig = Figure(size = fig_size)
ax = Axis(fig[1, 1];
          ax_style...,
          title = "Fp1",
          xlabel = "Time [s]",
          ylabel = "Amplitude [μV]")
lines!(ax, t[idx], fp1[idx]; color = (:gray, 0.6), label = "original")
lines!(ax, t[idx], fp1_c[idx]; color = :black, label = "EOG component removed")
axislegend(ax; position = :rb)
fig

The correlation with the EOG drops in all channels, most in Fp1, where the eye artifact was largest. In this recording the eye artifact in the EEG channels is small (Fp1: about 20 μV, compared with 100 μV in the EOG channel), so the difference between the signals is modest; it is largest during the eye movements at 29-36 s. ica_remove!() modifies the object in place.

Tip: Removing a component removes all its activity, including any brain activity mixed into it. Remove only components that are clearly artifacts, and check the result (e.g. by comparing spectra or ERPs before and after).

Limitations

  • ICA separates sources only as well as they are independent, non-Gaussian and stable in space; artifacts that change their topography (e.g. head movements) or are rare are not well separated.
  • The number of components is limited by the number of channels and the rank of the data: with 19 channels, many brain and artifact sources share components.
  • The decomposition depends on the data segment, the preprocessing and the random seed; the components have no natural order, sign or scale.
  • In NeuroAnalyzer, the components are computed from the continuous data passed to ica_decompose(); to clean data processed differently (e.g. filtered at 0.1 Hz instead of 1 Hz), decompose that data.