Unexpected behaviour in find_bads_eog

[eduardo (Eduardo Perez-Valero)]

I have applied the following workflow for automated artifact rejection (based on this example).

# Autoreject (local) to clean epochs before ICA
ar_pre_ica    = AutoReject(random_state = seed).fit(epochs[:20])
_, reject_log = ar_pre_ica.transform(epochs, return_log = True)

# Fit ICA on non-artifactual epochs 
ica = mne.preprocessing.ICA(random_state = seed)
ica.fit(epochs[~reject_log.bad_epochs])

# Exclude blink artifact components (use Fp1 as EOG proxy)
epochs_eog  = mne.preprocessing.create_eog_epochs(raw = raw, ch_name = 'Fp1') 
eog_idx, eog_scores = ica.find_bads_eog(epochs_eog, ch_name = 'Fp1', verbose = 2)
ica.exclude = eog_idx
            
# Apply ICA
epochs_clean = ica.apply(epochs.copy())
            
# Autoreject (local) on blink-artifact-free epochs
ar_post_ica  = AutoReject(random_state = seed).fit(epochs_clean[:20])
epochs_clean = ar_post_ica.transform(epochs_clean)

The thing is, find_bads_eog does not seem to catch the blink artifact components despite they are quite evident. Actually, the eog_scores of the two first ICA components are 0.81 and -0.77, however I get eog_idx = []

Below is the plot of the ICA components.

I have also tried using n_components = 0.99 to reduce the number of components. In this case, the two blink components are combined into a single component with an eog_score = 0.97, however, find_bads_eog still does not catch this component as eog_idx = [].

Does anyone know why this happens? Thanks.

[agramfort (Alexandre Gramfort)]

can you share the data file with a script to replicate what you report?

Alex

[eduardo (Eduardo Perez-Valero)]

Hi Alexandre.

I’ve created a folder with the necessary data. You just need to download the folder here and run the script inside.

Thank you.
Eduardo

[richard (Richard Höchenberger)]

Hello @eduardo,

simply lowering the detection threshold seems to do the trick for me.

eog_idx, eog_scores = ica.find_bads_eog(
    epochs_eog,
    ch_name='Fp1',
    threshold=2  # default is 3
)
ica.plot_scores(scores=eog_scores, exclude=eog_idx)

Btw there’s an undocumented feature of find_bads_eog: the paramter ch_name can actually also be a list of channels, so you could do:

eog_idx, eog_scores = ica.find_bads_eog(
    epochs_eog,
    ch_name=['Fp1', 'Fp2'],
    threshold=2
)
ica.plot_scores(scores=eog_scores, exclude=eog_idx)

and MNE would do the scoring based on those two channels:

… better don’t rely on this as, like I said, it’s undocumented; but in some situations it might come in handy. (In this specific example, though, both channels yield roughly the same scores, so there’s not much to gain by using both…)

PS:
You’re using ica.plot_components() in your example script. If you pass a Raw or Epochs instance, you can make the figure interactive: ica.plot_components(inst=epochs[~reject_log.bad_epochs]) Clicking on the topomaps will then open the properties plot of the respective component.

[richard (Richard Höchenberger)]

Oh, and another thing I’ve noticed: the ICA doesn’t converge on my computer with the default max_iter setting (1000 for fastica). So I had to increase the value to reach convergence:

ica = mne.preprocessing.ICA(random_state=seed, max_iter=5000)

You may wish to consider using Picard, which typically converges much faster than fastica. To use it, pass method='picard' to ICA().

[eduardo (Eduardo Perez-Valero)]

HI @richard, when I lower the threshold, I don’t seem to get the same components you get for exclusion. I modified the ICA subsection following your recommendations:

# Fit ICA on non-artifactual epochs 
ica = mne.preprocessing.ICA(method = 'picard', random_state = seed)
ica.fit(epochs[~reject_log.bad_epochs])

# Exclude blink artifact components (use Fp1 as EOG proxy)
epochs_eog  = mne.preprocessing.create_eog_epochs(raw = raw, ch_name = 'Fp1') 
eog_idx, eog_scores = ica.find_bads_eog(epochs_eog, ch_name = 'Fp1', threshold = 2)
ica.exclude = eog_idx
ica.plot_scores(eog_scores)

Now ICA converges, but find_bads_eog excludes three components.