ICA rejection for Raw and Epochs

[dengemann]

I would like to work on including ICA rejection functionality based on the FAST ICA from scikits-learn. Who else is interested in that? In the course of the next days I would like to start doing some first drafts on the API. One question certainly is whether to 'sex-up' the objects or whether to provide a more module like API, e.g. mne.ica or mne.signal.ica. Another issue would be how to handle the visual checking and component-(de)selection. I guess the numpy plottting is too slow for that. One way however might be to write fiffs in a way that men_browse_raw could be used for that. Channel names then would be components, and bad channels then would be bad components. The data dispatch could certainly handled by either modifying the object's data in place (interactive python session) or via modified Raw / Epochs constructors for interoperation with mne_browse_raw.

Wdyt?

[agramfort]

thanks for looking into this

I would put ica "denoising" in mne.artifacts.ica that takes as input a Raw and returns a the ICA rotation matrix and (if requested) a new Raw that can be visuallized in mne_browse_raw.

you can visualize topography in python for rejection based on topo [1]. For time series you will need mne_browse_raw.

[1] https://github.com/agramfort/mne-python/tree/topo2d

[dengemann]

Thanks for the pointers, that sounds good to me.

On 18.09.2012, at 16:43, Alexandre Gramfort notifications@github.com wrote:

thanks for looking into this

I would put ica "denoising" in mne.artifacts.ica that takes as input a Raw and returns a the ICA rotation matrix and (if requested) a new Raw that can be visuallized in mne_browse_raw.

you can visualize topography in python for rejection based on topo [1]. For time series you will need mne_browse_raw.

[1] https://github.com/agramfort/mne-python/tree/topo2d — Reply to this email directly or view it on GitHub.

[dengemann]

'you can visualize topography in python for rejection based on topo' ... wouldn't that be somewhat misleading? The resulting components aren't channels. Ok, but of course I see where your suggestion is going. One could just create some new kind of layout which prevents misleading interpretations.

[agramfort]

'you can visualize topography in python for rejection based on topo' ... wouldn't that be somewhat misleading? The resulting components aren't channels. Ok, but of course I see where your suggestion is going. One could just create some new kind of layout which prevents misleading interpretations.

I had something like this in mind:

http://lh4.ggpht.com/_BfUxGfm9eUY/S2yP2Ej3KSI/AAAAAAAAhtI/M9adaGj18VI/2DOR.jpg

(the way EEGLAB works)

for Neuromag that combines 3 types of sensors it's less obvious what's the way to go but showing one sensor type seems like a good start.

A

[dengemann]

Ok, got you, i like it!

As it seems to me there are several ways to use sklearn.decomposition.FastICA among those a clearly oo variant. Looking at the mne code however there's a clear preference for introducong classes at a late stage of development although for the ica workflow a -- let's say -- raw_ica object might be very useful. The other option would be to distribute the workflow across several fuctions for decomposing, investigating, deselecting and finally mixing. This would certainly bring the advantage of maximum flexibility for the user. On the other hand side it means more code and more parameter passing. For the sake of clarity and consistency I would at least for now opt for the rather fuctional variant. Wdyt?

On 23.09.2012, at 15:55, Alexandre Gramfort notifications@github.com wrote:

'you can visualize topography in python for rejection based on topo' ... wouldn't that be somewhat misleading? The resulting components aren't channels. Ok, but of course I see where your suggestion is going. One could just create some new kind of layout which prevents misleading interpretations.

I had something like this in mind:

http://lh4.ggpht.com/_BfUxGfm9eUY/S2yP2Ej3KSI/AAAAAAAAhtI/M9adaGj18VI/2DOR.jpg

(the way EEGLAB works)

for Neuromag that combines 3 types of sensors it's less obvious what's the way to go but showing one sensor type seems like a good start.

A — Reply to this email directly or view it on GitHub.

[agramfort]

As it seems to me there are several ways to use sklearn.decomposition.FastICA among those a clearly oo variant.

you have both the estimator object of just the fastica function. The object gives you access to the mixing matrix with get_mixing_matrix

