Cannot find ecg artifacts during ica

[YuZhou]

• MNE-Python version: 0.22.0
• operating system: MacOS Mojave version 10.14

Hey guys,
When I use find_bads_ecg() to find which ICs match the ECG pattern, it turns out that the ‘ecg_indices’ was empty no matter I modify the number of independent components as 15, 20, or 30. However, there did exist a component (ICA002) that looked like ecg in the ica components plot produced by ‘ica.plot_sources(raw_notch_filt_cp)’.

But the scalp polt of ICA002 didn’t look like the ecg artifacts should.

And when I create ecg epochs:

ecg_evoked = create_ecg_epochs(raw_notch_filt, reject_by_annotation=True).average()
ecg_evoked.apply_baseline(baseline=(None, -0.2))
fig_raw_evoked_ecg = ecg_evoked.plot_joint()

the output messages indicated as below, which suggested that the heartbeats were found:

Reconstructing ECG signal from Magnetometers
Setting up band-pass filter from 8 - 16 Hz

FIR filter parameters
---------------------
Designing a two-pass forward and reverse, zero-phase, non-causal bandpass filter:
- Windowed frequency-domain design (firwin2) method
- Hann window
- Lower passband edge: 8.00
- Lower transition bandwidth: 0.50 Hz (-12 dB cutoff frequency: 7.75 Hz)
- Upper passband edge: 16.00 Hz
- Upper transition bandwidth: 0.50 Hz (-12 dB cutoff frequency: 16.25 Hz)
- Filter length: 10000 samples (10.000 sec)

Number of ECG events detected : 603 (average pulse 90 / min.)
Not setting metadata
Not setting metadata
603 matching events found
No baseline correction applied
Loading data for 603 events and 1001 original time points ...
36 bad epochs dropped
Applying baseline correction (mode: mean)
No projector specified for this dataset. Please consider the method self.add_proj.

But the plots looked strange:

I really have no idea which part went wrong , is there anything I can do to fix this?

Thanks a lot!

Best

[agramfort (Alexandre Gramfort)]

yes you can guess ECG in ICA002 but the unmixing is clearly not good enough.
Your ICA002 is still dominated by other artifacts. Maybe try to highpass the data a bit more
before doing ICA. Or just mark manually ICA002 as bad as the topo suggests it’s really ECG

Alex

[richard (Richard Höchenberger)]

I’d like to second what @agramfort said. Can you show how you create your ICA object?

[YuZhou]

Hi Richard,

Thanks again for your prompt reply
Here is the code I use:

# concatenate 4 perception runs
raw_perception = mne.concatenate_raws([raw_per1, raw_per2, raw_per3, raw_per4])
meg_picks = mne.pick_types(raw_perception.info, meg=True)
freqs = (50, 100, 150, 200, 250, 300)
raw_notch_filt = raw_perception.copy().load_data().notch_filter(freqs=freqs, picks=meg_picks, method='spectrum_fit', filter_length='10s')
annot_muscle, scores_muscle = annotate_muscle_zscore(raw_notch_filt, ch_type="mag", threshold=4, min_length_good=0.1, filter_freq=[110, 140])
raw_notch_filt.set_annotations(annot_muscle)

# the ica part
raw_notch_filt_hp = raw_notch_filt.copy()
raw_notch_filt_hp.load_data().filter(l_freq=1., h_freq=None)
ica = ICA(n_components=30, random_state=98, method='picard')
ica.fit(raw_notch_filt_hp, reject_by_annotation=True)

Actually, I did the highpass filter with l_freq=1, should I set a higher threshold?

[richard (Richard Höchenberger)]

The 1 Hz high-pass filter is a good choice. I’m wondering if you’re trying to squeeze to much info into too few components. Try n_components=0.9 or n_components=0.8. How many components will this produce?

[YuZhou]

Hi Alex,

Thanks for your suggestions very much! Because I have many subjects to run, and I hope the processing workflow to be as automatic as possible, so it’s better to recognize the ecg and eog by the script itself, that’s why I choose the way of ’ Using a simulated channel to select ICA components’ in the tutorial

[YuZhou]

Hi,
first I tried n_components=0.9, it produced 24 components; then I tried n_components=0.8, it produced 14 components. And neither of these two options found the ecg component

[richard (Richard Höchenberger)]

@YuZhou

first I tried n_components=0.9, it produced 24 components; then I tried n_components=0.8, it produced 14 components.

This seems odd, how many sensors are there? I thought you were working with MEG data?

[YuZhou]

Yes, it’s MEG data, 306 sensors in total (mag + grad)

[richard (Richard Höchenberger)]

@agramfort What do you think of that? Why is so much variance explained by so few components?

@YuZhou Can you try again with n_components=0.999? How many components do you get?

[YuZhou]

