9.2: Exercise- A First Pass at ICA-Based Blink Correction

I find that it’s easiest to understand ICA by starting with an example and then explaining how it works after you’ve seen it in action. We’re going to start with the data from Subject 10 from the ERP CORE MMN experiment, which we looked at in the previous chapter in the section on eye movement detection.

First, reset everything with the Reset button so that any previous settings don't change how things look in this exercise. Then set Chapter_9 to be Matlab’s current folder using the Current Folder button in the EEGsets panel, and load the EEGset named 10_MMN_preprocessed.set. Make sure you’re getting the data from the Chapter_9 folder rather than the Chapter_8 folder, because the data are preprocessed a little differently but the filenames are the same.

The dataset you loaded has been high-pass filtered at 0.1 Hz and referenced to the average of P9 and P10 (except for the VEOG-bipolar and HEOG-bipolar channels). Scroll through the EEG to remind yourself what it looks like. You’ll see lots of blinks and eye movements and a fair amount of EMG noise.

When we looked at this participant’s data in the previous chapter, the EEG had been low-pass filtered. Without the low-pass filtering, you can now more easily see the spike potential at the beginning of each eye movement, which is a result of the muscle activity that sets the eyes in motion. This is illustrated in Figure 9.1, which shows the first 5 seconds of data from Subject 10. You can see a brief spike in several channels at the moment when the HEOG voltage starts to change from one level to another (which indicates a change in gaze location from one location to another). You may want to adjust the vertical scale and which channels are displayed so that you can see these voltage changes more clearly.

Now let's start the artifact correction process. Go to the EEGLAB ICA panel, click the Decompose data button, and set the parameters as shown in Screenshot 9.1. Everything is probably already set correctly by default, except for the channels, which you should set to 1:31. With this channel list, we’re including all the channels that have the same reference, and we’re excluding VEOG-bipolar and HEOG-bipolar. ICA assumes that the channels are linearly independent of each other, and the VEOG-bipolar and HEOG-bipolar were created by recombining other channels, so they would violate this assumption. However, they are useful channels to have later so that we can tell when blinks and eye movements occurred even after we’ve corrected the EEG data.

Click OK to start the ICA decomposition process, which figures out the independent components that best represent the current EEGset. You’ll see some text appear in the Matlab command window, followed by a long series of lines that look like this:

step 1 - lrate 0.001000, wchange 16.88431280, angledelta  0.0 deg

step 2 - lrate 0.001000, wchange 0.90426626, angledelta  0.0 deg

ICA works by training a machine learning algorithm (much like a neural network), and it can take the algorithm a very long time to converge on a solution. The amount of time depends on the number of channels, the noise level of the data, and the speed of your computer. It took over 5 minutes to reach a solution for this dataset on an old laptop, but only about 30 seconds on my current laptop (which has an Apple M4 Pro processor). If you don’t want to wait for it to finish, you can kill the process by clicking the Interrupt button that should be visible while ICA is running. If that doesn’t work, you can usually kill a Matlab process by typing Ctrl-C.

If you waited for the decomposition process to finish, you can name the resulting EEGset 10_MMN_preprocessed_ica. If you didn’t wait for the process to finish, you can just load the file 10_MMN_preprocessed_ica.set, which has the ICA weights that I got when I ran the decomposition.

ICA decomposes the EEG data into a set of independent components or ICs. These are statistically defined components, not ERP components, and they may not be biologically meaningful. I will refer to them as ICs rather than as components to maintain this important distinction.

Each IC is a spatial filter. When we apply an IC to the distribution of voltages over the scalp at a given moment in time, the output is the magnitude of that IC at that moment in time. Conveniently, ICA creates scalp maps of the ICs showing how each channel is weighted by the IC. Take a look at the scalp maps now by clicking Inspect/label ICs in the EEGLAB ICA panel. Tell the tiny window that appears that you want to plot components 1:31. You should then see a new window that looks like Screenshot 9.2. You'll also see a window that says While you are using the Inspect/Label ICs tool, the main ERPLAB Studio window will be frozen. Click "Okay" when you have closed the Inspect/Label ICs tool.

If you loaded the EEGset that already had the weights in it, the scalp maps you see should look exactly like those in the screenshot. However, if you ran the component decomposition process yourself, the maps might look a little different. This is because the ICA machine learning algorithm starts from a set of random weights and therefore doesn’t yield exactly the same results every time. If things are working well, then the results should be similar if you repeat the decomposition process multiple times. If you get very different results, then something is wrong (usually either very noisy data, not enough data, or a linear dependency among your channels).

