MRG+1: ENH: functionality for doing ICA on MEG reference channels

[jshanna100]

I recently had some success removing a pattern of intermittent external noise using ICA decompositions of the reference channel array in a 4D MEG system. The procedure went as follows:

1. Use ICA to decompose the MEG reference channels into components. I used n_components=0.99, which generally yielded one or two components.
2. Decompose the MEG data channels into components.
3. Use preprocessing.ICA.scores_sources to get pearsonr for each individual reference component against the data channel components.
4. Eliminate data channel components that have a z-transformed pearsonr with reference component beyond a threshold.

MNE as it currently stands does not allow ICA decomposition of reference channels. The suggested changes are a first step in providing this functionality. They allow the decomposition itself, as well as permitting visualisation functions like e.g. preprocessing.ICA.plot_properties. The topographies are nonsense, but plotting functions at least now do not fail entirely, allowing inspection of other types of visualisations.

One possible substitute for the topographies might be quiver plots, where the arrows orientations are determined by the orientations of the reference channels, projected onto 2D space. I could also put together a proper component identifier function that carries out steps 3 and 4.

[larsoner]

Naive question, why do all of this rather than simply doing a joint decomposition of the reference and non-reference channels together, and reject components that are most clearly in the reference channels?

I know mechanically in MNE master you can't do this because we don't allow ref channels to be in the decomposition. But it looks like this branch would allow it. Have you tried this simpler approach and it does not work?

Admittedly I don't know why we don't allow ref channels in ICA decompositions in the first place...

[jshanna100]

I've just tried this approach, and it also works well. The one advantage to doing it the original way is having a unified algorithmic approach: you're essentially doing the same thing as with automatic EOG- or EKG-related component removal; the only difference is instead of using the EOG or EKG channels as your ground truth, you're using the reference channel components.

With the approach you've suggested, you'll ultimately need to also decide an algorithm which determines "most clearly in the reference channels," which probably wouldn't be so challenging, but it would also be different than the established tools already built into MNE.

[nbara]

Apologies if this is slightly OT, but this reminds me of the following two papers :

1. De Cheveigné, A., & Simon, J. Z. (2007). Denoising based on time-shift PCA. Journal of neuroscience methods, 165(2), 297-305. https://www.sciencedirect.com/science/article/pii/S0165027007002701
2. De Cheveigné, A., & Simon, J. Z. (2008). Sensor noise suppression. Journal of neuroscience methods, 168(1), 195-202.
   https://www.sciencedirect.com/science/article/pii/S0165027007004621

If this seems relevant, I coded those in python a while back (but then stopped working with MEG so I never took the time to fully test it on real data).

[codecov]

Codecov Report

Merging #5807 into master will decrease coverage by 0.01%.
The diff coverage is 72.22%.

@@            Coverage Diff             @@
##           master    #5807      +/-   ##
==========================================
- Coverage   88.69%   88.67%   -0.02%     
==========================================
  Files         396      396              
  Lines       71870    71905      +35     
  Branches    12071    12081      +10     
==========================================
+ Hits        63744    63762      +18     
- Misses       5198     5213      +15     
- Partials     2928     2930       +2

[jshanna100]

It should be working now. 1) ICA fitting accepts reference channels 2) any attempt to visualise a topography of components derived only from reference channels will simply return and:
3) there is now a new function in class ICA, find_bads_ref. This is directly adapted from analogous functions find_bads_eog and find_bads_ecg. The only difference is that it searches for channels starting with REF_ICA, instead of EOG/EKG channels. These of course must first be appended to the inst object which is passed to the function.

I've found that this approach is good not only for removing intermittent environmental noise, but to my astonishment, EKG components are in nearly every case detected, and once or twice I've even picked up EOG components.

[larsoner]

@jshanna100 can you look into why Travis is unhappy? It's better to review after CIs are green. But what you say sounds useful!

[jshanna100]

The Travis problems seem to come from the fact that the ICA tests use a threshold of .9 of the variance. Now that reference channels are allowed, it attempts an ICA on those, and can find no more than one component under that threshold, which results in an error.

[larsoner]

I see you have not modified test_ica.py, so this suggests that you have made some backward-incompatible changes. Looking at the code, ref_meg is always included now (right?), meaning ICA done on 0.17 or earlier will not be the same on 0.18 or later. We have to be careful with such changes. And provide a gradual upgrade path for people.

How about:

1. adding a ref_meg=False (by default for backward compat) kwarg to the ICA constructor
2. add a test with ref_meg=True for the dataset that currently fails, wrapped in a with pytest.raises(...): context manager
3. add a test with ref_meg=True with a higher threshold (or integer component number) so you get a more reasonable component count

[larsoner]

How much of this code is shared with find_bads_eog / find_bads_ecg? If it's 95%+, then we should refactor so that this and those functions each just call a return self._find_bads_type(...) to reduce code duplication.

[jshanna100]

It's maybe about 75% between eog and ref. ecg is fairly different. I've left them as is, for now.

[agramfort]

@jshanna100 you should then refactor to have a private function eg _find_bads_ch_type that takes the channel type you want to correlate with and make find_bads_ecg and find_bads_ref_meg use this private function. This would avoid the code duplication.

also you can run:

make pep

from the mne-python folder. It will show you the code formatting errors that prevent travisci to be green.

thanks a lot for your efforts.

[jshanna100]

Ok, this is done, I think. It works on some of my own data with reference channels. I've run some quick tests on the MNE sample data, and as far as I can tell, the ECG and EOG functions also seem to work as they should.

[jshanna100]

Sorry, "make pep," passed on my system, but now the ica_test.py is broken. Will sort out. Please bear with me; this is my first time sorting out a PR here...

[larsoner]

CircleCI:

/home/circleci/project/mne/preprocessing/ica.py:docstring of mne.preprocessing.ICA.find_bads_eog:63: WARNING: py:obj reference target not found: _find_bads_ch
/home/circleci/project/mne/preprocessing/ica.py:docstring of mne.preprocessing.ICA.find_bads_ref:63: WARNING: py:obj reference target not found: _find_bads_ch

You can't reference private functions in public docstrings (users should not use or think about them)

[larsoner]

The refactoring looks nice!

[jshanna100]

Alright, good! Got a deadline this week, so will get to all this first half of the next one.

[agramfort]

ok cool

[jshanna100]

I'm having a hard time figuring out what the problem with travis is now. All I changed from last commit was some commented lines. Does this have to do with something else, somehow?

[larsoner]

Just ensure sure you did not add errors to the files you've edited, #5892 will fix the flake problems with the others.

[larsoner]

If you rebase, it should fix the CircleCI build problem (I think)

[jshanna100]

Alright, that appears to be sorted now.

I don't know how to interpret these codecov errors. Are they a problem?

[larsoner]

codecov can be a bit unreliable sometimes, it's probably okay. I see a few more minor things to change, okay for me to push a commit to your branch?

[jshanna100]

Yep, go for it.

[jshanna100]

Is there anything left to be done on this?

[larsoner]

LGTM +1 for merge

@agramfort feel free to merge if you're happy

[agramfort]

thx @jshanna100