When n_components=0.999, it shows that ‘Selecting by explained variance: 66 components’, and there was a warning saying

/anaconda3/envs/mne/lib/python3.7/site-packages/picard/solver.py:216: UserWarning: Picard did not converge. Final gradient norm : 3.031e-05. Requested tolerance : 1e-07. Consider increasing the number of iterations or the tolerance.
  % (gradient_norm, tol))

And there still no ecg component being found

# find which ICs match the ECG pattern
ecg_indices, ecg_scores = ica.find_bads_ecg(raw_notch_filt, reject_by_annotation=True,  
                                            method='ctps', threshold='auto')
ica.exclude = ecg_indices

[richard (Richard Höchenberger)]

@YuZhou

When n_components=0.999, it shows that ‘ Selecting by explained variance: 66 components ’, and there was a warning saying

/anaconda3/envs/mne/lib/python3.7/site-packages/picard/solver.py:216: UserWarning: Picard did not converge. Final gradient norm : 3.031e-05. Requested tolerance : 1e-07. Consider increasing the number of iterations or the tolerance.
  % (gradient_norm, tol))

The number of components sounds much more sensible to me now, so let’s stick with this value for n_components for now.

The warning you see already suggests a solution When creating your ICA object, set max_iter to a higher value, e.g. max_iter=1000. The warning should hopefully disappear then.

[YuZhou]

You are right When n_components=0.999 and max_iter=1000, the warning disappeared. But sadly, although I can see that the ICA008 is heartbeats:

the scalp map of ICA008 does still not look like the ecg should:

and still, no ecg component has been found

[richard (Richard Höchenberger)]

@YuZhou

ecg_indices, ecg_scores = ica.find_bads_ecg(raw_notch_filt, reject_by_annotation=True,  
                                            method='ctps', threshold='auto')

Whoops. Can you try to run find_bads_ecg on your ECG epochs, not the raw data? Not that it should really make a difference, but … might be worth a try?

Also can you try method='correlation' to see if this gives any better results?

[YuZhou]

First try:

ecg_epochs = create_ecg_epochs(raw_notch_filt, reject_by_annotation=True)
ica.exclude = []
# find which ICs match the ECG pattern
ecg_indices, ecg_scores = ica.find_bads_ecg(ecg_epochs, reject_by_annotation=True,  
                                            method='ctps', threshold='auto')
ica.exclude = ecg_indices

Returns:

Reconstructing ECG signal from Magnetometers
Using threshold: 0.16 for CTPS ECG detection

and ‘ecg_indices’ returns []

Second try:

ecg_indices, ecg_scores = ica.find_bads_ecg(ecg_epochs, reject_by_annotation=True,  
                                            method='correlation', threshold='auto')
ica.exclude = ecg_indices

Returns:
Reconstructing ECG signal from Magnetometers
and happily, it finally gave me an index 36! But……it’s not 8, which is a little bit odd. Below is what ICA036 looks like and the scores plot:
the raw plot of ICA036 looks less like heartbeats than ICA008:

but the scalp map of ICA036 looks more like heartbeats than ICA008:

it’s really hard to decide whether to adopt the ecg components the MNE gave me, can you decide for me

[richard (Richard Höchenberger)]

This doesn’t look like ECG to me at all…
Have you tried not using the notch filter before ICA? Maybe this is causing some crazy things to happen? (It really shouldn’t, but who knows…?)

[YuZhou]

I tried to use the raw data only with annotations of muscle artifacts, not using notch filtering, the ecg component couldn’t be found again.
I’m considering that whether it’s the problem of the original heartbeats because when I draw the heartbeats evoked plots, it really seems not like heartbeat waves as shown in the tutorial:

Here is the drawing code:

ecg_evoked = create_ecg_epochs(raw_perception, reject_by_annotation=True).average()
ecg_evoked.apply_baseline(baseline=(None, -0.2))
fig_raw_evoked_ecg = ecg_evoked.plot_joint()

[richard (Richard Höchenberger)]

I would suggest you try to extract the heartbeat events via find_ecg_events and adjust the qrs_threshold parameter and see what that gets you… you can then use those events to create ECG epochs manually, which you can average and visualize.

[YuZhou]

Hi Richard,
Thanks a lot for your patient help .
I tried the way you suggested, but no matter how I tried the different options for qrs_threshold, the ecg plots just didn’t get much better. I almost want to give up the way of automatically finding the ecg and eog, and turn to manually select the components to be reject.
Thank you again!

[richard (Richard Höchenberger)]

@YuZhou

but no matter how I tried the different options for qrs_threshold, the ecg plots just didn’t get much better.

Since you don’t have a proper ECG channel in your data, find_ecg_events() synthesizes one. You can ask the function to return the synthesized signal by passing return_ecg=True (see the docs) I’d do just that and plot the synthesized ECG signal to get an idea whether it even makes sense?