Applying ICA with No Components Selected out but Signal Changes Using Default Parameters

Hi,

    In an analysis, I am running:

ica = ICA(method='fastica', n_components=n_components, # n_components=None
          random_state=seed)
ica.fit(inst2)
...
inst2 = ica.apply(inst2, exclude=ica.exclude)

    and when I skip all intermediate steps and just fit the ICA and apply it with an empty list for ica.exclude the signal still changes, quite a bit. I thought if no components were selected out and all the max PCA components were used the signal would be unchanged or basically unchanged. Is this a bug or something with my implementation?

Thanks,

Alex

Translational NeuroEngineering Laboratory
Division of Neurotherapeutics, Department of Psychiatry
Massachusetts General Hospital, Martinos Center
149 13th St Charlestown #2301, Boston, MA 02129
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.nmr.mgh.harvard.edu/pipermail/mne_analysis/attachments/20190702/f627fb36/attachment.html

Also, of note the ica scale is off by quite a lot in the plot_sources plot, it is way too zoomed in.

Alex

Translational NeuroEngineering Laboratory
Division of Neurotherapeutics, Department of Psychiatry
Massachusetts General Hospital, Martinos Center
149 13th St Charlestown #2301, Boston, MA 02129

External Email - Use Caution

Hi Alex,

Could you provide us a full script on the MNE sample data that we can run?

Mainak

External Email - Use Caution

hi,

can you check that the number of components you fit is equal to the number
of channels?
If it's less you have a dimensionality reduction step.

Alex

Hi Alex and Mainak,

    The n_components argument is given None which yields 57 components and there are 64 channels with 7 bad channels which are not included so no I don't think it's because of the dimensionality reduction. Maybe it's some whitening.

To see something similar to what I'm looking at as far as scaling you can use the script below but I haven't been able to replicate the changes after ICA with sample data. I filedropped you both test epochs to the emails you responded to the thread with that does show that.

from time import time
import matplotlib.pyplot as plt
import mne
from mne.preprocessing import ICA
from mne.datasets import sample

'''data_path = sample.data_path()
raw_fname = data_path + '/MEG/sample/sample_audvis-raw.fif'

raw = mne.io.Raw(raw_fname, preload=True)
events = mne.find_events(raw)
epochs = mne.Epochs(raw, events, preload=True)

'''
epochs = mne.read_epochs('test-epo.fif', preload=True)

epochs = epochs.pick_types(meg=False, eeg=True)

fig, (ax0, ax1) = plt.subplots(1,2)
epochs.average().plot(axes=ax0, show=False)

ica = ICA(method='fastica', random_state=0)
t0 = time()
ica.fit(epochs)
fit_time = time() - t0
epochs = ica.apply(epochs, exclude=ica.exclude)
epochs.average().plot(axes=ax1, show=False)
ica.plot_sources(epochs)

Thanks,

Alex

Translational NeuroEngineering Laboratory
Division of Neurotherapeutics, Department of Psychiatry
Massachusetts General Hospital, Martinos Center
149 13th St Charlestown #2301, Boston, MA 02129

External Email - Use Caution

Hi Alex,

I could replicate the pb and start to investigate.

You have a strong transient artifact in your data (stim artifact I suspect)
and this affects the conditioning of the mixing matrix. So when computing
the pinv here: https://github.com/mne-tools/mne-python/blob/master/mne/preprocessing/ica.py#L677
the mixing is not numerically the inverse of the unmixing matrix.

I don't know exactly how to fix this but that's a starting point.
He someone can look into this it's great.

Alex

External Email - Use Caution

I can confirm what Alex said. Here?s a simpler script to play with:

import mnefrom mne.preprocessing import ICA
from numpy.testing import assert_array_equal

epochs = mne.read_epochs('test-epo.fif', preload=True)

epochs.pick_types(meg=False, eeg=True)
epochs.crop(0.2, None)

epochs2 = epochs.copy()

reject = dict(eeg=200e-6)
ica = ICA(method='fastica', random_state=42, max_iter=500, n_components=0.95)
ica.fit(epochs.copy().drop_bad(reject=reject), decim=10)

ica.apply(epochs2, exclude=[])

epochs.average().plot()
epochs2.average().plot()
# from mne import combine_evoked# evoked_diff =
combine_evoked([epochs.average(), -epochs2.average()], 'equal')#
evoked_diff.plot()

assert_array_equal(epochs.get_data(), epochs2.get_data())

Using n_components=0.95 vs n_components=1.0 actually seems to make a huge
difference, perhaps for the reasons that Alex outlines. Ultimately, what
you reconstruct depends on the quality of your decomposition -- and it's
not perfect.

Best,
Mainak

Thanks Alex and Mainak,

    I?ll try changing the number of components and interpolating a larger section of the data to abolish the stim artifact (the data I sent was TMS-EEG data).

Alex

Translational NeuroEngineering Laboratory
Division of Neurotherapeutics, Department of Psychiatry
Massachusetts General Hospital, Martinos Center
149 13th St Charlestown #2301, Boston, MA 02129