There are some things that may differ across repetitions that don’t indicate a problem. First, the ordering of the ICs might change a little. For example, the scalp maps for ICs 13 and 14 might be swapped in your data. Second, there might be small changes in the weights for a given scalp map. Third, the polarity of one or more maps might be inverted (which we will discuss in more detail later). These are normal and no cause for concern. However, if the maps you obtained are radically different from those in Screenshot 9.2, then something is wrong. Perhaps you forgot to specify channels 1:31 in the window of parameters that appeared when you clicked Decompose data.

There are 31 ICs because we gave the algorithm data from 31 channels. Ordinarily, the number of ICs is equal to the number of channels. This is necessary to make the math work. The fact that the number of ICs equals the number of channels reflects an important way in which ICA is an imperfect method for identifying the true components underlying your data: Putting more or fewer electrodes on the scalp doesn’t change the number of sources of brain activity, and yet it changes the number of ICs. We will discuss this more later.

After the decomposition algorithm finishes, it reorders the components in terms of the amount of variance in the data they account for, with IC 1 accounting for the greatest variance. If the participant blinked frequently, the IC corresponding to the blinks is usually in the top five ICs (and is often IC 1). In our current example, IC 1 has the kind of frontally focused scalp distribution that you’d expect for a blink. But is this really a blink IC, or is it some other kind of artifact, or maybe even frontally distributed brain activity?

There are a few things we can check to determine whether IC 1 reflects blinks. If you click the number above IC 1, a new window will pop up showing you more details about this IC. In the upper right of the window, you’ll see a heat map showing how the magnitude of this IC varies over time. The X axis is time within a “trial” ("trials" are arbitrary time periods when continuous data are used). The Y axis is time over the course of the recording. It’s labeled Trials but actually corresponds to time in continuous data. You can see many small red blobs. Each of those is a brief time period in which the magnitude of IC 1 was large. This is exactly what we’d expect for blinks: many brief blips of voltage. The window also shows the frequency spectrum. The spectrum for IC 1 shows a gradual decline as the frequency increases, leveling out around 20 Hz, which is also characteristic of blinks (but is also seen for many types of non-blink activity). Close this window, then close the window with the IC scalp maps, and click Okay in the window that says While you are using the Inspect/Label ICs tool... to unfreeze the main eStudio window.

The next thing you can do (which I find the most informative) is to directly compare the time course of the ICs with the time course of the EEG and EOG signals. Start by going to the Plot Settings panel, checking the box for Display ICs, setting the Vertical Scale to 150 for the channels and 60 for the ICs, and clicking Apply. Then go to the Channel and IC Selection panel, select FP1, FP2, and all the EOG channels from the Channels column, select ICs 1-8 from the ICs column, and click Apply. The result should look like Screenshot 9.4. Now you can look at how the activations for ICs 1-8 vary over time and compare them to the key EOG and EEG signals at the corresponding time points. Try scrolling through the first 60 seconds of data to see whether IC 1 matches the time course of the blinks in the VEOG and frontal EOG channels.

The first thing to notice is that each blink that you can see in the EOG and EEG channels is accompanied by a deflection with the same shape in IC 1 (e.g., at 2.2 and 4.4 seconds in Screenshot 9.4). You should also notice that there are occasional step-like deflections in both VEOG-bipolar and IC 1 (e.g., between 16 and 17 seconds). These are vertical eye movements. Often, but not always, blinks and vertical eye movements will be captured by the same IC. Given the scalp distribution of IC 1 and its close correspondence in time with the blinks and eye movements in the VEOG-bipolar channel, we can be quite sure that this IC reflects vertical EOG activity.

Now let’s remove IC 1 from the EEG/EOG data to eliminate the blinks and vertical eye movements. Ordinarily, we would also remove other artifact-related ICs at this point (e.g., horizontal eye movements), but we’ll just focus on blinks and vertical eye movements for now. To do this, go to the EEGLAB ICA panel and click the Remove ICs button. A small window will appear; put 1 in the text box labeled List of component(s) to remove from data so that we remove activity corresponding to IC 1, and click OK. Name the resulting EEGset 10_MMN_preprocessed_ica_pruned and save it to disk with this name.

As shown in Screenshot 9.5, the blinks should now be eliminated from all channels except for the VEOG-bipolar channel (which we excluded from ICA when we specified channels 1:31 in the ICA decomposition step). Scroll through the first 60 seconds of data. You can still see ordinary EEG deflection and horizontal eye movements, but the blinks are largely gone. It's like a miracle!

You can also see that the vertical eye movements are largely gone as well. For example, if you go back to the previous EEGset, you'll see a sudden step in voltage at 5.7 seconds, but this step is gone after correction. 

This really is a minor miracle. We’ve nearly eliminated the contribution of blinks and vertical eye movements to the EEG at every single time point. Now we don't need to reject trials with blinks, which should increase the data quality substantially.

---