Looking at the mne code however there's a clear preference for introducong classes at a late stage of development although for the ica workflow a -- let's say -- raw_ica object might be very useful. The other option would be to distribute the workflow across several fuctions for decomposing, investigating, deselecting and finally mixing. This would certainly bring the advantage of maximum flexibility for the user. On the other hand side it means more code and more parameter passing. For the sake of clarity and consistency I would at least for now opt for the rather fuctional variant.

can you write a code snippet going over the use case, the way you would like to use it without thinking about the internals of ICA?

[dengemann]

Sure, I already started doing some usage examples. I will get back to you as soon as I'm happy with it. Shouldn't take long anymore.

[dengemann]

Ok, just to keep it simple and stupid here a first implementation basically drawing on the sklearn examples

(it works)

from sklearn.decomposition import FastICA
import numpy as np

def decompose_raw(raw, picks, n_components=None, start=None, stop=None,
                  *args, **kwargs):
    """ Run ICA decomposition on instance of Raw
    Paramerters
    -----------
    raw : instance of Raw
        MEG raw data
    picks : array-like
        indices for channel selection as returned from mne.fiff.pick_channels
    n_components : integer
        number of components to extract. If none, no dimension
        reduciton will be applied.
    start : integer
        starting time slice
    stop : integer
        final time slice

    Return Values
    -------------
    latentent_sources : ndarray
        time slices x n_components array
    mix_matrix : ndarray
        n_components x n_components
    """
    start = raw.first_samp if start == None else start
    stop = raw.last_samp if stop == None else stop
    data, _ = raw[picks, start:stop]
    data = data.T  # sklearn expects column vectors / matrixes

    print ('\nCalculating signal decomposition.'
           '\n    Please be patient. This may take some time')
    ica = FastICA(n_components=n_components, *args, **kwargs)
    latent_sources = ica.fit(data).transform(data)
    mix_matrix = ica.get_mixing_matrix()
    assert np.allclose(data, np.dot(latent_sources, mix_matrix.T))

    print '\nDone.'
    return latent_sources, mix_matrix

[dengemann]

And here the very first usage example.

Wdyt?

import mne
from mne.fiff import Raw
from mne.artifacts.ica import decompose_raw

from mne.datasets import sample
data_path = sample.data_path()
raw_fname = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw.fif'

raw = Raw(raw_fname)

picks = mne.fiff.pick_types(raw.info, meg=True, eeg=False, stim=False)

start, stop = raw.time_to_index(100, 115)

sources, mix = decompose_raw(raw, picks, start=start, stop=stop)

# then do something with that...

[agramfort]

besides the problem of handling the whitening process it seems like a very natural approach.

we could have:

sources, mix = ica_decompose_raw(raw, cov, picks, start=start, stop=stop)
# a process to select the sources / components we want to keep
# and then recompose
raw_denoised = ica_recompose_raw(raw, cov, mix, bads=[1, ...], ...)

where bads contains the index of the bad components

[dengemann]

Makes perfect sense to me. What would the cov specifically do here? Wouldn't sources be sufficient? The whitening process at least can be controlled via the args and kwargs passed to the FastICA inside ica_decompose_raw. For the example case also dimensionality reduction makes sense, with n_components = 25 it goes really fast. Another thing to be added is the projection handling. On / off? ICA before / after SSP and so on.

As to the 'process to select the sources / components we want to keep'. Besides visual inspection one could also think about a similarity metric for automatically rejecting components similar to the ECG / EOG channel. In addition to the fieldmap plots one could also think about a trellis plot. I was thinking of something like:

http://slendrmeans.wordpress.com/2012/04/27/will-it-python-machine-learning-for-hackers-chapter-1-part-5-trellis-graphs/

http://slendrmeans.files.wordpress.com/2012/04/ufo_ts_bystate.png

For all the components. Wdyt?

[agramfort]

Makes perfect sense to me. What would the cov specifically do here? Wouldn't sources be sufficient?

you need a noise covariance to whiten the data before calling ICA when you have different types of sensors.

The whitening process at least can be controlled via the args and kwargs passed to the FastICA inside ica_decompose_raw. For the example case also dimensionality reduction makes sense, with n_components = 25 it goes really fast. Another thing to be added is the projection handling. On / off? ICA before / after SSP and so on.

FastICA won't know about the type of sensors.

As to the 'process to select the sources / components we want to keep'. Besides visual inspection one could also think about a similarity metric for automatically rejecting components similar to the ECG / EOG channel.

indeed.

In addition to the fieldmap plots one could also think about a trellis plot. I was thinking of something like:

http://slendrmeans.wordpress.com/2012/04/27/will-it-python-machine-learning-for-hackers-chapter-1-part-5-trellis-graphs/

http://slendrmeans.files.wordpress.com/2012/04/ufo_ts_bystate.png

For all the components. Wdyt?

looks good.

[dengemann]

So I would do something like

# ... module

from ..cov import compute_whitener

# ... e.g. ica_decomompose

whitener, _ = compute_whitener(noise_cov, info, picks)

And then pass the already whitened data to the ICA with whiten=False?

Also I would allow to order the components by variance or kurtosis. This could be accessed via a dict similar to reject in Epochs.

[agramfort]

from ..cov import compute_whitener

... e.g. ica_decomompose

whitener, _ = compute_whitener(noise_cov, info, picks)

yes exactly. Once you apply it to the data then you can feed everyting to PCA-ICA

Also I would allow to order the components by variance or curtosis. This could be accessed via a dict similar to reject in Epochs.

yes :)

[dengemann]

you need a noise covariance to whiten the data before calling ICA when you have different types of sensors.

As we're talking about raw data decomposition event related covs wouldn't make sense, right? Is there an empty room data set available for calculating a noise covariance matrix? I looked into the sample data but could not find anything looking like that.

[agramfort]

point taken. Yes empty room would work. In case you don't have it and you have different sensor types, what I suggest is to z-score each channels (equivalent of whitening with a diagonal covariance) or z-score blocks of channels of the same type. You basically want to make EEG around 1e-6 comparable with MEG around 1e-13. ok?

[dengemann]

2012/9/25 Alexandre Gramfort notifications@github.com

point taken. Yes empty room would work. In case you don't have it and you have different sensor types, what I suggest is to z-score each channels (equivalent of whitening with a diagonal covariance) or z-score blocks of channels of the same type. You basically want to make EEG around 1e-6 comparable with MEG around 1e-13. ok?

I mean, we are trying to make a nice example using the sample data, right? Of course this should work. If there was an empty-room measuremnt it would be nicer though ;-) But right, this way we will create some options for different situations. For some data sets there simply might not be an empty room... To make this fit into the current API sketch, we could say, if noise_cov == None z-scores for channel groups will be applied. Wdyt?

—

Reply to this email directly or view it on GitHubhttps://github.com/mne-tools/mne-python/issues/97#issuecomment-8866313.

[agramfort]

+1 I like that

[dengemann]

Ok, just to make sure we are on the right track: i now wrote some lines for separately z-scoring the chan types.

I then would pass the resulting matrix to ica with whiten=True (with noise_cov != None whiten would be false)

On 25.09.2012, at 21:19, Alexandre Gramfort notifications@github.com wrote:

+1 I like that

— Reply to this email directly or view it on GitHub.

[agramfort]

always use whiten=True for ICA even after using the noise_cov. Besides this we are on the same page

[dengemann]

So here's som updates:

• selective z-scoreing in absence of noise_cov works (but I'm not sure whether it does what we want) and the ICA runs through with all channel types included
• however the assert command trying to state the identity between data and re-composed data fails if i reduce the dimensionality via ICA/PCA. This, however, only after the z-scoring.
• interestingly the assert fails without z-scoring and EEG data included

I'm somewhat stuck here, any pointers / suggestions?

See:

https://github.com/dengemann/mne-python/blob/ica/mne/artifacts/ica.py https://github.com/dengemann/mne-python/blob/ica/examples/artifacts/compute_ica_components.py

[agramfort]

Can you open a WIP pull request to see the diff more easily?

[dengemann]

Closing this to continue on the related WIP PR.